heavy_tails.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='heavy-tails'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heavy-Tailed Distributions\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Heavy-Tailed Distributions](#Heavy-Tailed-Distributions)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Visual Comparisons](#Visual-Comparisons)  \n",
    "  - [Failure of the LLN](#Failure-of-the-LLN)  \n",
    "  - [Classifying Tail Properties](#Classifying-Tail-Properties)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon\n",
    "!pip install --upgrade yfinance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Most commonly used probability distributions in classical statistics and\n",
    "the natural sciences have either bounded support or light tails.\n",
    "\n",
    "When a distribution is light-tailed, extreme observations are rare and\n",
    "draws tend not to deviate too much from the mean.\n",
    "\n",
    "Having internalized these kinds of distributions, many researchers and\n",
    "practitioners use rules of thumb such as “outcomes more than four or five\n",
    "standard deviations from the mean can safely be ignored.”\n",
    "\n",
    "However, some distributions encountered in economics have far more probability\n",
    "mass in the tails than distributions like the normal distribution.\n",
    "\n",
    "With such **heavy-tailed** distributions, what would be regarded as extreme\n",
    "outcomes for someone accustomed to thin tailed distributions occur relatively\n",
    "frequently.\n",
    "\n",
    "Examples of heavy-tailed distributions observed in economic and financial\n",
    "settings include\n",
    "\n",
    "- the income distributions and the wealth distribution (see, e.g., [[Vil96]](https://python-programming.quantecon.org/zreferences.html#pareto1896cours), [[BB18]](https://python-programming.quantecon.org/zreferences.html#benhabib2018skewed)),  \n",
    "- the firm size distribution ([[Axt01]](https://python-programming.quantecon.org/zreferences.html#axtell2001zipf), [[Gab16]](https://python-programming.quantecon.org/zreferences.html#gabaix2016power)}),  \n",
    "- the distribution of returns on holding assets over short time horizons ([[Man63]](https://python-programming.quantecon.org/zreferences.html#mandelbrot1963variation), [[Rac03]](https://python-programming.quantecon.org/zreferences.html#rachev2003handbook)), and  \n",
    "- the distribution of city sizes ([[RRGM11]](https://python-programming.quantecon.org/zreferences.html#rozenfeld2011area), [[Gab16]](https://python-programming.quantecon.org/zreferences.html#gabaix2016power)).  \n",
    "\n",
    "\n",
    "These heavy tails turn out to be important for our understanding of economic outcomes.\n",
    "\n",
    "As one example, the heaviness of the tail in the wealth distribution is one\n",
    "natural measure of inequality.\n",
    "\n",
    "It matters for taxation and redistribution\n",
    "policies, as well as for flow-on effects for productivity growth, business\n",
    "cycles, and political economy\n",
    "\n",
    "- see, e.g., [[AR02]](https://python-programming.quantecon.org/zreferences.html#acemoglu2002political), [[GSS03]](https://python-programming.quantecon.org/zreferences.html#glaeser2003injustice), [[BEGS18]](https://python-programming.quantecon.org/zreferences.html#bhandari2018inequality) or [[AKM+18]](https://python-programming.quantecon.org/zreferences.html#ahn2018inequality).  \n",
    "\n",
    "\n",
    "This lecture formalizes some of the concepts introduced above and reviews the\n",
    "key ideas.\n",
    "\n",
    "Let’s start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import quantecon as qe\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following two lines can be added to avoid an annoying FutureWarning, and prevent a specific compatibility issue between pandas and matplotlib from causing problems down the line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from pandas.plotting import register_matplotlib_converters\n",
    "register_matplotlib_converters()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visual Comparisons\n",
    "\n",
    "One way to build intuition on the difference between light and heavy tails is\n",
    "to plot independent draws and compare them side-by-side."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Simulation\n",
    "\n",
    "The figure below shows a simulation.  (You will be asked to replicate it in\n",
    "the exercises.)\n",
    "\n",
    "The top two subfigures each show 120 independent draws from the normal distribution, which is light-tailed.\n",
    "\n",
    "The bottom subfigure shows 120 independent draws from [the Cauchy distribution](https://en.wikipedia.org/wiki/Cauchy_distribution), which is heavy-tailed.\n",
    "\n",
    "\n",
    "<a id='light-heavy-fig1'></a>\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/python-programming.quantecon.org/_static/lecture_specific/heavy_tails/light_heavy_fig1.png\" style=\"\">\n",
    "\n",
    "  \n",
    "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 expected, the amount of dispersion rises.\n",
    "\n",
    "The bottom subfigure, with the Cauchy draws, shows a\n",
    "different pattern: tight clustering around the mean for the great majority of\n",
    "observations, combined with a few sudden large deviations from the mean.\n",
    "\n",
    "This is typical of a heavy-tailed distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Heavy Tails in Asset Returns\n",
    "\n",
    "Next let’s look at some financial data.\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 November 2019.\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,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import yfinance as yf\n",
    "import pandas as pd\n",
    "\n",
    "s = yf.download('AMZN', '2015-1-1', '2019-11-1')['Adj Close']\n",
    "\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",
    "\n",
    "ax.set_ylabel('returns', fontsize=12)\n",
    "ax.set_xlabel('date', fontsize=12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Five of the 1217 observations are more than 5 standard\n",
    "deviations from the mean.\n",
    "\n",
    "Overall, the figure is suggestive of heavy tails,\n",
    "although not to the same degree as the Cauchy distribution the\n",
    "figure above.\n",
    "\n",
    "If, however, one takes tick-by-tick data rather\n",
    "daily data, the heavy-tailedness of the distribution increases further."
   ]
  },
  {
   "cell_type": "markdown",
   "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://python-programming.quantecon.org/lln_clt.html) of the Law of Large Numbers, which considered IID $ X_1,\n",
    "\\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{1}\n",
    "$$\n",
    "\n",
    "where $ \\mu := \\mathbb E X_i = \\int x F(x) $ is the common mean of the sample.\n",
    "\n",
    "The condition $ \\mathbb E | X_i | = \\int |x| F(x) < \\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,
   "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",
    "\n",
    "ax.plot(range(N), np.zeros(N), 'k--', lw=0.5)\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sequence shows no sign of converging.\n",
    "\n",
    "Will convergence occur if we take $ n $ even larger?\n",
    "\n",
    "The answer is no.\n",
    "\n",
    "To see this, recall that 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{2}\n",
    "$$\n",
    "\n",
    "Using 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 [(2)](#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 $!\n",
    "\n",
    "In particular, the sequence $ \\bar X_n $ does not converge to any point.\n",
    "\n",
    "\n",
    "<a id='cltail'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classifying Tail Properties\n",
    "\n",
    "To keep our discussion precise, we need some definitions concerning tail\n",
    "properties.\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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Light and Heavy Tails\n",
    "\n",
    "A distribution $ F $ on $ \\mathbb R_+ $ is called **heavy-tailed** if\n",
    "\n",
    "\n",
    "<a id='equation-defht'></a>\n",
    "$$\n",
    "\\int_0^\\infty \\exp(tx) F(dx) = \\infty \\; \\text{ for all } t > 0. \\tag{3}\n",
    "$$\n",
    "\n",
    "We say that a nonnegative random variable $ X $ is **heavy-tailed** if its distribution $ F(x) := \\mathbb P\\{X \\leq x\\} $ is heavy-tailed.\n",
    "\n",
    "This is equivalent to stating that its **moment generating function**\n",
    "$ m(t) := \\mathbb E \\exp(t X) $ is infinite for all $ t > 0 $.\n",
    "\n",
    "- For example, the lognormal distribution is heavy-tailed because its\n",
    "  moment generating function is infinite everywhere on $ (0, \\infty) $.  \n",
    "\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",
    "- Example: Every random variable with bounded support is light-tailed. (Why?)  \n",
    "- Example: If $ X $ has the exponential distribution, with cdf $ F(x) = 1 - \\exp(-\\lambda x) $ for some $ \\lambda > 0 $, then its moment generating function is finite whenever $ t < \\lambda $.  Hence $ X $ is light-tailed.  \n",
    "\n",
    "\n",
    "One can show that if $ X $ is light-tailed, then all of its moments are finite.\n",
    "\n",
    "The contrapositive is that if some moment is infinite, then $ X $ is heavy-tailed.\n",
    "\n",
    "The latter condition is not necessary, however.\n",
    "\n",
    "- Example: the lognormal distribution is heavy-tailed but every moment is finite.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pareto Tails\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",
    "Specifically, given $ \\alpha > 0 $, a nonnegative random variable $ X $ is said to have a **Pareto tail** with **tail index** $ \\alpha $ if\n",
    "\n",
    "\n",
    "<a id='equation-plrt'></a>\n",
    "$$\n",
    "\\lim_{x \\to \\infty} x^\\alpha \\, \\mathbb P\\{X > x\\} = c. \\tag{4}\n",
    "$$\n",
    "\n",
    "Evidently [(4)](#equation-plrt) implies the existence of positive constants $ b $ and $ \\bar x $ such that $ \\mathbb P\\{X > x\\} \\geq b x^{- \\alpha} $ whenever $ x \\geq \\bar x $.\n",
    "\n",
    "The implication is that $ \\mathbb P\\{X > x\\} $ converges to zero no faster than $ x^{-\\alpha} $.\n",
    "\n",
    "In some sources, a random variable obeying [(4)](#equation-plrt) is said to have a **power law tail**.\n",
    "\n",
    "The primary example is the **Pareto distribution**, which has distribution\n",
    "\n",
    "\n",
    "<a id='equation-pareto'></a>\n",
    "$$\n",
    "F(x) =\n",
    "\\begin{cases}\n",
    "    1 - \\left( \\bar x/x \\right)^{\\alpha}\n",
    "        & \\text{ if } x \\geq \\bar x\n",
    "    \\\\\n",
    "    0\n",
    "        & \\text{ if } x < \\bar x\n",
    "\\end{cases} \\tag{5}\n",
    "$$\n",
    "\n",
    "for some positive constants $ \\bar x $ and $ \\alpha $.\n",
    "\n",
    "It is easy to see that if $ X \\sim F $, then $ \\mathbb P\\{X > x\\} $ satisfies [(4)](#equation-plrt).\n",
    "\n",
    "Thus, in line with the terminology, Pareto distributed random variables have a Pareto tail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rank-Size Plots\n",
    "\n",
    "One graphical technique for investigating Pareto tails and power laws is the so-called **rank-size plot**.\n",
    "\n",
    "This kind of figure plots\n",
    "log size against log rank of the population (i.e., location in the population\n",
    "when sorted from smallest to largest).\n",
    "\n",
    "Often just the largest 5 or 10% of observations are plotted.\n",
    "\n",
    "For a sufficiently large number of draws from a Pareto distribution, the plot generates a straight line. For distributions with thinner tails, the data points are concave.\n",
    "\n",
    "A discussion of why this occurs can be found in [[NOM04]](https://python-programming.quantecon.org/zreferences.html#nishiyama2004estimation).\n",
    "\n",
    "The figure below provides one example, using simulated data.\n",
    "\n",
    "The rank-size plots shows draws from three different distributions: folded normal, chi-squared with 1 degree of freedom and Pareto.\n",
    "\n",
    "The Pareto sample produces a straight line, while the lines produced by the other samples are concave.\n",
    "\n",
    "You are asked to reproduce this figure in the exercises.\n",
    "\n",
    "\n",
    "<a id='rank-size-fig1'></a>\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/python-programming.quantecon.org/_static/lecture_specific/heavy_tails/rank_size_fig1.png\" style=\"\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Replicate [the figure presented above](#light-heavy-fig1) that compares normal and Cauchy draws.\n",
    "\n",
    "Use `np.random.seed(11)` to set the seed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 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",
   "metadata": {},
   "source": [
    "### Exercise 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",
   "metadata": {},
   "source": [
    "### Exercise 4\n",
    "\n",
    "Replicate the rank-size plot figure [presented above](#rank-size-fig1).\n",
    "\n",
    "If you like you can use the function `qe.rank_size_plot` from the `quantecon` library to generate the plots.\n",
    "\n",
    "Use `np.random.seed(13)` to set the seed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5\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",
    "[[FDGA+04]](https://python-programming.quantecon.org/zreferences.html#fujiwara2004pareto), [[KLS18]](https://python-programming.quantecon.org/zreferences.html#kondo2018us) or [[ST19a]](https://python-programming.quantecon.org/zreferences.html#schluter2019size)).\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 [(5)](#equation-pareto) with $ \\bar x = 1 $ and $ \\alpha = 1.05 $.\n",
    "\n",
    "(The value the tail index $ \\alpha $ is plausible given the data [[Gab16]](https://python-programming.quantecon.org/zreferences.html#gabaix2016power).)\n",
    "\n",
    "To make the lognormal option as similar as possible to the Pareto option, choose 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 because it is based on a finite number of draws.\n",
    "\n",
    "To take this into account, generate 100 replications (evaluations of tax revenue) 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",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
   "metadata": {},
   "source": [
    "### Exercise 2\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",
    "As discussed after [(4)](#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 \\} x\n",
    "\\geq\n",
    "r \\int_0^{\\bar x} x^{r-1} \\mathbb P\\{ X > x \\} x\n",
    "+ r \\int_{\\bar x}^\\infty  x^{r-1} b x^{-\\alpha} x.\n",
    "$$\n",
    "\n",
    "We know that $ \\int_{\\bar x}^\\infty x^{r-\\alpha-1} x = \\infty $ whenever $ r - \\alpha - 1 \\geq -1 $.\n",
    "\n",
    "Since $ r \\geq \\alpha $, we have $ \\mathbb E X^r = \\infty $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
   "metadata": {},
   "source": [
    "### Exercise 4\n",
    "\n",
    "First let’s generate the data for the plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "sample_size = 1000\n",
    "np.random.seed(13)\n",
    "z = np.random.randn(sample_size)\n",
    "\n",
    "data_1 = np.abs(z)\n",
    "data_2 = np.exp(z)\n",
    "data_3 = np.exp(np.random.exponential(scale=1.0, size=sample_size))\n",
    "\n",
    "data_list = [data_1, data_2, data_3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, 1, figsize=(6, 8))\n",
    "axes = axes.flatten()\n",
    "labels = ['$|z|$', '$\\exp(z)$', 'Pareto with tail index $1.0$']\n",
    "\n",
    "for data, label, ax in zip(data_list, labels, axes):\n",
    "\n",
    "    qe.rank_size_plot(data, ax, label=label)\n",
    "    ax.legend()\n",
    "\n",
    "fig.subplots_adjust(hspace=0.4)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5\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 [(5)](#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 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,
   "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",
   "metadata": {},
   "source": [
    "Let’s compute the lognormal parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
   "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,
   "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",
   "metadata": {},
   "source": [
    "Now let’s generate the violin plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
   "metadata": {},
   "source": [
    "Finally, let’s print the means and standard deviations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "tax_rev_pareto.mean(), tax_rev_pareto.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "tax_rev_lognorm.mean(), tax_rev_lognorm.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "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."
   ]
  }
 ],
 "metadata": {
  "date": 1591697894.3928902,
  "filename": "heavy_tails.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Heavy-Tailed Distributions"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}