heavy_tails.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7a44c6ce",
   "metadata": {},
   "source": [
    "\n",
    "<a id='heavy-tail'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fa58672",
   "metadata": {},
   "source": [
    "# Heavy-Tailed Distributions\n",
    "\n",
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "975d6dae",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade yfinance pandas_datareader"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98f57f91",
   "metadata": {},
   "source": [
    "We use the following imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65a20dc1",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import yfinance as yf\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "\n",
    "from pandas_datareader import wb\n",
    "from scipy.stats import norm, cauchy\n",
    "from pandas.plotting import register_matplotlib_converters\n",
    "register_matplotlib_converters()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de08c0ad",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Heavy-tailed distributions are a class of distributions that generate “extreme” outcomes.\n",
    "\n",
    "In the natural sciences (and in more traditional economics courses), heavy-tailed distributions are seen as quite exotic and non-standard.\n",
    "\n",
    "However, it turns out that heavy-tailed distributions play a crucial role in economics.\n",
    "\n",
    "In fact many – if not most – of the important distributions in economics are heavy-tailed.\n",
    "\n",
    "In this lecture we explain what heavy tails are and why they are – or at least\n",
    "why they should be – central to economic analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd0eb776",
   "metadata": {},
   "source": [
    "### Introduction: light tails\n",
    "\n",
    "Most [commonly used probability distributions](https://intro.quantecon.org/prob_dist.html) in classical statistics and\n",
    "the natural sciences have “light tails.”\n",
    "\n",
    "To explain this concept, let’s look first at examples."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10de71f5",
   "metadata": {},
   "source": [
    "### \n",
    "\n",
    "The classic example is the [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution), which has density\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{1}{\\sqrt{2\\pi}\\sigma} \n",
    "\\exp\\left( -\\frac{(x-\\mu)^2}{2 \\sigma^2} \\right)\n",
    "\\qquad\n",
    "(-\\infty < x < \\infty)\n",
    "$$\n",
    "\n",
    "The two parameters $ \\mu $ and $ \\sigma $ are the mean and standard deviation\n",
    "respectively.\n",
    "\n",
    "As $ x $ deviates from $ \\mu $, the value of $ f(x) $ goes to zero extremely\n",
    "quickly.\n",
    "\n",
    "We can see this when we plot the density and show a histogram of observations,\n",
    "as with the following code (which assumes $ \\mu=0 $ and $ \\sigma=1 $)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7112b58",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "X = norm.rvs(size=1_000_000)\n",
    "ax.hist(X, bins=40, alpha=0.4, label='histogram', density=True)\n",
    "x_grid = np.linspace(-4, 4, 400)\n",
    "ax.plot(x_grid, norm.pdf(x_grid), label='density')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e317a5d6",
   "metadata": {},
   "source": [
    "Notice how\n",
    "\n",
    "- the density’s tails converge quickly to zero in both directions and  \n",
    "- even with 1,000,000 draws, we get no very large or very small observations.  \n",
    "\n",
    "\n",
    "We can see the last point more clearly by executing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c71c0f7f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "X.min(), X.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b60f36e",
   "metadata": {},
   "source": [
    "Here’s another view of draws from the same distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f499373",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n = 2000\n",
    "fig, ax = plt.subplots()\n",
    "data = norm.rvs(size=n)\n",
    "ax.plot(list(range(n)), data, linestyle='', marker='o', alpha=0.5, ms=4)\n",
    "ax.vlines(list(range(n)), 0, data, lw=0.2)\n",
    "ax.set_ylim(-15, 15)\n",
    "ax.set_xlabel('$i$')\n",
    "ax.set_ylabel('$X_i$', rotation=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "226848e4",
   "metadata": {},
   "source": [
    "We have plotted each individual draw $ X_i $ against $ i $.\n",
    "\n",
    "None are very large or very small.\n",
    "\n",
    "In other words, extreme observations are rare and draws tend not to deviate\n",
    "too much from the mean.\n",
    "\n",
    "Putting this another way, light-tailed distributions are those that\n",
    "rarely generate extreme values.\n",
    "\n",
    "(A more formal definition is given [below](#heavy-tail-formal-definition).)\n",
    "\n",
    "Many statisticians and econometricians\n",
    "use rules of thumb such as “outcomes more than four or five\n",
    "standard deviations from the mean can safely be ignored.”\n",
    "\n",
    "But this is only true when distributions have light tails."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7491ab42",
   "metadata": {},
   "source": [
    "### When are light tails valid?\n",
    "\n",
    "In probability theory and in the real world, many distributions are\n",
    "light-tailed.\n",
    "\n",
    "For example, human height is light-tailed.\n",
    "\n",
    "Yes, it’s true that we see some very tall people.\n",
    "\n",
    "- For example, basketballer [Sun Mingming](https://en.wikipedia.org/wiki/Sun_Mingming) is 2.32 meters tall  \n",
    "\n",
    "\n",
    "But have you ever heard of someone who is 20 meters tall?  Or 200?  Or 2000?\n",
    "\n",
    "Have you ever wondered why not?\n",
    "\n",
    "After all, there are 8 billion people in the world!\n",
    "\n",
    "In essence, the reason we don’t see such draws is that the distribution of\n",
    "human height has very light tails.\n",
    "\n",
    "In fact the distribution of human height obeys a bell-shaped curve similar to the normal distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b92d5ad1",
   "metadata": {},
   "source": [
    "### Returns on assets\n",
    "\n",
    "But what about economic data?\n",
    "\n",
    "Let’s look at some financial data first.\n",
    "\n",
    "Our aim is to plot the daily change in the price of Amazon (AMZN) stock for\n",
    "the period from 1st January 2015 to 1st July 2022.\n",
    "\n",
    "This equates to daily returns if we set dividends aside.\n",
    "\n",
    "The code below produces the desired plot using Yahoo financial data via the `yfinance` library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e690b00c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "data = yf.download('AMZN', '2015-1-1', '2022-7-1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0afc4678",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "s = data['Adj Close']\n",
    "r = s.pct_change()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(r, linestyle='', marker='o', alpha=0.5, ms=4)\n",
    "ax.vlines(r.index, 0, r.values, lw=0.2)\n",
    "ax.set_ylabel('returns', fontsize=12)\n",
    "ax.set_xlabel('date', fontsize=12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fe4a3f6",
   "metadata": {},
   "source": [
    "This data looks different to the draws from the normal distribution we saw above.\n",
    "\n",
    "Several of observations are quite extreme.\n",
    "\n",
    "We get a similar picture if we look at other assets, such as Bitcoin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91fd0989",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "data = yf.download('BTC-USD', '2015-1-1', '2022-7-1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb7ae356",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "s = data['Adj Close']\n",
    "r = s.pct_change()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(r, linestyle='', marker='o', alpha=0.5, ms=4)\n",
    "ax.vlines(r.index, 0, r.values, lw=0.2)\n",
    "ax.set_ylabel('returns', fontsize=12)\n",
    "ax.set_xlabel('date', fontsize=12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85b2f389",
   "metadata": {},
   "source": [
    "The histogram also looks different to the histogram of the normal\n",
    "distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "975aaf6c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "r = np.random.standard_t(df=5, size=1000)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(r, bins=60, alpha=0.4, label='bitcoin returns', density=True)\n",
    "\n",
    "xmin, xmax = plt.xlim()\n",
    "x = np.linspace(xmin, xmax, 100)\n",
    "p = norm.pdf(x, np.mean(r), np.std(r))\n",
    "ax.plot(x, p, linewidth=2, label='normal distribution')\n",
    "\n",
    "ax.set_xlabel('returns', fontsize=12)\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1ad084a",
   "metadata": {},
   "source": [
    "If we look at higher frequency returns data (e.g., tick-by-tick), we often see\n",
    "even more extreme observations.\n",
    "\n",
    "See, for example, [[Mandelbrot, 1963](https://intro.quantecon.org/zreferences.html#id86)] or [[Rachev, 2003](https://intro.quantecon.org/zreferences.html#id85)]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f60fdfff",
   "metadata": {},
   "source": [
    "### Other data\n",
    "\n",
    "The data we have just seen is said to be “heavy-tailed”.\n",
    "\n",
    "With heavy-tailed distributions, extreme outcomes occur relatively\n",
    "frequently."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5114a36",
   "metadata": {},
   "source": [
    "### \n",
    "\n",
    "Importantly, there are many examples of heavy-tailed distributions\n",
    "observed in economic and financial settings!\n",
    "\n",
    "For example, the income and the wealth distributions are heavy-tailed\n",
    "\n",
    "- You can imagine this: most people have low or modest wealth but some people\n",
    "  are extremely rich.  \n",
    "\n",
    "\n",
    "The firm size distribution is also heavy-tailed\n",
    "\n",
    "- You can imagine this too: most firms are small but some firms are enormous.  \n",
    "\n",
    "\n",
    "The distribution of town and city sizes is heavy-tailed\n",
    "\n",
    "- Most towns and cities are small but some are very large.  \n",
    "\n",
    "\n",
    "Later in this lecture, we examine heavy tails in these distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4c3db82",
   "metadata": {},
   "source": [
    "### Why should we care?\n",
    "\n",
    "Heavy tails are common in economic data but does that mean they are important?\n",
    "\n",
    "The answer to this question is affirmative!\n",
    "\n",
    "When distributions are heavy-tailed, we need to think carefully about issues\n",
    "like\n",
    "\n",
    "- diversification and risk  \n",
    "- forecasting  \n",
    "- taxation (across a heavy-tailed income distribution), etc.  \n",
    "\n",
    "\n",
    "We return to these points [below](#heavy-tail-application)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0f079f1",
   "metadata": {},
   "source": [
    "## Visual comparisons\n",
    "\n",
    "In this section, we will introduce important concepts such as the Pareto distribution, Counter CDFs, and Power laws, which aid in recognizing heavy-tailed distributions.\n",
    "\n",
    "Later we will provide a mathematical definition of the difference between\n",
    "light and heavy tails.\n",
    "\n",
    "But for now let’s do some visual comparisons to help us build intuition on the\n",
    "difference between these two types of distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8cdda7e",
   "metadata": {},
   "source": [
    "### Simulations\n",
    "\n",
    "The figure below shows a simulation.\n",
    "\n",
    "The top two subfigures each show 120 independent draws from the normal\n",
    "distribution, which is light-tailed.\n",
    "\n",
    "The bottom subfigure shows 120 independent draws from [the Cauchy\n",
    "distribution](https://en.wikipedia.org/wiki/Cauchy_distribution), which is\n",
    "heavy-tailed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5377661a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n = 120\n",
    "np.random.seed(11)\n",
    "\n",
    "fig, axes = plt.subplots(3, 1, figsize=(6, 12))\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_ylim((-120, 120))\n",
    "\n",
    "s_vals = 2, 12\n",
    "\n",
    "for ax, s in zip(axes[:2], s_vals):\n",
    "    data = np.random.randn(n) * s\n",
    "    ax.plot(list(range(n)), data, linestyle='', marker='o', alpha=0.5, ms=4)\n",
    "    ax.vlines(list(range(n)), 0, data, lw=0.2)\n",
    "    ax.set_title(f\"draws from $N(0, \\sigma^2)$ with $\\sigma = {s}$\", fontsize=11)\n",
    "\n",
    "ax = axes[2]\n",
    "distribution = cauchy()\n",
    "data = distribution.rvs(n)\n",
    "ax.plot(list(range(n)), data, linestyle='', marker='o', alpha=0.5, ms=4)\n",
    "ax.vlines(list(range(n)), 0, data, lw=0.2)\n",
    "ax.set_title(f\"draws from the Cauchy distribution\", fontsize=11)\n",
    "\n",
    "plt.subplots_adjust(hspace=0.25)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "849c827f",
   "metadata": {},
   "source": [
    "In the top subfigure, the standard deviation of the normal distribution is 2,\n",
    "and the draws are clustered around the mean.\n",
    "\n",
    "In the middle subfigure, the standard deviation is increased to 12 and, as\n",
    "expected, the amount of dispersion rises.\n",
    "\n",
    "The bottom subfigure, with the Cauchy draws, shows a different pattern: tight\n",
    "clustering around the mean for the great majority of observations, combined\n",
    "with a few sudden large deviations from the mean.\n",
    "\n",
    "This is typical of a heavy-tailed distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54bcf26d",
   "metadata": {},
   "source": [
    "### Nonnegative distributions\n",
    "\n",
    "Let’s compare some distributions that only take nonnegative values.\n",
    "\n",
    "One is the exponential distribution, which we discussed in [our lecture on probability and distributions](https://intro.quantecon.org/prob_dist.html).\n",
    "\n",
    "The exponential distribution is a light-tailed distribution.\n",
    "\n",
    "Here are some draws from the exponential distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c02b230",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n = 120\n",
    "np.random.seed(11)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_ylim((0, 50))\n",
    "\n",
    "data = np.random.exponential(size=n)\n",
    "ax.plot(list(range(n)), data, linestyle='', marker='o', alpha=0.5, ms=4)\n",
    "ax.vlines(list(range(n)), 0, data, lw=0.2)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cac9feaa",
   "metadata": {},
   "source": [
    "Another nonnegative distribution is the [Pareto distribution](https://en.wikipedia.org/wiki/Pareto_distribution).\n",
    "\n",
    "If $ X $ has the Pareto distribution, then there are positive constants $ \\bar x $\n",
    "and $ \\alpha $ such that\n",
    "\n",
    "\n",
    "<a id='equation-pareto'></a>\n",
    "$$\n",
    "\\mathbb P\\{X > x\\} =\n",
    "\\begin{cases}\n",
    "    \\left( \\bar x/x \\right)^{\\alpha}\n",
    "        & \\text{ if } x \\geq \\bar x\n",
    "    \\\\\n",
    "    1\n",
    "        & \\text{ if } x < \\bar x\n",
    "\\end{cases} \\tag{21.1}\n",
    "$$\n",
    "\n",
    "The parameter $ \\alpha $ is called the **tail index** and $ \\bar x $ is called the\n",
    "**minimum**.\n",
    "\n",
    "The Pareto distribution is a heavy-tailed distribution.\n",
    "\n",
    "One way that the Pareto distribution arises is as the exponential of an\n",
    "exponential random variable.\n",
    "\n",
    "In particular, if $ X $ is exponentially distributed with rate parameter $ \\alpha $, then\n",
    "\n",
    "$$\n",
    "Y = \\bar x \\exp(X)\n",
    "$$\n",
    "\n",
    "is Pareto-distributed with minimum $ \\bar x $ and tail index $ \\alpha $.\n",
    "\n",
    "Here are some draws from the Pareto distribution with tail index $ 1 $ and minimum\n",
    "$ 1 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dff65ed",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n = 120\n",
    "np.random.seed(11)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_ylim((0, 80))\n",
    "exponential_data = np.random.exponential(size=n)\n",
    "pareto_data = np.exp(exponential_data)\n",
    "ax.plot(list(range(n)), pareto_data, linestyle='', marker='o', alpha=0.5, ms=4)\n",
    "ax.vlines(list(range(n)), 0, pareto_data, lw=0.2)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c6d1060",
   "metadata": {},
   "source": [
    "Notice how extreme outcomes are more common."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aa3e0d0",
   "metadata": {},
   "source": [
    "### Counter CDFs\n",
    "\n",
    "For nonnegative random variables, one way to visualize the difference between\n",
    "light and heavy tails is to look at the\n",
    "**counter CDF** (CCDF).\n",
    "\n",
    "For a random variable $ X $ with CDF $ F $, the CCDF is the function\n",
    "\n",
    "$$\n",
    "G(x) := 1 - F(x) = \\mathbb P\\{X > x\\}\n",
    "$$\n",
    "\n",
    "(Some authors call $ G $ the “survival” function.)\n",
    "\n",
    "The CCDF shows how fast the upper tail goes to zero as $ x \\to \\infty $.\n",
    "\n",
    "If $ X $ is exponentially distributed with rate parameter $ \\alpha $, then the CCDF is\n",
    "\n",
    "$$\n",
    "G_E(x) = \\exp(- \\alpha x)\n",
    "$$\n",
    "\n",
    "This function goes to zero relatively quickly as $ x $ gets large.\n",
    "\n",
    "The standard Pareto distribution, where $ \\bar x = 1 $, has CCDF\n",
    "\n",
    "$$\n",
    "G_P(x) = x^{- \\alpha}\n",
    "$$\n",
    "\n",
    "This function goes to zero as $ x \\to \\infty $, but much slower than $ G_E $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e2f9a0a",
   "metadata": {},
   "source": [
    "### Exercise 21.1\n",
    "\n",
    "Show how the CCDF of the standard Pareto distribution can be derived from the CCDF of the exponential distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cad39043",
   "metadata": {},
   "source": [
    "### Solution to[ Exercise 21.1](https://intro.quantecon.org/#ht_ex_x1)\n",
    "\n",
    "Letting $ G_E $ and $ G_P $ be defined as above, letting $ X $ be exponentially\n",
    "distributed with rate parameter $ \\alpha $, and letting $ Y = \\exp(X) $, we have\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    " G_P(y) & = \\mathbb P\\{Y > y\\} \\\\\n",
    "         & = \\mathbb P\\{\\exp(X) > y\\} \\\\\n",
    "         & = \\mathbb P\\{X > \\ln y\\} \\\\\n",
    "         & = G_E(\\ln y) \\\\\n",
    "         & = \\exp( - \\alpha \\ln y) \\\\\n",
    "        & = y^{-\\alpha}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Here’s a plot that illustrates how $ G_E $ goes to zero faster than $ G_P $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fb82858",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "x = np.linspace(1.5, 100, 1000)\n",
    "fig, ax = plt.subplots()\n",
    "alpha = 1.0\n",
    "ax.plot(x, np.exp(- alpha * x), label='exponential', alpha=0.8)\n",
    "ax.plot(x, x**(- alpha), label='Pareto', alpha=0.8)\n",
    "ax.set_xlabel('X value')\n",
    "ax.set_ylabel('CCDF')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e333b38",
   "metadata": {},
   "source": [
    "Here’s a log-log plot of the same functions, which makes visual comparison\n",
    "easier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2fa8df2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "alpha = 1.0\n",
    "ax.loglog(x, np.exp(- alpha * x), label='exponential', alpha=0.8)\n",
    "ax.loglog(x, x**(- alpha), label='Pareto', alpha=0.8)\n",
    "ax.set_xlabel('log value')\n",
    "ax.set_ylabel('log prob')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bca35bf",
   "metadata": {},
   "source": [
    "In the log-log plot, the Pareto CCDF is linear, while the exponential one is\n",
    "concave.\n",
    "\n",
    "This idea is often used to separate light- and heavy-tailed distributions in\n",
    "visualisations — we return to this point below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4442980e",
   "metadata": {},
   "source": [
    "### Empirical CCDFs\n",
    "\n",
    "The sample counterpart of the CCDF function is the **empirical CCDF**.\n",
    "\n",
    "Given a sample $ x_1, \\ldots, x_n $, the empirical CCDF is given by\n",
    "\n",
    "$$\n",
    "\\hat G(x) = \\frac{1}{n} \\sum_{i=1}^n \\mathbb 1\\{x_i > x\\}\n",
    "$$\n",
    "\n",
    "Thus, $ \\hat G(x) $ shows the fraction of the sample that exceeds $ x $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94c607f8",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def eccdf(x, data):\n",
    "    \"Simple empirical CCDF function.\"\n",
    "    return np.mean(data > x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48b86856",
   "metadata": {},
   "source": [
    "Here’s a figure containing some empirical CCDFs from simulated data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1610a145",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Parameters and grid\n",
    "x_grid = np.linspace(1, 1000, 1000)\n",
    "sample_size = 1000\n",
    "np.random.seed(13)\n",
    "z = np.random.randn(sample_size)\n",
    "\n",
    "# Draws\n",
    "data_exp = np.random.exponential(size=sample_size)\n",
    "data_logn = np.exp(z)\n",
    "data_pareto = np.exp(np.random.exponential(size=sample_size))\n",
    "\n",
    "data_list = [data_exp, data_logn, data_pareto]\n",
    "\n",
    "# Build figure\n",
    "fig, axes = plt.subplots(3, 1, figsize=(6, 8))\n",
    "axes = axes.flatten()\n",
    "labels = ['exponential', 'lognormal', 'Pareto']\n",
    "\n",
    "for data, label, ax in zip(data_list, labels, axes):\n",
    "\n",
    "    ax.loglog(x_grid, [eccdf(x, data) for x in x_grid], \n",
    "        'o', markersize=3.0, alpha=0.5, label=label)\n",
    "    ax.set_xlabel(\"log value\")\n",
    "    ax.set_ylabel(\"log prob\")\n",
    "    \n",
    "    ax.legend()\n",
    "    \n",
    "    \n",
    "fig.subplots_adjust(hspace=0.4)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4facf7f9",
   "metadata": {},
   "source": [
    "As with the CCDF, the empirical CCDF from the Pareto distributions is\n",
    "approximately linear in a log-log plot.\n",
    "\n",
    "We will use this idea [below](https://intro.quantecon.org/heavy_tails.html#heavy-tails-in-economic-cross-sections) when we look at real data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8767cbf",
   "metadata": {},
   "source": [
    "#### Q-Q Plots\n",
    "\n",
    "We can also use a [qq plot](https://en.wikipedia.org/wiki/Q%E2%80%93Q_plot) to do a visual comparison between two probability distributions.\n",
    "\n",
    "The [statsmodels](https://www.statsmodels.org/stable/index.html) package provides a convenient [qqplot](https://www.statsmodels.org/stable/generated/statsmodels.graphics.gofplots.qqplot.html) function that, by default, compares sample data to the quintiles of the normal distribution.\n",
    "\n",
    "If the data is drawn from a normal distribution, the plot would look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3e91f95",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "data_normal = np.random.normal(size=sample_size)\n",
    "sm.qqplot(data_normal, line='45')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abbd8135",
   "metadata": {},
   "source": [
    "We can now compare this with the exponential, log-normal, and Pareto distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ce68310",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Build figure\n",
    "fig, axes = plt.subplots(1, 3, figsize=(12, 4))\n",
    "axes = axes.flatten()\n",
    "labels = ['exponential', 'lognormal', 'Pareto']\n",
    "for data, label, ax in zip(data_list, labels, axes):\n",
    "    sm.qqplot(data, line='45', ax=ax, )\n",
    "    ax.set_title(label)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e7dcc12",
   "metadata": {},
   "source": [
    "### Power laws\n",
    "\n",
    "One specific class of heavy-tailed distributions has been found repeatedly in\n",
    "economic and social phenomena: the class of so-called power laws.\n",
    "\n",
    "A random variable $ X $ is said to have a **power law** if, for some $ \\alpha > 0 $,\n",
    "\n",
    "$$\n",
    "\\mathbb P\\{X > x\\} \\approx  x^{-\\alpha}\n",
    "\\quad \\text{when \\$x\\$ is large}\n",
    "$$\n",
    "\n",
    "We can write this more mathematically as\n",
    "\n",
    "\n",
    "<a id='equation-plrt'></a>\n",
    "$$\n",
    "\\lim_{x \\to \\infty} x^\\alpha \\, \\mathbb P\\{X > x\\} = c\n",
    "\\quad \\text{for some \\$c > 0\\$} \\tag{21.2}\n",
    "$$\n",
    "\n",
    "It is also common to say that a random variable $ X $ with this property\n",
    "has a **Pareto tail** with **tail index** $ \\alpha $.\n",
    "\n",
    "Notice that every Pareto distribution with tail index $ \\alpha $\n",
    "has a **Pareto tail** with **tail index** $ \\alpha $.\n",
    "\n",
    "We can think of power laws as a generalization of Pareto distributions.\n",
    "\n",
    "They are distributions that resemble Pareto distributions in their upper right\n",
    "tail.\n",
    "\n",
    "Another way to think of power laws is a set of distributions with a specific\n",
    "kind of (very) heavy tail."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1044b1d9",
   "metadata": {},
   "source": [
    "## Heavy tails in economic cross-sections\n",
    "\n",
    "As mentioned above, heavy tails are pervasive in economic data.\n",
    "\n",
    "In fact power laws seem to be very common as well.\n",
    "\n",
    "We now illustrate this by showing the empirical CCDF of heavy tails.\n",
    "\n",
    "All plots are in log-log, so that a power law shows up as a linear log-log\n",
    "plot, at least in the upper tail.\n",
    "\n",
    "We hide the code that generates the figures, which is somewhat complex, but\n",
    "readers are of course welcome to explore the code (perhaps after examining the figures)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d614e5b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def empirical_ccdf(data, \n",
    "                   ax, \n",
    "                   aw=None,   # weights\n",
    "                   label=None,\n",
    "                   xlabel=None,\n",
    "                   add_reg_line=False, \n",
    "                   title=None):\n",
    "    \"\"\"\n",
    "    Take data vector and return prob values for plotting.\n",
    "    Upgraded empirical_ccdf\n",
    "    \"\"\"\n",
    "    y_vals = np.empty_like(data, dtype='float64')\n",
    "    p_vals = np.empty_like(data, dtype='float64')\n",
    "    n = len(data)\n",
    "    if aw is None:\n",
    "        for i, d in enumerate(data):\n",
    "            # record fraction of sample above d\n",
    "            y_vals[i] = np.sum(data >= d) / n\n",
    "            p_vals[i] = np.sum(data == d) / n\n",
    "    else:\n",
    "        fw = np.empty_like(aw, dtype='float64')\n",
    "        for i, a in enumerate(aw):\n",
    "            fw[i] = a / np.sum(aw)\n",
    "        pdf = lambda x: np.interp(x, data, fw)\n",
    "        data = np.sort(data)\n",
    "        j = 0\n",
    "        for i, d in enumerate(data):\n",
    "            j += pdf(d)\n",
    "            y_vals[i] = 1- j\n",
    "\n",
    "    x, y = np.log(data), np.log(y_vals)\n",
    "    \n",
    "    results = sm.OLS(y, sm.add_constant(x)).fit()\n",
    "    b, a = results.params\n",
    "    \n",
    "    kwargs = [('alpha', 0.3)]\n",
    "    if label:\n",
    "        kwargs.append(('label', label))\n",
    "    kwargs = dict(kwargs)\n",
    "\n",
    "    ax.scatter(x, y, **kwargs)\n",
    "    if add_reg_line:\n",
    "        ax.plot(x, x * a + b, 'k-', alpha=0.6, label=f\"slope = ${a: 1.2f}$\")\n",
    "    if not xlabel:\n",
    "        xlabel='log value'\n",
    "    ax.set_xlabel(xlabel, fontsize=12)\n",
    "    ax.set_ylabel(\"log prob\", fontsize=12)\n",
    "        \n",
    "    if label:\n",
    "        ax.legend(loc='lower left', fontsize=12)\n",
    "        \n",
    "    if title:\n",
    "        ax.set_title(title)\n",
    "        \n",
    "    return np.log(data), y_vals, p_vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc3422b4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def extract_wb(varlist=['NY.GDP.MKTP.CD'], \n",
    "               c='all_countries', \n",
    "               s=1900, \n",
    "               e=2021, \n",
    "               varnames=None):\n",
    "    if c == \"all_countries\":\n",
    "        # Keep countries only (no aggregated regions)\n",
    "        countries = wb.get_countries()\n",
    "        countries_name = countries[countries['region'] != 'Aggregates']['name'].values\n",
    "        c = \"all\"\n",
    "    \n",
    "    df = wb.download(indicator=varlist, country=c, start=s, end=e).stack().unstack(0).reset_index()\n",
    "    df = df.drop(['level_1'], axis=1).transpose()\n",
    "    if varnames is not None:\n",
    "        df.columns = varnames\n",
    "        df = df[1:]\n",
    "    \n",
    "    df1 =df[df.index.isin(countries_name)]\n",
    "    return df1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d418659",
   "metadata": {},
   "source": [
    "### Firm size\n",
    "\n",
    "Here is a plot of the firm size distribution for the largest 500 firms in 2020 taken from Forbes Global 2000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7283578",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "df_fs = pd.read_csv('https://media.githubusercontent.com/media/QuantEcon/high_dim_data/main/cross_section/forbes-global2000.csv')\n",
    "df_fs = df_fs[['Country', 'Sales', 'Profits', 'Assets', 'Market Value']]\n",
    "fig, ax = plt.subplots(figsize=(6.4, 3.5))\n",
    "\n",
    "label=\"firm size (market value)\"\n",
    "top = 500 # set the cutting for top\n",
    "d = df_fs.sort_values('Market Value', ascending=False)\n",
    "empirical_ccdf(np.asarray(d['Market Value'])[:top], ax, label=label, add_reg_line=True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5202735",
   "metadata": {},
   "source": [
    "### City size\n",
    "\n",
    "Here are plots of the city size distribution for the US and Brazil in 2023 from the World Population Review.\n",
    "\n",
    "The size is measured by population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b630f06e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# import population data of cities in 2023 United States and 2023 Brazil from world population review\n",
    "df_cs_us = pd.read_csv('https://media.githubusercontent.com/media/QuantEcon/high_dim_data/main/cross_section/cities_us.csv')\n",
    "df_cs_br = pd.read_csv('https://media.githubusercontent.com/media/QuantEcon/high_dim_data/main/cross_section/cities_brazil.csv')\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(8.8, 3.6))\n",
    "\n",
    "empirical_ccdf(np.asarray(df_cs_us[\"pop2023\"]), axes[0], label=\"US\", add_reg_line=True)\n",
    "empirical_ccdf(np.asarray(df_cs_br['pop2023']), axes[1], label=\"Brazil\", add_reg_line=True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4a6af76",
   "metadata": {},
   "source": [
    "### Wealth\n",
    "\n",
    "Here is a plot of the upper tail (top 500) of the wealth distribution.\n",
    "\n",
    "The data is from the Forbes Billionaires list in 2020."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c4f7cd5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "df_w = pd.read_csv('https://media.githubusercontent.com/media/QuantEcon/high_dim_data/main/cross_section/forbes-billionaires.csv')\n",
    "df_w = df_w[['country', 'realTimeWorth', 'realTimeRank']].dropna()\n",
    "df_w = df_w.astype({'realTimeRank': int})\n",
    "df_w = df_w.sort_values('realTimeRank', ascending=True).copy()\n",
    "countries = ['United States', 'Japan', 'India', 'Italy']  \n",
    "N = len(countries)\n",
    "\n",
    "fig, axs = plt.subplots(2, 2, figsize=(8, 6))\n",
    "axs = axs.flatten()\n",
    "\n",
    "for i, c in enumerate(countries):\n",
    "    df_w_c = df_w[df_w['country'] == c].reset_index()\n",
    "    z = np.asarray(df_w_c['realTimeWorth'])\n",
    "    # print('number of the global richest 2000 from '+ c, len(z))\n",
    "    top = 500           # cut-off number: top 500\n",
    "    if len(z) <= top:    \n",
    "        z = z[:top]\n",
    "\n",
    "    empirical_ccdf(z[:top], axs[i], label=c, xlabel='log wealth', add_reg_line=True)\n",
    "    \n",
    "fig.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c23260d",
   "metadata": {},
   "source": [
    "### GDP\n",
    "\n",
    "Of course, not all cross-sectional distributions are heavy-tailed.\n",
    "\n",
    "Here we show cross-country per capita GDP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d45c058",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# get gdp and gdp per capita for all regions and countries in 2021\n",
    "\n",
    "variable_code = ['NY.GDP.MKTP.CD', 'NY.GDP.PCAP.CD']\n",
    "variable_names = ['GDP', 'GDP per capita']\n",
    "\n",
    "df_gdp1 = extract_wb(varlist=variable_code, \n",
    "                     c=\"all_countries\", \n",
    "                     s=\"2021\", \n",
    "                     e=\"2021\", \n",
    "                     varnames=variable_names)\n",
    "df_gdp1.dropna(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5fbdbb3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(8.8, 3.6))\n",
    "\n",
    "for name, ax in zip(variable_names, axes):\n",
    "    empirical_ccdf(np.asarray(df_gdp1[name]).astype(\"float64\"), ax, add_reg_line=False, label=name)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b48c126e",
   "metadata": {},
   "source": [
    "The plot is concave rather than linear, so the distribution has light tails.\n",
    "\n",
    "One reason is that this is data on an aggregate variable, which involves some\n",
    "averaging in its definition.\n",
    "\n",
    "Averaging tends to eliminate extreme outcomes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75c6f9d7",
   "metadata": {},
   "source": [
    "## Failure of the LLN\n",
    "\n",
    "One impact of heavy tails is that sample averages can be poor estimators of\n",
    "the underlying mean of the distribution.\n",
    "\n",
    "To understand this point better, recall [our earlier discussion](https://intro.quantecon.org/lln_clt.html)\n",
    "of the law of large numbers, which considered IID $ X_1, \\ldots, X_n $ with common distribution $ F $\n",
    "\n",
    "If $ \\mathbb E |X_i| $ is finite, then\n",
    "the sample mean $ \\bar X_n := \\frac{1}{n} \\sum_{i=1}^n X_i $ satisfies\n",
    "\n",
    "\n",
    "<a id='equation-lln-as2'></a>\n",
    "$$\n",
    "\\mathbb P \\left\\{ \\bar X_n \\to \\mu \\text{ as } n \\to \\infty \\right\\} = 1 \\tag{21.3}\n",
    "$$\n",
    "\n",
    "where $ \\mu := \\mathbb E X_i = \\int x F(dx) $ is the common mean of the sample.\n",
    "\n",
    "The condition $ \\mathbb E | X_i | = \\int |x| F(dx) < \\infty $ holds\n",
    "in most cases but can fail if the distribution $ F $ is very heavy-tailed.\n",
    "\n",
    "For example, it fails for the Cauchy distribution.\n",
    "\n",
    "Let’s have a look at the behavior of the sample mean in this case, and see\n",
    "whether or not the LLN is still valid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c88c65b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import cauchy\n",
    "\n",
    "np.random.seed(1234)\n",
    "N = 1_000\n",
    "\n",
    "distribution = cauchy()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "data = distribution.rvs(N)\n",
    "\n",
    "# Compute sample mean at each n\n",
    "sample_mean = np.empty(N)\n",
    "for n in range(1, N):\n",
    "    sample_mean[n] = np.mean(data[:n])\n",
    "\n",
    "# Plot\n",
    "ax.plot(range(N), sample_mean, alpha=0.6, label='$\\\\bar{X}_n$')\n",
    "ax.plot(range(N), np.zeros(N), 'k--', lw=0.5)\n",
    "ax.set_xlabel(r\"$n$\")\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "567efe4e",
   "metadata": {},
   "source": [
    "The sequence shows no sign of converging.\n",
    "\n",
    "We return to this point in the exercises.\n",
    "\n",
    "\n",
    "<a id='heavy-tail-application'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eed61b57",
   "metadata": {},
   "source": [
    "## Why do heavy tails matter?\n",
    "\n",
    "We have now seen that\n",
    "\n",
    "1. heavy tails are frequent in economics and  \n",
    "1. the law of large numbers fails when tails are very heavy.  \n",
    "\n",
    "\n",
    "But what about in the real world?  Do heavy tails matter?\n",
    "\n",
    "Let’s briefly discuss why they do."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90d330a5",
   "metadata": {},
   "source": [
    "### Diversification\n",
    "\n",
    "One of the most important ideas in investing is using diversification to\n",
    "reduce risk.\n",
    "\n",
    "This is a very old idea — consider, for example, the expression “don’t put all your eggs in one basket”.\n",
    "\n",
    "To illustrate, consider an investor with one dollar of wealth and a choice over\n",
    "$ n $ assets with payoffs $ X_1, \\ldots, X_n $.\n",
    "\n",
    "Suppose that returns on distinct  assets are\n",
    "independent and each return has  mean $ \\mu $ and variance $ \\sigma^2 $.\n",
    "\n",
    "If the investor puts all wealth in one asset, say, then the expected payoff of the\n",
    "portfolio is $ \\mu $ and the variance is $ \\sigma^2 $.\n",
    "\n",
    "If instead the investor puts share $ 1/n $ of her wealth in each asset, then the portfolio payoff is\n",
    "\n",
    "$$\n",
    "Y_n = \\sum_{i=1}^n \\frac{X_i}{n} = \\frac{1}{n} \\sum_{i=1}^n X_i.\n",
    "$$\n",
    "\n",
    "Try computing the mean and variance.\n",
    "\n",
    "You will find that\n",
    "\n",
    "- The mean is unchanged at $ \\mu $, while  \n",
    "- the variance of the portfolio has fallen to $ \\sigma^2 / n $.  \n",
    "\n",
    "\n",
    "Diversification reduces risk, as expected.\n",
    "\n",
    "But there is a hidden assumption here: the variance of returns is finite.\n",
    "\n",
    "If the distribution is heavy-tailed and the variance is infinite, then this\n",
    "logic is incorrect.\n",
    "\n",
    "For example, we saw above that if every $ X_i $ is Cauchy, then so is $ Y_n $.\n",
    "\n",
    "This means that diversification doesn’t help at all!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e06d3b7",
   "metadata": {},
   "source": [
    "### Fiscal policy\n",
    "\n",
    "The heaviness of the tail in the wealth distribution matters for taxation and redistribution policies.\n",
    "\n",
    "The same is true for the income distribution.\n",
    "\n",
    "For example, the heaviness of the tail of the income distribution helps\n",
    "determine [how much revenue a given tax policy will raise](https://intro.quantecon.org/mle.html).\n",
    "\n",
    "\n",
    "<a id='cltail'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5816cece",
   "metadata": {},
   "source": [
    "## Classifying tail properties\n",
    "\n",
    "Up until now we have discussed light and heavy tails without any mathematical\n",
    "definitions.\n",
    "\n",
    "Let’s now rectify this.\n",
    "\n",
    "We will focus our attention on the right hand tails of\n",
    "nonnegative random variables and their distributions.\n",
    "\n",
    "The definitions for\n",
    "left hand tails are very similar and we omit them to simplify the exposition.\n",
    "\n",
    "\n",
    "<a id='heavy-tail-formal-definition'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4640e77e",
   "metadata": {},
   "source": [
    "### Light and heavy tails\n",
    "\n",
    "A distribution $ F $ with density $ f $ on $ \\mathbb R_+ $ is called [heavy-tailed](https://en.wikipedia.org/wiki/Heavy-tailed_distribution) if\n",
    "\n",
    "\n",
    "<a id='equation-defht'></a>\n",
    "$$\n",
    "\\int_0^\\infty \\exp(tx) f(x) dx = \\infty \\; \\text{ for all } t > 0. \\tag{21.4}\n",
    "$$\n",
    "\n",
    "We say that a nonnegative random variable $ X $ is **heavy-tailed** if its density is heavy-tailed.\n",
    "\n",
    "This is equivalent to stating that its **moment generating function** $ m(t) :=\n",
    "\\mathbb E \\exp(t X) $ is infinite for all $ t > 0 $.\n",
    "\n",
    "For example, the [log-normal\n",
    "distribution](https://en.wikipedia.org/wiki/Log-normal_distribution) is\n",
    "heavy-tailed because its moment generating function is infinite everywhere on\n",
    "$ (0, \\infty) $.\n",
    "\n",
    "The Pareto distribution is also heavy-tailed.\n",
    "\n",
    "Less formally, a heavy-tailed distribution is one that is not exponentially bounded (i.e. the tails are heavier than the exponential distribution).\n",
    "\n",
    "A distribution $ F $ on $ \\mathbb R_+ $ is called **light-tailed** if it is not heavy-tailed.\n",
    "\n",
    "A nonnegative random variable $ X $ is **light-tailed** if its distribution $ F $ is light-tailed.\n",
    "\n",
    "For example, every random variable with bounded support is light-tailed. (Why?)\n",
    "\n",
    "As another example, if $ X $ has the [exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution), with cdf $ F(x) = 1 - \\exp(-\\lambda x) $ for some $ \\lambda > 0 $, then its moment generating function is\n",
    "\n",
    "$$\n",
    "m(t) = \\frac{\\lambda}{\\lambda - t} \\quad \\text{when } t < \\lambda\n",
    "$$\n",
    "\n",
    "In particular, $ m(t) $ is finite whenever $ t < \\lambda $, so $ X $ is light-tailed.\n",
    "\n",
    "One can show that if $ X $ is light-tailed, then all of its\n",
    "[moments](https://en.wikipedia.org/wiki/Moment_%28mathematics%29) are finite.\n",
    "\n",
    "Conversely, if some moment is infinite, then $ X $ is heavy-tailed.\n",
    "\n",
    "The latter condition is not necessary, however.\n",
    "\n",
    "For example, the lognormal distribution is heavy-tailed but every moment is finite."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5ea838a",
   "metadata": {},
   "source": [
    "## Further reading\n",
    "\n",
    "For more on heavy tails in the wealth distribution, see e.g., [[Vilfredo, 1896](https://intro.quantecon.org/zreferences.html#id90)] and [[Benhabib and Bisin, 2018](https://intro.quantecon.org/zreferences.html#id89)].\n",
    "\n",
    "For more on heavy tails in the firm size distribution, see e.g., [[Axtell, 2001](https://intro.quantecon.org/zreferences.html#id88)], [[Gabaix, 2016](https://intro.quantecon.org/zreferences.html#id87)].\n",
    "\n",
    "For more on heavy tails in the city size distribution, see e.g., [[Rozenfeld *et al.*, 2011](https://intro.quantecon.org/zreferences.html#id84)], [[Gabaix, 2016](https://intro.quantecon.org/zreferences.html#id87)].\n",
    "\n",
    "There are other important implications of heavy tails, aside from those\n",
    "discussed above.\n",
    "\n",
    "For example, heavy tails in income and wealth affect productivity growth, business cycles, and political economy.\n",
    "\n",
    "For further reading, see, for example, [[Acemoglu and Robinson, 2002](https://intro.quantecon.org/zreferences.html#id83)], [[Glaeser *et al.*, 2003](https://intro.quantecon.org/zreferences.html#id82)], [[Bhandari *et al.*, 2018](https://intro.quantecon.org/zreferences.html#id81)] or [[Ahn *et al.*, 2018](https://intro.quantecon.org/zreferences.html#id80)]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b7f59ab",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "423847ce",
   "metadata": {},
   "source": [
    "## Exercise 21.2\n",
    "\n",
    "Prove: If $ X $ has a Pareto tail with tail index $ \\alpha $, then\n",
    "$ \\mathbb E[X^r] = \\infty $ for all $ r \\geq \\alpha $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "219807be",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 21.2](https://intro.quantecon.org/#ht_ex2)\n",
    "\n",
    "Let $ X $ have a Pareto tail with tail index $ \\alpha $ and let $ F $ be its cdf.\n",
    "\n",
    "Fix $ r \\geq \\alpha $.\n",
    "\n",
    "In view of [(21.2)](#equation-plrt), we can take positive constants $ b $ and $ \\bar x $ such that\n",
    "\n",
    "$$\n",
    "\\mathbb P\\{X > x\\} \\geq b x^{- \\alpha} \\text{ whenever } x \\geq \\bar x\n",
    "$$\n",
    "\n",
    "But then\n",
    "\n",
    "$$\n",
    "\\mathbb E X^r = r \\int_0^\\infty x^{r-1} \\mathbb P\\{ X > x \\} dx\n",
    "\\geq\n",
    "r \\int_0^{\\bar x} x^{r-1} \\mathbb P\\{ X > x \\} dx\n",
    "+ r \\int_{\\bar x}^\\infty  x^{r-1} b x^{-\\alpha} dx.\n",
    "$$\n",
    "\n",
    "We know that $ \\int_{\\bar x}^\\infty x^{r-\\alpha-1} dx = \\infty $ whenever $ r - \\alpha - 1 \\geq -1 $.\n",
    "\n",
    "Since $ r \\geq \\alpha $, we have $ \\mathbb E X^r = \\infty $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "930cd1a9",
   "metadata": {},
   "source": [
    "## Exercise 21.3\n",
    "\n",
    "Repeat exercise 1, but replace the three distributions (two normal, one\n",
    "Cauchy) with three Pareto distributions using different choices of\n",
    "$ \\alpha $.\n",
    "\n",
    "For $ \\alpha $, try 1.15, 1.5 and 1.75.\n",
    "\n",
    "Use `np.random.seed(11)` to set the seed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f4b117d",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 21.3](https://intro.quantecon.org/#ht_ex3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1a128cd",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import pareto\n",
    "\n",
    "np.random.seed(11)\n",
    "\n",
    "n = 120\n",
    "alphas = [1.15, 1.50, 1.75]\n",
    "\n",
    "fig, axes = plt.subplots(3, 1, figsize=(6, 8))\n",
    "\n",
    "for (a, ax) in zip(alphas, axes):\n",
    "    ax.set_ylim((-5, 50))\n",
    "    data = pareto.rvs(size=n, scale=1, b=a)\n",
    "    ax.plot(list(range(n)), data, linestyle='', marker='o', alpha=0.5, ms=4)\n",
    "    ax.vlines(list(range(n)), 0, data, lw=0.2)\n",
    "    ax.set_title(f\"Pareto draws with $\\\\alpha = {a}$\", fontsize=11)\n",
    "\n",
    "plt.subplots_adjust(hspace=0.4)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed79e492",
   "metadata": {},
   "source": [
    "## Exercise 21.4\n",
    "\n",
    "There is an ongoing argument about whether the firm size distribution should\n",
    "be modeled as a Pareto distribution or a lognormal distribution (see, e.g.,\n",
    "[[Fujiwara *et al.*, 2004](https://intro.quantecon.org/zreferences.html#id73)], [[Kondo *et al.*, 2018](https://intro.quantecon.org/zreferences.html#id71)] or [[Schluter and Trede, 2019](https://intro.quantecon.org/zreferences.html#id72)]).\n",
    "\n",
    "This sounds esoteric but has real implications for a variety of economic\n",
    "phenomena.\n",
    "\n",
    "To illustrate this fact in a simple way, let us consider an economy with\n",
    "100,000 firms, an interest rate of `r = 0.05` and a corporate tax rate of\n",
    "15%.\n",
    "\n",
    "Your task is to estimate the present discounted value of projected corporate\n",
    "tax revenue over the next 10 years.\n",
    "\n",
    "Because we are forecasting, we need a model.\n",
    "\n",
    "We will suppose that\n",
    "\n",
    "1. the number of firms and the firm size distribution (measured in profits) remain fixed and  \n",
    "1. the firm size distribution is either lognormal or Pareto.  \n",
    "\n",
    "\n",
    "Present discounted value of tax revenue will be estimated by\n",
    "\n",
    "1. generating 100,000 draws of firm profit from the firm size distribution,  \n",
    "1. multiplying by the tax rate, and  \n",
    "1. summing the results with discounting to obtain present value.  \n",
    "\n",
    "\n",
    "The Pareto distribution is assumed to take the form [(21.1)](#equation-pareto) with $ \\bar x = 1 $ and $ \\alpha = 1.05 $.\n",
    "\n",
    "(The value of the tail index $ \\alpha $ is plausible given the data [[Gabaix, 2016](https://intro.quantecon.org/zreferences.html#id87)].)\n",
    "\n",
    "To make the lognormal option as similar as possible to the Pareto option, choose\n",
    "its parameters such that the mean and median of both distributions are the same.\n",
    "\n",
    "Note that, for each distribution, your estimate of tax revenue will be random\n",
    "because it is based on a finite number of draws.\n",
    "\n",
    "To take this into account, generate 100 replications (evaluations of tax revenue)\n",
    "for each of the two distributions and compare the two samples by\n",
    "\n",
    "- producing a [violin plot](https://en.wikipedia.org/wiki/Violin_plot) visualizing the two samples side-by-side and  \n",
    "- printing the mean and standard deviation of both samples.  \n",
    "\n",
    "\n",
    "For the seed use `np.random.seed(1234)`.\n",
    "\n",
    "What differences do you observe?\n",
    "\n",
    "(Note: a better approach to this problem would be to model firm dynamics and\n",
    "try to track individual firms given the current distribution.  We will discuss\n",
    "firm dynamics in later lectures.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "118628ff",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 21.4](https://intro.quantecon.org/#ht_ex5)\n",
    "\n",
    "To do the exercise, we need to choose the parameters $ \\mu $\n",
    "and $ \\sigma $ of the lognormal distribution to match the mean and median\n",
    "of the Pareto distribution.\n",
    "\n",
    "Here we understand the lognormal distribution as that of the random variable\n",
    "$ \\exp(\\mu + \\sigma Z) $ when $ Z $ is standard normal.\n",
    "\n",
    "The mean and median of the Pareto distribution [(21.1)](#equation-pareto) with\n",
    "$ \\bar x = 1 $ are\n",
    "\n",
    "$$\n",
    "\\text{mean } = \\frac{\\alpha}{\\alpha - 1}\n",
    "\\quad \\text{and} \\quad\n",
    "\\text{median } = 2^{1/\\alpha}\n",
    "$$\n",
    "\n",
    "Using the corresponding expressions for the lognormal distribution leads us to\n",
    "the equations\n",
    "\n",
    "$$\n",
    "\\frac{\\alpha}{\\alpha - 1} = \\exp(\\mu + \\sigma^2/2)\n",
    "\\quad \\text{and} \\quad\n",
    "2^{1/\\alpha} = \\exp(\\mu)\n",
    "$$\n",
    "\n",
    "which we solve for $ \\mu $ and $ \\sigma $ given $ \\alpha = 1.05 $.\n",
    "\n",
    "Here is the code that generates the two samples, produces the violin plot and\n",
    "prints the mean and standard deviation of the two samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f59ce130",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "num_firms = 100_000\n",
    "num_years = 10\n",
    "tax_rate = 0.15\n",
    "r = 0.05\n",
    "\n",
    "β = 1 / (1 + r)    # discount factor\n",
    "\n",
    "x_bar = 1.0\n",
    "α = 1.05\n",
    "\n",
    "def pareto_rvs(n):\n",
    "    \"Uses a standard method to generate Pareto draws.\"\n",
    "    u = np.random.uniform(size=n)\n",
    "    y = x_bar / (u**(1/α))\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa3917f4",
   "metadata": {},
   "source": [
    "Let’s compute the lognormal parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aea2f606",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "μ = np.log(2) / α\n",
    "σ_sq = 2 * (np.log(α/(α - 1)) - np.log(2)/α)\n",
    "σ = np.sqrt(σ_sq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41572041",
   "metadata": {},
   "source": [
    "Here’s a function to compute a single estimate of tax revenue for a particular\n",
    "choice of distribution `dist`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b20ffb2c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def tax_rev(dist):\n",
    "    tax_raised = 0\n",
    "    for t in range(num_years):\n",
    "        if dist == 'pareto':\n",
    "            π = pareto_rvs(num_firms)\n",
    "        else:\n",
    "            π = np.exp(μ + σ * np.random.randn(num_firms))\n",
    "        tax_raised += β**t * np.sum(π * tax_rate)\n",
    "    return tax_raised"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7bd53f3",
   "metadata": {},
   "source": [
    "Now let’s generate the violin plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96bd7690",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "num_reps = 100\n",
    "np.random.seed(1234)\n",
    "\n",
    "tax_rev_lognorm = np.empty(num_reps)\n",
    "tax_rev_pareto = np.empty(num_reps)\n",
    "\n",
    "for i in range(num_reps):\n",
    "    tax_rev_pareto[i] = tax_rev('pareto')\n",
    "    tax_rev_lognorm[i] = tax_rev('lognorm')\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "data = tax_rev_pareto, tax_rev_lognorm\n",
    "\n",
    "ax.violinplot(data)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fdd52c8",
   "metadata": {},
   "source": [
    "Finally, let’s print the means and standard deviations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "453be2c5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "tax_rev_pareto.mean(), tax_rev_pareto.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7b7eb34",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "tax_rev_lognorm.mean(), tax_rev_lognorm.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2a9007c",
   "metadata": {},
   "source": [
    "Looking at the output of the code, our main conclusion is that the Pareto\n",
    "assumption leads to a lower mean and greater dispersion."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8482997",
   "metadata": {},
   "source": [
    "## Exercise 21.5\n",
    "\n",
    "The [characteristic function](https://en.wikipedia.org/wiki/Characteristic_function_%28probability_theory%29) of the Cauchy distribution is\n",
    "\n",
    "\n",
    "<a id='equation-lln-cch'></a>\n",
    "$$\n",
    "\\phi(t) = \\mathbb E e^{itX} = \\int e^{i t x} f(x) dx = e^{-|t|} \\tag{21.5}\n",
    "$$\n",
    "\n",
    "Prove that the sample mean $ \\bar X_n $ of $ n $ independent draws $ X_1, \\ldots,\n",
    "X_n $ from the Cauchy distribution has the same characteristic function as\n",
    "$ X_1 $.\n",
    "\n",
    "(This means that the sample mean never converges.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e682f2c",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 21.5](https://intro.quantecon.org/#ht_ex_cauchy)\n",
    "\n",
    "By independence, the characteristic function of the sample mean becomes\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\mathbb E e^{i t \\bar X_n }\n",
    "    & = \\mathbb E \\exp \\left\\{ i \\frac{t}{n} \\sum_{j=1}^n X_j \\right\\}\n",
    "    \\\\\n",
    "    & = \\mathbb E \\prod_{j=1}^n \\exp \\left\\{ i \\frac{t}{n} X_j \\right\\}\n",
    "    \\\\\n",
    "    & = \\prod_{j=1}^n \\mathbb E \\exp \\left\\{ i \\frac{t}{n} X_j \\right\\}\n",
    "    = [\\phi(t/n)]^n\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "In view of [(21.5)](#equation-lln-cch), this is just $ e^{-|t|} $.\n",
    "\n",
    "Thus, in the case of the Cauchy distribution, the sample mean itself has the very same Cauchy distribution, regardless of $ n $!"
   ]
  }
 ],
 "metadata": {
  "date": 1733269016.1359477,
  "filename": "heavy_tails.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Heavy-Tailed Distributions"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}