lln_clt.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='lln-clt'></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": [
    "# LLN and CLT\n",
    "\n",
    "\n",
    "<a id='index-3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [LLN and CLT](#LLN-and-CLT)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Relationships](#Relationships)  \n",
    "  - [LLN](#LLN)  \n",
    "  - [CLT](#CLT)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture illustrates two of the most important theorems of probability and statistics: The\n",
    "law of large numbers (LLN) and the central limit theorem (CLT).\n",
    "\n",
    "These beautiful theorems lie behind many of the most fundamental results in econometrics and quantitative economic modeling.\n",
    "\n",
    "The lecture is based around simulations that show the LLN and CLT in action.\n",
    "\n",
    "We also demonstrate how the LLN and CLT break down when the assumptions they are based on do not hold.\n",
    "\n",
    "In addition, we examine several useful extensions of the classical theorems, such as\n",
    "\n",
    "- The delta method, for smooth functions of random variables.  \n",
    "- The multivariate case.  \n",
    "\n",
    "\n",
    "Some of these extensions are presented as exercises.\n",
    "\n",
    "We’ll need the following imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy.stats import t, beta, lognorm, expon, gamma, uniform, cauchy\n",
    "from scipy.stats import gaussian_kde, poisson, binom, norm, chi2\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.collections import PolyCollection\n",
    "from scipy.linalg import inv, sqrtm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relationships\n",
    "\n",
    "The CLT refines the LLN.\n",
    "\n",
    "The LLN gives conditions under which sample moments converge to population moments as sample size increases.\n",
    "\n",
    "The CLT provides information about the rate at which sample moments converge to population moments as sample size increases.\n",
    "\n",
    "\n",
    "<a id='lln-mr'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LLN\n",
    "\n",
    "\n",
    "<a id='index-4'></a>\n",
    "We begin with the law of large numbers, which tells us when sample averages\n",
    "will converge to their population means.\n",
    "\n",
    "\n",
    "<a id='lln-ksl'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Classical LLN\n",
    "\n",
    "The classical law of large numbers concerns independent and\n",
    "identically distributed (IID) random variables.\n",
    "\n",
    "Here is the strongest version of the classical LLN, known as *Kolmogorov’s strong law*.\n",
    "\n",
    "Let $ X_1, \\ldots, X_n $ be independent and identically\n",
    "distributed scalar random variables, with common distribution $ F $.\n",
    "\n",
    "When it exists, let $ \\mu $ denote the common mean of this sample:\n",
    "\n",
    "$$\n",
    "\\mu := \\mathbb E X = \\int x F(dx)\n",
    "$$\n",
    "\n",
    "In addition, let\n",
    "\n",
    "$$\n",
    "\\bar X_n := \\frac{1}{n} \\sum_{i=1}^n X_i\n",
    "$$\n",
    "\n",
    "Kolmogorov’s strong law states that, if $ \\mathbb E |X| $ is finite, then\n",
    "\n",
    "\n",
    "<a id='equation-lln-as'></a>\n",
    "$$\n",
    "\\mathbb P \\left\\{ \\bar X_n \\to \\mu \\text{ as } n \\to \\infty \\right\\} = 1 \\tag{1}\n",
    "$$\n",
    "\n",
    "What does this last expression mean?\n",
    "\n",
    "Let’s think about it from a simulation perspective, imagining for a moment that\n",
    "our computer can generate perfect random samples (which of course [it can’t](https://en.wikipedia.org/wiki/Pseudorandom_number_generator)).\n",
    "\n",
    "Let’s also imagine that we can generate infinite sequences so that the\n",
    "statement $ \\bar X_n \\to \\mu $ can be evaluated.\n",
    "\n",
    "In this setting, [(1)](#equation-lln-as) should be interpreted as meaning that the\n",
    "probability of the computer producing a sequence where $ \\bar X_n \\to \\mu $ fails to occur\n",
    "is zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proof\n",
    "\n",
    "\n",
    "<a id='index-5'></a>\n",
    "The proof of Kolmogorov’s strong law is nontrivial – see, for example, theorem 8.3.5 of [[Dud02]](https://python.quantecon.org/zreferences.html#dudley2002).\n",
    "\n",
    "On the other hand, we can prove a weaker version of the LLN very easily and\n",
    "still get most of the intuition.\n",
    "\n",
    "The version we prove is as follows: If $ X_1, \\ldots, X_n $ is IID with $ \\mathbb E X_i^2 < \\infty $,\n",
    "then, for any $ \\epsilon > 0 $, we have\n",
    "\n",
    "\n",
    "<a id='equation-lln-ip'></a>\n",
    "$$\n",
    "\\mathbb P \\left\\{ | \\bar X_n - \\mu | \\geq \\epsilon \\right\\} \\to 0\n",
    "\\quad \\text{as} \\quad\n",
    "n \\to \\infty \\tag{2}\n",
    "$$\n",
    "\n",
    "(This version is weaker because we claim only [convergence in probability](https://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_probability) rather than [almost sure convergence](https://en.wikipedia.org/wiki/Convergence_of_random_variables#Almost_sure_convergence), and assume a finite second moment)\n",
    "\n",
    "To see that this is so, fix $ \\epsilon > 0 $, and let $ \\sigma^2 $ be the variance of each $ X_i $.\n",
    "\n",
    "Recall the [Chebyshev inequality](https://en.wikipedia.org/wiki/Chebyshev%27s_inequality), which tells us that\n",
    "\n",
    "\n",
    "<a id='equation-lln-cheb'></a>\n",
    "$$\n",
    "\\mathbb P \\left\\{ | \\bar X_n - \\mu | \\geq \\epsilon \\right\\}\n",
    "\\leq \\frac{\\mathbb E [ (\\bar X_n - \\mu)^2]}{\\epsilon^2} \\tag{3}\n",
    "$$\n",
    "\n",
    "Now observe that\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\mathbb E [ (\\bar X_n - \\mu)^2 ]\n",
    "    & = \\mathbb E \\left\\{ \\left[\n",
    "    \\frac{1}{n} \\sum_{i=1}^n (X_i - \\mu)\n",
    "    \\right]^2 \\right\\}\n",
    "    \\\\\n",
    "    & = \\frac{1}{n^2} \\sum_{i=1}^n \\sum_{j=1}^n \\mathbb E (X_i - \\mu)(X_j - \\mu) \\nonumber\n",
    "    \\\\\n",
    "    & = \\frac{1}{n^2} \\sum_{i=1}^n \\mathbb E (X_i - \\mu)^2 \\nonumber\n",
    "    \\\\\n",
    "    & = \\frac{\\sigma^2}{n} \\nonumber\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Here the crucial step is at the third equality, which follows from\n",
    "independence.\n",
    "\n",
    "Independence means that if $ i \\not= j $, then the covariance term $ \\mathbb E (X_i - \\mu)(X_j - \\mu) $ drops out.\n",
    "\n",
    "As a result, $ n^2 - n $ terms vanish, leading us to a final expression that goes to zero in $ n $.\n",
    "\n",
    "Combining our last result with [(3)](#equation-lln-cheb), we come to the estimate\n",
    "\n",
    "\n",
    "<a id='equation-lln-cheb2'></a>\n",
    "$$\n",
    "\\mathbb P \\left\\{ | \\bar X_n - \\mu | \\geq \\epsilon \\right\\}\n",
    "\\leq \\frac{\\sigma^2}{n \\epsilon^2} \\tag{4}\n",
    "$$\n",
    "\n",
    "The claim in [(2)](#equation-lln-ip) is now clear.\n",
    "\n",
    "Of course, if the sequence $ X_1, \\ldots, X_n $ is correlated, then the cross-product terms\n",
    "$ \\mathbb E (X_i - \\mu)(X_j - \\mu) $ are not necessarily zero.\n",
    "\n",
    "While this doesn’t mean that the same line of argument is impossible, it does mean\n",
    "that if we want a similar result then the covariances should be “almost zero”\n",
    "for “most” of these terms.\n",
    "\n",
    "In a long sequence, this would be true if, for example, $ \\mathbb E (X_i - \\mu)(X_j - \\mu) $\n",
    "approached zero when the difference between $ i $ and $ j $ became\n",
    "large.\n",
    "\n",
    "In other words, the LLN can still work if the sequence $ X_1, \\ldots, X_n $ has a kind of “asymptotic independence”, in the sense that correlation falls to zero as variables become further apart in the sequence.\n",
    "\n",
    "This idea is very important in time series analysis, and we’ll come across it again soon enough."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Illustration\n",
    "\n",
    "\n",
    "<a id='index-6'></a>\n",
    "Let’s now illustrate the classical IID law of large numbers using simulation.\n",
    "\n",
    "In particular, we aim to generate some sequences of IID random variables and plot the evolution\n",
    "of $ \\bar X_n $ as $ n $ increases.\n",
    "\n",
    "Below is a figure that does just this (as usual, you can click on it to expand it).\n",
    "\n",
    "It shows IID observations from three different distributions and plots $ \\bar X_n $ against $ n $ in each case.\n",
    "\n",
    "The dots represent the underlying observations $ X_i $ for $ i = 1, \\ldots, 100 $.\n",
    "\n",
    "In each of the three cases, convergence of $ \\bar X_n $ to $ \\mu $ occurs as predicted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n = 100\n",
    "\n",
    "# Arbitrary collection of distributions\n",
    "distributions = {\"student's t with 10 degrees of freedom\": t(10),\n",
    "                 \"β(2, 2)\": beta(2, 2),\n",
    "                 \"lognormal LN(0, 1/2)\": lognorm(0.5),\n",
    "                 \"γ(5, 1/2)\": gamma(5, scale=2),\n",
    "                 \"poisson(4)\": poisson(4),\n",
    "                 \"exponential with λ = 1\": expon(1)}\n",
    "\n",
    "# Create a figure and some axes\n",
    "num_plots = 3\n",
    "fig, axes = plt.subplots(num_plots, 1, figsize=(10, 20))\n",
    "\n",
    "# Set some plotting parameters to improve layout\n",
    "bbox = (0., 1.02, 1., .102)\n",
    "legend_args = {'ncol': 2,\n",
    "               'bbox_to_anchor': bbox,\n",
    "               'loc': 3,\n",
    "               'mode': 'expand'}\n",
    "plt.subplots_adjust(hspace=0.5)\n",
    "\n",
    "for ax in axes:\n",
    "    # Choose a randomly selected distribution\n",
    "    name = random.choice(list(distributions.keys()))\n",
    "    distribution = distributions.pop(name)\n",
    "\n",
    "    # Generate n draws from the distribution\n",
    "    data = distribution.rvs(n)\n",
    "\n",
    "    # Compute sample mean at each n\n",
    "    sample_mean = np.empty(n)\n",
    "    for i in range(n):\n",
    "        sample_mean[i] = np.mean(data[:i+1])\n",
    "\n",
    "    # Plot\n",
    "    ax.plot(list(range(n)), data, 'o', color='grey', alpha=0.5)\n",
    "    axlabel = '$\\\\bar X_n$ for $X_i \\sim$' + name\n",
    "    ax.plot(list(range(n)), sample_mean, 'g-', lw=3, alpha=0.6, label=axlabel)\n",
    "    m = distribution.mean()\n",
    "    ax.plot(list(range(n)), [m] * n, 'k--', lw=1.5, label='$\\mu$')\n",
    "    ax.vlines(list(range(n)), m, data, lw=0.2)\n",
    "    ax.legend(**legend_args, fontsize=12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The three distributions are chosen at random from a selection stored in the dictionary `distributions`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CLT\n",
    "\n",
    "\n",
    "<a id='index-7'></a>\n",
    "Next, we turn to the central limit theorem, which tells us about the distribution of the deviation between sample averages and population means."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statement of the Theorem\n",
    "\n",
    "The central limit theorem is one of the most remarkable results in all of mathematics.\n",
    "\n",
    "In the classical IID setting, it tells us the following:\n",
    "\n",
    "\n",
    "<a id='statement-clt'></a>\n",
    "If the sequence $ X_1, \\ldots, X_n $ is IID, with common mean\n",
    "$ \\mu $ and common variance $ \\sigma^2 \\in (0, \\infty) $, then\n",
    "\n",
    "\n",
    "<a id='equation-lln-clt'></a>\n",
    "$$\n",
    "\\sqrt{n} ( \\bar X_n - \\mu ) \\stackrel { d } {\\to} N(0, \\sigma^2)\n",
    "\\quad \\text{as} \\quad\n",
    "n \\to \\infty \\tag{5}\n",
    "$$\n",
    "\n",
    "Here $ \\stackrel { d } {\\to} N(0, \\sigma^2) $ indicates [convergence in distribution](https://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_distribution) to a centered (i.e, zero mean) normal with standard deviation $ \\sigma $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intuition\n",
    "\n",
    "\n",
    "<a id='index-8'></a>\n",
    "The striking implication of the CLT is that for **any** distribution with\n",
    "finite second moment, the simple operation of adding independent\n",
    "copies **always** leads to a Gaussian curve.\n",
    "\n",
    "A relatively simple proof of the central limit theorem can be obtained by\n",
    "working with characteristic functions (see, e.g., theorem 9.5.6 of [[Dud02]](https://python.quantecon.org/zreferences.html#dudley2002)).\n",
    "\n",
    "The proof is elegant but almost anticlimactic, and it provides surprisingly little intuition.\n",
    "\n",
    "In fact, all of the proofs of the CLT that we know are similar in this respect.\n",
    "\n",
    "Why does adding independent copies produce a bell-shaped distribution?\n",
    "\n",
    "Part of the answer can be obtained by investigating the addition of independent Bernoulli\n",
    "random variables.\n",
    "\n",
    "In particular, let $ X_i $ be binary, with $ \\mathbb P\\{X_i = 0\\} = \\mathbb P\\{X_i =\n",
    "1 \\} = 0.5 $, and let $ X_1, \\ldots, X_n $ be independent.\n",
    "\n",
    "Think of $ X_i = 1 $ as a “success”, so that $ Y_n = \\sum_{i=1}^n X_i $ is the number of successes in $ n $ trials.\n",
    "\n",
    "The next figure plots the probability mass function of $ Y_n $ for $ n = 1, 2, 4, 8 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(10, 6))\n",
    "plt.subplots_adjust(hspace=0.4)\n",
    "axes = axes.flatten()\n",
    "ns = [1, 2, 4, 8]\n",
    "dom = list(range(9))\n",
    "\n",
    "for ax, n in zip(axes, ns):\n",
    "    b = binom(n, 0.5)\n",
    "    ax.bar(dom, b.pmf(dom), alpha=0.6, align='center')\n",
    "    ax.set(xlim=(-0.5, 8.5), ylim=(0, 0.55),\n",
    "           xticks=list(range(9)), yticks=(0, 0.2, 0.4),\n",
    "           title=f'$n = {n}$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When $ n = 1 $, the distribution is flat — one success or no successes\n",
    "have the same probability.\n",
    "\n",
    "When $ n = 2 $ we can either have 0, 1 or 2 successes.\n",
    "\n",
    "Notice the peak in probability mass at the mid-point $ k=1 $.\n",
    "\n",
    "The reason is that there are more ways to get 1 success (“fail then succeed”\n",
    "or “succeed then fail”) than to get zero or two successes.\n",
    "\n",
    "Moreover, the two trials are independent, so the outcomes “fail then succeed” and “succeed then\n",
    "fail” are just as likely as the outcomes “fail then fail” and “succeed then succeed”.\n",
    "\n",
    "(If there was positive correlation, say, then “succeed then fail” would be less likely than “succeed then succeed”)\n",
    "\n",
    "Here, already we have the essence of the CLT: addition under independence leads probability mass to pile up in the middle and thin out at the tails.\n",
    "\n",
    "For $ n = 4 $ and $ n = 8 $ we again get a peak at the “middle” value (halfway between the minimum and the maximum possible value).\n",
    "\n",
    "The intuition is the same — there are simply more ways to get these middle outcomes.\n",
    "\n",
    "If we continue, the bell-shaped curve becomes even more pronounced.\n",
    "\n",
    "We are witnessing the [binomial approximation of the normal distribution](https://en.wikipedia.org/wiki/De_Moivre%E2%80%93Laplace_theorem)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation 1\n",
    "\n",
    "Since the CLT seems almost magical, running simulations that verify its implications is one good way to build intuition.\n",
    "\n",
    "To this end, we now perform the following simulation\n",
    "\n",
    "1. Choose an arbitrary distribution $ F $ for the underlying observations $ X_i $.  \n",
    "1. Generate independent draws of $ Y_n := \\sqrt{n} ( \\bar X_n - \\mu ) $.  \n",
    "1. Use these draws to compute some measure of their distribution — such as a histogram.  \n",
    "1. Compare the latter to $ N(0, \\sigma^2) $.  \n",
    "\n",
    "\n",
    "Here’s some code that does exactly this for the exponential distribution\n",
    "$ F(x) = 1 - e^{- \\lambda x} $.\n",
    "\n",
    "(Please experiment with other choices of $ F $, but remember that, to conform with the conditions of the CLT, the distribution must have a finite second moment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Set parameters\n",
    "n = 250                  # Choice of n\n",
    "k = 100000               # Number of draws of Y_n\n",
    "distribution = expon(2)  # Exponential distribution, λ = 1/2\n",
    "μ, s = distribution.mean(), distribution.std()\n",
    "\n",
    "# Draw underlying RVs. Each row contains a draw of X_1,..,X_n\n",
    "data = distribution.rvs((k, n))\n",
    "# Compute mean of each row, producing k draws of \\bar X_n\n",
    "sample_means = data.mean(axis=1)\n",
    "# Generate observations of Y_n\n",
    "Y = np.sqrt(n) * (sample_means - μ)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "xmin, xmax = -3 * s, 3 * s\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.hist(Y, bins=60, alpha=0.5, density=True)\n",
    "xgrid = np.linspace(xmin, xmax, 200)\n",
    "ax.plot(xgrid, norm.pdf(xgrid, scale=s), 'k-', lw=2, label='$N(0, \\sigma^2)$')\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the absence of for loops — every operation is vectorized, meaning that the major calculations are all shifted to highly optimized C code.\n",
    "\n",
    "The fit to the normal density is already tight and can be further improved by increasing `n`.\n",
    "\n",
    "You can also experiment with other specifications of $ F $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation 2\n",
    "\n",
    "Our next simulation is somewhat like the first, except that we aim to track the distribution of $ Y_n := \\sqrt{n} ( \\bar X_n - \\mu ) $ as $ n $ increases.\n",
    "\n",
    "In the simulation, we’ll be working with random variables having $ \\mu = 0 $.\n",
    "\n",
    "Thus, when $ n=1 $, we have $ Y_1 = X_1 $, so the first distribution is just\n",
    "the distribution of the underlying random variable.\n",
    "\n",
    "For $ n=2 $, the distribution of $ Y_2 $ is that of $ (X_1 + X_2) / \\sqrt{2} $, and so on.\n",
    "\n",
    "What we expect is that, regardless of the distribution of the underlying\n",
    "random variable, the distribution of $ Y_n $ will smooth out into a bell-shaped curve.\n",
    "\n",
    "The next figure shows this process for $ X_i \\sim f $, where $ f $ was\n",
    "specified as the convex combination of three different beta densities.\n",
    "\n",
    "(Taking a convex combination is an easy way to produce an irregular shape for $ f $)\n",
    "\n",
    "In the figure, the closest density is that of $ Y_1 $, while the furthest is that of\n",
    "$ Y_5 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "beta_dist = beta(2, 2)\n",
    "\n",
    "def gen_x_draws(k):\n",
    "    \"\"\"\n",
    "    Returns a flat array containing k independent draws from the\n",
    "    distribution of X, the underlying random variable.  This distribution\n",
    "    is itself a convex combination of three beta distributions.\n",
    "    \"\"\"\n",
    "    bdraws = beta_dist.rvs((3, k))\n",
    "    # Transform rows, so each represents a different distribution\n",
    "    bdraws[0, :] -= 0.5\n",
    "    bdraws[1, :] += 0.6\n",
    "    bdraws[2, :] -= 1.1\n",
    "    # Set X[i] = bdraws[j, i], where j is a random draw from {0, 1, 2}\n",
    "    js = np.random.randint(0, 2, size=k)\n",
    "    X = bdraws[js, np.arange(k)]\n",
    "    # Rescale, so that the random variable is zero mean\n",
    "    m, sigma = X.mean(), X.std()\n",
    "    return (X - m) / sigma\n",
    "\n",
    "nmax = 5\n",
    "reps = 100000\n",
    "ns = list(range(1, nmax + 1))\n",
    "\n",
    "# Form a matrix Z such that each column is reps independent draws of X\n",
    "Z = np.empty((reps, nmax))\n",
    "for i in range(nmax):\n",
    "    Z[:, i] = gen_x_draws(reps)\n",
    "# Take cumulative sum across columns\n",
    "S = Z.cumsum(axis=1)\n",
    "# Multiply j-th column by sqrt j\n",
    "Y = (1 / np.sqrt(ns)) * S\n",
    "\n",
    "# Plot\n",
    "fig = plt.figure(figsize = (10, 6))\n",
    "ax = fig.gca(projection='3d')\n",
    "\n",
    "a, b = -3, 3\n",
    "gs = 100\n",
    "xs = np.linspace(a, b, gs)\n",
    "\n",
    "# Build verts\n",
    "greys = np.linspace(0.3, 0.7, nmax)\n",
    "verts = []\n",
    "for n in ns:\n",
    "    density = gaussian_kde(Y[:, n-1])\n",
    "    ys = density(xs)\n",
    "    verts.append(list(zip(xs, ys)))\n",
    "\n",
    "poly = PolyCollection(verts, facecolors=[str(g) for g in greys])\n",
    "poly.set_alpha(0.85)\n",
    "ax.add_collection3d(poly, zs=ns, zdir='x')\n",
    "\n",
    "ax.set(xlim3d=(1, nmax), xticks=(ns), ylabel='$Y_n$', zlabel='$p(y_n)$',\n",
    "       xlabel=(\"n\"), yticks=((-3, 0, 3)), ylim3d=(a, b),\n",
    "       zlim3d=(0, 0.4), zticks=((0.2, 0.4)))\n",
    "ax.invert_xaxis()\n",
    "# Rotates the plot 30 deg on z axis and 45 deg on x axis\n",
    "ax.view_init(30, 45)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the distribution smooths out into a bell curve as $ n $\n",
    "increases.\n",
    "\n",
    "We leave you to investigate its contents if you wish to know more.\n",
    "\n",
    "If you run the file from the ordinary IPython shell, the figure should pop up in a\n",
    "window that you can rotate with your mouse, giving different views on the\n",
    "density sequence.\n",
    "\n",
    "\n",
    "<a id='multivariate-clt'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Multivariate Case\n",
    "\n",
    "\n",
    "<a id='index-10'></a>\n",
    "The law of large numbers and central limit theorem work just as nicely in multidimensional settings.\n",
    "\n",
    "To state the results, let’s recall some elementary facts about random vectors.\n",
    "\n",
    "A random vector $ \\mathbf X $ is just a sequence of $ k $ random variables $ (X_1, \\ldots, X_k) $.\n",
    "\n",
    "Each realization of $ \\mathbf X $ is an element of $ \\mathbb R^k $.\n",
    "\n",
    "A collection of random vectors $ \\mathbf X_1, \\ldots, \\mathbf X_n $ is called independent if, given any $ n $ vectors $ \\mathbf x_1, \\ldots, \\mathbf x_n $ in $ \\mathbb R^k $, we have\n",
    "\n",
    "$$\n",
    "\\mathbb P\\{\\mathbf X_1 \\leq \\mathbf x_1,\\ldots, \\mathbf X_n \\leq \\mathbf x_n \\}\n",
    "= \\mathbb P\\{\\mathbf X_1 \\leq \\mathbf x_1 \\}\n",
    "\\times \\cdots \\times \\mathbb P\\{ \\mathbf X_n \\leq \\mathbf x_n \\}\n",
    "$$\n",
    "\n",
    "(The vector inequality $ \\mathbf X \\leq \\mathbf x $ means that $ X_j \\leq x_j $ for $ j = 1,\\ldots,k $)\n",
    "\n",
    "Let $ \\mu_j := \\mathbb E [X_j] $ for all $ j =1,\\ldots,k $.\n",
    "\n",
    "The expectation $ \\mathbb E [\\mathbf X] $ of $ \\mathbf X $ is defined to be the vector of expectations:\n",
    "\n",
    "$$\n",
    "\\mathbb E [\\mathbf X] :=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "    \\mathbb E [X_1] \\\\\n",
    "    \\mathbb E [X_2] \\\\\n",
    "    \\vdots \\\\\n",
    "    \\mathbb E [X_k]\n",
    "\\end{array}\n",
    "\\right) =\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "    \\mu_1 \\\\\n",
    "    \\mu_2\\\\\n",
    "    \\vdots \\\\\n",
    "    \\mu_k\n",
    "\\end{array}\n",
    "\\right) =: \\boldsymbol \\mu\n",
    "$$\n",
    "\n",
    "The *variance-covariance matrix* of random vector $ \\mathbf X $ is defined as\n",
    "\n",
    "$$\n",
    "\\mathop{\\mathrm{Var}}[\\mathbf X]\n",
    ":= \\mathbb E\n",
    "[ (\\mathbf X - \\boldsymbol \\mu) (\\mathbf X - \\boldsymbol \\mu)']\n",
    "$$\n",
    "\n",
    "Expanding this out, we get\n",
    "\n",
    "$$\n",
    "\\mathop{\\mathrm{Var}}[\\mathbf X] =\n",
    "\\left(\n",
    "\\begin{array}{ccc}\n",
    "    \\mathbb E [(X_1 - \\mu_1)(X_1 - \\mu_1)]\n",
    "        & \\cdots & \\mathbb E [(X_1 - \\mu_1)(X_k - \\mu_k)] \\\\\n",
    "    \\mathbb E [(X_2 - \\mu_2)(X_1 - \\mu_1)]\n",
    "        & \\cdots & \\mathbb E [(X_2 - \\mu_2)(X_k - \\mu_k)] \\\\\n",
    "    \\vdots & \\vdots & \\vdots \\\\\n",
    "    \\mathbb E [(X_k - \\mu_k)(X_1 - \\mu_1)]\n",
    "        & \\cdots & \\mathbb E [(X_k - \\mu_k)(X_k - \\mu_k)] \\\\\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "The $ j,k $-th term is the scalar covariance between $ X_j $ and\n",
    "$ X_k $.\n",
    "\n",
    "With this notation, we can proceed to the multivariate LLN and CLT.\n",
    "\n",
    "Let $ \\mathbf X_1, \\ldots, \\mathbf X_n $ be a sequence of independent and\n",
    "identically distributed random vectors, each one taking values in\n",
    "$ \\mathbb R^k $.\n",
    "\n",
    "Let $ \\boldsymbol \\mu $ be the vector $ \\mathbb E [\\mathbf X_i] $, and let $ \\Sigma $\n",
    "be the variance-covariance matrix of $ \\mathbf X_i $.\n",
    "\n",
    "Interpreting vector addition and scalar multiplication in the usual way (i.e., pointwise), let\n",
    "\n",
    "$$\n",
    "\\bar{\\mathbf X}_n := \\frac{1}{n} \\sum_{i=1}^n \\mathbf X_i\n",
    "$$\n",
    "\n",
    "In this setting, the LLN tells us that\n",
    "\n",
    "\n",
    "<a id='equation-lln-asmv'></a>\n",
    "$$\n",
    "\\mathbb P \\left\\{ \\bar{\\mathbf X}_n \\to \\boldsymbol \\mu \\text{ as } n \\to \\infty \\right\\} = 1 \\tag{6}\n",
    "$$\n",
    "\n",
    "Here $ \\bar{\\mathbf X}_n \\to \\boldsymbol \\mu $ means that $ \\| \\bar{\\mathbf X}_n - \\boldsymbol \\mu \\| \\to 0 $, where $ \\| \\cdot \\| $ is the standard Euclidean norm.\n",
    "\n",
    "The CLT tells us that, provided $ \\Sigma $ is finite,\n",
    "\n",
    "\n",
    "<a id='equation-lln-cltmv'></a>\n",
    "$$\n",
    "\\sqrt{n} ( \\bar{\\mathbf X}_n - \\boldsymbol \\mu ) \\stackrel { d } {\\to} N(\\mathbf 0, \\Sigma)\n",
    "\\quad \\text{as} \\quad\n",
    "n \\to \\infty \\tag{7}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "\n",
    "<a id='lln-ex1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "One very useful consequence of the central limit theorem is as follows.\n",
    "\n",
    "Assume the conditions of the CLT as [stated above](#statement-clt).\n",
    "\n",
    "If $ g \\colon \\mathbb R \\to \\mathbb R $ is differentiable at $ \\mu $ and $ g'(\\mu) \\not= 0 $, then\n",
    "\n",
    "\n",
    "<a id='equation-lln-dm'></a>\n",
    "$$\n",
    "\\sqrt{n} \\{ g(\\bar X_n) - g(\\mu) \\}\n",
    "\\stackrel { d } {\\to} N(0, g'(\\mu)^2 \\sigma^2)\n",
    "\\quad \\text{as} \\quad\n",
    "n \\to \\infty \\tag{8}\n",
    "$$\n",
    "\n",
    "This theorem is used frequently in statistics to obtain the asymptotic distribution of estimators — many of which can be expressed as functions of sample means.\n",
    "\n",
    "(These kinds of results are often said to use the “delta method”)\n",
    "\n",
    "The proof is based on a Taylor expansion of $ g $ around the point $ \\mu $.\n",
    "\n",
    "Taking the result as given, let the distribution $ F $ of each $ X_i $ be uniform on $ [0, \\pi / 2] $ and let $ g(x) = \\sin(x) $.\n",
    "\n",
    "Derive the asymptotic distribution of $ \\sqrt{n} \\{ g(\\bar X_n) - g(\\mu) \\} $ and illustrate convergence in the same spirit as the program `illustrate_clt.py` discussed above.\n",
    "\n",
    "What happens when you replace $ [0, \\pi / 2] $ with $ [0, \\pi] $?\n",
    "\n",
    "What is the source of the problem?\n",
    "\n",
    "\n",
    "<a id='lln-ex2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "Here’s a result that’s often used in developing statistical tests, and is connected to the multivariate central limit theorem.\n",
    "\n",
    "If you study econometric theory, you will see this result used again and again.\n",
    "\n",
    "Assume the setting of the multivariate CLT [discussed above](#multivariate-clt), so that\n",
    "\n",
    "1. $ \\mathbf X_1, \\ldots, \\mathbf X_n $ is a sequence of IID random vectors, each taking values in $ \\mathbb R^k $.  \n",
    "1. $ \\boldsymbol \\mu := \\mathbb E [\\mathbf X_i] $, and $ \\Sigma $ is the variance-covariance matrix of $ \\mathbf X_i $.  \n",
    "1. The convergence  \n",
    "\n",
    "\n",
    "\n",
    "<a id='equation-lln-cltmv2'></a>\n",
    "$$\n",
    "\\sqrt{n} ( \\bar{\\mathbf X}_n - \\boldsymbol \\mu ) \\stackrel { d } {\\to} N(\\mathbf 0, \\Sigma) \\tag{9}\n",
    "$$\n",
    "\n",
    "is valid.\n",
    "\n",
    "In a statistical setting, one often wants the right-hand side to be **standard** normal so that confidence intervals are easily computed.\n",
    "\n",
    "This normalization can be achieved on the basis of three observations.\n",
    "\n",
    "First, if $ \\mathbf X $ is a random vector in $ \\mathbb R^k $ and $ \\mathbf A $ is constant and $ k \\times k $, then\n",
    "\n",
    "$$\n",
    "\\mathop{\\mathrm{Var}}[\\mathbf A \\mathbf X]\n",
    "= \\mathbf A \\mathop{\\mathrm{Var}}[\\mathbf X] \\mathbf A'\n",
    "$$\n",
    "\n",
    "Second, by the [continuous mapping theorem](https://en.wikipedia.org/wiki/Continuous_mapping_theorem), if $ \\mathbf Z_n \\stackrel{d}{\\to} \\mathbf Z $ in $ \\mathbb R^k $ and $ \\mathbf A $ is constant and $ k \\times k $, then\n",
    "\n",
    "$$\n",
    "\\mathbf A \\mathbf Z_n\n",
    "\\stackrel{d}{\\to} \\mathbf A \\mathbf Z\n",
    "$$\n",
    "\n",
    "Third, if $ \\mathbf S $ is a $ k \\times k $ symmetric positive definite matrix, then there\n",
    "exists a symmetric positive definite matrix $ \\mathbf Q $, called the inverse\n",
    "[square root](https://en.wikipedia.org/wiki/Square_root_of_a_matrix) of $ \\mathbf S $, such that\n",
    "\n",
    "$$\n",
    "\\mathbf Q \\mathbf S\\mathbf Q' = \\mathbf I\n",
    "$$\n",
    "\n",
    "Here $ \\mathbf I $ is the $ k \\times k $ identity matrix.\n",
    "\n",
    "Putting these things together, your first exercise is to show that if\n",
    "$ \\mathbf Q $ is the inverse square root of $ \\mathbf \\Sigma $, then\n",
    "\n",
    "$$\n",
    "\\mathbf Z_n := \\sqrt{n} \\mathbf Q ( \\bar{\\mathbf X}_n - \\boldsymbol \\mu )\n",
    "\\stackrel{d}{\\to}\n",
    "\\mathbf Z \\sim N(\\mathbf 0, \\mathbf I)\n",
    "$$\n",
    "\n",
    "Applying the continuous mapping theorem one more time tells us that\n",
    "\n",
    "$$\n",
    "\\| \\mathbf Z_n \\|^2\n",
    "\\stackrel{d}{\\to}\n",
    "\\| \\mathbf Z \\|^2\n",
    "$$\n",
    "\n",
    "Given the distribution of $ \\mathbf Z $, we conclude that\n",
    "\n",
    "\n",
    "<a id='equation-lln-ctc'></a>\n",
    "$$\n",
    "n \\| \\mathbf Q ( \\bar{\\mathbf X}_n - \\boldsymbol \\mu ) \\|^2\n",
    "\\stackrel{d}{\\to}\n",
    "\\chi^2(k) \\tag{10}\n",
    "$$\n",
    "\n",
    "where $ \\chi^2(k) $ is the chi-squared distribution with $ k $ degrees\n",
    "of freedom.\n",
    "\n",
    "(Recall that $ k $ is the dimension of $ \\mathbf X_i $, the underlying random vectors)\n",
    "\n",
    "Your second exercise is to illustrate the convergence in [(10)](#equation-lln-ctc) with a simulation.\n",
    "\n",
    "In doing so, let\n",
    "\n",
    "$$\n",
    "\\mathbf X_i :=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "    W_i \\\\\n",
    "    U_i + W_i\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- each $ W_i $ is an IID draw from the uniform distribution on $ [-1, 1] $.  \n",
    "- each $ U_i $ is an IID draw from the uniform distribution on $ [-2, 2] $.  \n",
    "- $ U_i $ and $ W_i $ are independent of each other.  \n",
    "\n",
    "\n",
    "Hints:\n",
    "\n",
    "1. `scipy.linalg.sqrtm(A)` computes the square root of `A`.  You still need to invert it.  \n",
    "1. You should be able to work out $ \\Sigma $ from the preceding information.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Here is one solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Illustrates the delta method, a consequence of the central limit theorem.\n",
    "\"\"\"\n",
    "\n",
    "# Set parameters\n",
    "n = 250\n",
    "replications = 100000\n",
    "distribution = uniform(loc=0, scale=(np.pi / 2))\n",
    "μ, s = distribution.mean(), distribution.std()\n",
    "\n",
    "g = np.sin\n",
    "g_prime = np.cos\n",
    "\n",
    "# Generate obs of sqrt{n} (g(X_n) - g(μ))\n",
    "data = distribution.rvs((replications, n))\n",
    "sample_means = data.mean(axis=1)  # Compute mean of each row\n",
    "error_obs = np.sqrt(n) * (g(sample_means) - g(μ))\n",
    "\n",
    "# Plot\n",
    "asymptotic_sd = g_prime(μ) * s\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "xmin = -3 * g_prime(μ) * s\n",
    "xmax = -xmin\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.hist(error_obs, bins=60, alpha=0.5, density=True)\n",
    "xgrid = np.linspace(xmin, xmax, 200)\n",
    "lb = \"$N(0, g'(\\mu)^2  \\sigma^2)$\"\n",
    "ax.plot(xgrid, norm.pdf(xgrid, scale=asymptotic_sd), 'k-', lw=2, label=lb)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens when you replace $ [0, \\pi / 2] $ with\n",
    "$ [0, \\pi] $?\n",
    "\n",
    "In this case, the mean $ \\mu $ of this distribution is\n",
    "$ \\pi/2 $, and since $ g' = \\cos $, we have $ g'(\\mu) = 0 $.\n",
    "\n",
    "Hence the conditions of the delta theorem are not satisfied."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "First we want to verify the claim that\n",
    "\n",
    "$$\n",
    "\\sqrt{n} \\mathbf Q ( \\bar{\\mathbf X}_n - \\boldsymbol \\mu )\n",
    "\\stackrel{d}{\\to}\n",
    "N(\\mathbf 0, \\mathbf I)\n",
    "$$\n",
    "\n",
    "This is straightforward given the facts presented in the exercise.\n",
    "\n",
    "Let\n",
    "\n",
    "$$\n",
    "\\mathbf Y_n := \\sqrt{n} ( \\bar{\\mathbf X}_n - \\boldsymbol \\mu )\n",
    "\\quad \\text{and} \\quad\n",
    "\\mathbf Y \\sim N(\\mathbf 0, \\Sigma)\n",
    "$$\n",
    "\n",
    "By the multivariate CLT and the continuous mapping theorem, we have\n",
    "\n",
    "$$\n",
    "\\mathbf Q \\mathbf Y_n\n",
    "\\stackrel{d}{\\to}\n",
    "\\mathbf Q \\mathbf Y\n",
    "$$\n",
    "\n",
    "Since linear combinations of normal random variables are normal, the\n",
    "vector $ \\mathbf Q \\mathbf Y $ is also normal.\n",
    "\n",
    "Its mean is clearly $ \\mathbf 0 $, and its variance-covariance\n",
    "matrix is\n",
    "\n",
    "$$\n",
    "\\mathrm{Var}[\\mathbf Q \\mathbf Y]\n",
    "= \\mathbf Q \\mathrm{Var}[\\mathbf Y] \\mathbf Q'\n",
    "= \\mathbf Q \\Sigma \\mathbf Q'\n",
    "= \\mathbf I\n",
    "$$\n",
    "\n",
    "In conclusion,\n",
    "$ \\mathbf Q \\mathbf Y_n \\stackrel{d}{\\to} \\mathbf Q \\mathbf Y \\sim N(\\mathbf 0, \\mathbf I) $,\n",
    "which is what we aimed to show.\n",
    "\n",
    "Now we turn to the simulation exercise.\n",
    "\n",
    "Our solution is as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Set parameters\n",
    "n = 250\n",
    "replications = 50000\n",
    "dw = uniform(loc=-1, scale=2)  # Uniform(-1, 1)\n",
    "du = uniform(loc=-2, scale=4)  # Uniform(-2, 2)\n",
    "sw, su = dw.std(), du.std()\n",
    "vw, vu = sw**2, su**2\n",
    "Σ = ((vw, vw), (vw, vw + vu))\n",
    "Σ = np.array(Σ)\n",
    "\n",
    "# Compute Σ^{-1/2}\n",
    "Q = inv(sqrtm(Σ))\n",
    "\n",
    "# Generate observations of the normalized sample mean\n",
    "error_obs = np.empty((2, replications))\n",
    "for i in range(replications):\n",
    "    # Generate one sequence of bivariate shocks\n",
    "    X = np.empty((2, n))\n",
    "    W = dw.rvs(n)\n",
    "    U = du.rvs(n)\n",
    "    # Construct the n observations of the random vector\n",
    "    X[0, :] = W\n",
    "    X[1, :] = W + U\n",
    "    # Construct the i-th observation of Y_n\n",
    "    error_obs[:, i] = np.sqrt(n) * X.mean(axis=1)\n",
    "\n",
    "# Premultiply by Q and then take the squared norm\n",
    "temp = Q @ error_obs\n",
    "chisq_obs = np.sum(temp**2, axis=0)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "xmax = 8\n",
    "ax.set_xlim(0, xmax)\n",
    "xgrid = np.linspace(0, xmax, 200)\n",
    "lb = \"Chi-squared with 2 degrees of freedom\"\n",
    "ax.plot(xgrid, chi2.pdf(xgrid, 2), 'k-', lw=2, label=lb)\n",
    "ax.legend()\n",
    "ax.hist(chisq_obs, bins=50, density=True)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "date": 1584334743.4384282,
  "filename": "lln_clt.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "LLN and CLT"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}