{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "77bf0a64",
   "metadata": {},
   "source": [
    "# LLN and CLT"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e39a6f33",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture illustrates two of the most important results in probability and statistics:\n",
    "\n",
    "1. the law of large numbers (LLN) and  \n",
    "1. the central limit theorem (CLT).  \n",
    "\n",
    "\n",
    "These beautiful theorems lie behind many of the most fundamental results in\n",
    "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\n",
    "are based on do not hold.\n",
    "\n",
    "This lecture will focus on the univariate case (the multivariate case is treated [in a more advanced lecture](https://python.quantecon.org/lln_clt.html#the-multivariate-case)).\n",
    "\n",
    "We’ll need the following imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2c83c73",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats as st"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ca8baa7",
   "metadata": {},
   "source": [
    "\n",
    "<a id='lln-mr'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa7d8b3a",
   "metadata": {},
   "source": [
    "## The law of large numbers\n",
    "\n",
    "\n",
    "<a id='index-0'></a>\n",
    "We begin with the law of large numbers, which tells us when sample averages\n",
    "will converge to their population means."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d71dc73",
   "metadata": {},
   "source": [
    "### The LLN in action\n",
    "\n",
    "Let’s see an example of the LLN in action before we go further."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f3e2c54",
   "metadata": {},
   "source": [
    "### \n",
    "\n",
    "Consider a [Bernoulli random variable](https://en.wikipedia.org/wiki/Bernoulli_distribution) $ X $ with parameter $ p $.\n",
    "\n",
    "This means that $ X $ takes values in $ \\{0,1\\} $ and $ \\mathbb P\\{X=1\\} = p $.\n",
    "\n",
    "We can think of drawing $ X $ as tossing a biased coin where\n",
    "\n",
    "- the coin falls on “heads” with probability $ p $ and  \n",
    "- the coin falls on “tails” with probability $ 1-p $  \n",
    "\n",
    "\n",
    "We set $ X=1 $ if the coin is “heads” and zero otherwise.\n",
    "\n",
    "The (population) mean of $ X $ is\n",
    "\n",
    "$$\n",
    "\\mathbb E X \n",
    "    = 0 \\cdot \\mathbb P\\{X=0\\} + 1 \\cdot \\mathbb P\\{X=1\\} = \\mathbb P\\{X=1\\} = p\n",
    "$$\n",
    "\n",
    "We can generate a draw of $ X $ with `scipy.stats` (imported as `st`) as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d666e2e3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "p = 0.8\n",
    "X = st.bernoulli.rvs(p)\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37cd8195",
   "metadata": {},
   "source": [
    "In this setting, the LLN tells us if we flip the coin many times, the fraction\n",
    "of heads that we see will be close to the mean $ p $.\n",
    "\n",
    "We use $ n $ to represent the number of times the coin is flipped.\n",
    "\n",
    "Let’s check this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4a2495b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n = 1_000_000\n",
    "X_draws = st.bernoulli.rvs(p, size=n)\n",
    "print(X_draws.mean()) # count the number of 1's and divide by n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e08caa76",
   "metadata": {},
   "source": [
    "If we change $ p $ the claim still holds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7ee2ddf",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "p = 0.3\n",
    "X_draws = st.bernoulli.rvs(p, size=n)\n",
    "print(X_draws.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbdaa771",
   "metadata": {},
   "source": [
    "Let’s connect this to the discussion above, where we said the sample average\n",
    "converges to the “population mean”.\n",
    "\n",
    "Think of $ X_1, \\ldots, X_n $ as independent flips of the coin.\n",
    "\n",
    "The population mean is the mean in an infinite sample, which equals the\n",
    "expectation $ \\mathbb E X $.\n",
    "\n",
    "The sample mean of the draws $ X_1, \\ldots, X_n $ is\n",
    "\n",
    "$$\n",
    "\\bar X_n := \\frac{1}{n} \\sum_{i=1}^n X_i\n",
    "$$\n",
    "\n",
    "In this case, it is the fraction of draws that equal one (the number of heads divided by $ n $).\n",
    "\n",
    "Thus, the LLN tells us that for the Bernoulli trials above\n",
    "\n",
    "\n",
    "<a id='equation-exp'></a>\n",
    "$$\n",
    "\\bar X_n \\to \\mathbb E X = p\n",
    "    \\qquad (n \\to \\infty) \\tag{20.1}\n",
    "$$\n",
    "\n",
    "This is exactly what we illustrated in the code.\n",
    "\n",
    "\n",
    "<a id='lln-ksl'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3590674f",
   "metadata": {},
   "source": [
    "### Statement of the LLN\n",
    "\n",
    "Let’s state the LLN more carefully.\n",
    "\n",
    "Let $ X_1, \\ldots, X_n $ be random variables, all of which have the same\n",
    "distribution.\n",
    "\n",
    "These random variables can be continuous or discrete.\n",
    "\n",
    "For simplicity we will\n",
    "\n",
    "- assume they are continuous and  \n",
    "- let $ f $ denote their common density function  \n",
    "\n",
    "\n",
    "The last statement means that for any $ i $ in $ \\{1, \\ldots, n\\} $ and any\n",
    "numbers $ a, b $,\n",
    "\n",
    "$$\n",
    "\\mathbb P\\{a \\leq X_i \\leq b\\} = \\int_a^b f(x) dx\n",
    "$$\n",
    "\n",
    "(For the discrete case, we need to replace densities with probability mass\n",
    "functions and integrals with sums.)\n",
    "\n",
    "Let $ \\mu $ denote the common mean of this sample.\n",
    "\n",
    "Thus, for each $ i $,\n",
    "\n",
    "$$\n",
    "\\mu := \\mathbb E X_i = \\int_{-\\infty}^{\\infty} x f(x) dx\n",
    "$$\n",
    "\n",
    "The sample mean is\n",
    "\n",
    "$$\n",
    "\\bar X_n := \\frac{1}{n} \\sum_{i=1}^n X_i\n",
    "$$\n",
    "\n",
    "The next theorem is called Kolmogorov’s strong law of large numbers.\n",
    "\n",
    "\n",
    "<a id='iid-theorem'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48b55ffc",
   "metadata": {},
   "source": [
    "### \n",
    "\n",
    "If $ X_1, \\ldots, X_n $ are IID and $ \\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{20.2}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- IID means independent and identically distributed and  \n",
    "- $ \\mathbb E |X| = \\int_{-\\infty}^\\infty |x| f(x) dx $  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "273a81ce",
   "metadata": {},
   "source": [
    "### Comments on the theorem\n",
    "\n",
    "What does the probability one statement in the theorem 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 (although this [isn’t strictly true](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, [(20.2)](#equation-lln-as) should be interpreted as meaning that the\n",
    "probability of the computer producing a sequence where $ \\bar X_n \\to \\mu $\n",
    "fails to occur is zero."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b25c57d",
   "metadata": {},
   "source": [
    "### Illustration\n",
    "\n",
    "\n",
    "<a id='index-1'></a>\n",
    "Let’s illustrate the LLN using simulation.\n",
    "\n",
    "When we illustrate it, we will use a key idea: the sample mean $ \\bar X_n $ is\n",
    "itself a random variable.\n",
    "\n",
    "The reason $ \\bar X_n $ is a random variable is that it’s a function of the\n",
    "random variables $ X_1, \\ldots, X_n $.\n",
    "\n",
    "What we are going to do now is\n",
    "\n",
    "1. pick some fixed distribution to draw each $ X_i $ from  \n",
    "1. set $ n $ to some large number  \n",
    "\n",
    "\n",
    "and then repeat the following three instructions.\n",
    "\n",
    "1. generate the draws $ X_1, \\ldots, X_n $  \n",
    "1. calculate the sample mean $ \\bar X_n $ and record its value in an array `sample_means`  \n",
    "1. go to step 1.  \n",
    "\n",
    "\n",
    "We will loop over these three steps $ m $ times, where $ m $ is some large integer.\n",
    "\n",
    "The array `sample_means` will now contain $ m $ draws of the random variable $ \\bar X_n $.\n",
    "\n",
    "If we histogram these observations of $ \\bar X_n $, we should see that they are clustered around the population mean $ \\mathbb E X $.\n",
    "\n",
    "Moreover, if we repeat the exercise with a larger value of $ n $, we should see that the observations are even more tightly clustered around the population mean.\n",
    "\n",
    "This is, in essence, what the LLN is telling us.\n",
    "\n",
    "To implement these steps, we will use functions.\n",
    "\n",
    "Our first function generates a sample mean of size $ n $ given a distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a9b9915",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def draw_means(X_distribution,  # The distribution of each X_i\n",
    "               n):              # The size of the sample mean\n",
    "\n",
    "    # Generate n draws: X_1, ..., X_n\n",
    "    X_samples = X_distribution.rvs(size=n)\n",
    "\n",
    "    # Return the sample mean\n",
    "    return np.mean(X_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b831740e",
   "metadata": {},
   "source": [
    "Now we write a function to generate $ m $ sample means and histogram them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e01ac3ed",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def generate_histogram(X_distribution, n, m): \n",
    "\n",
    "    # Compute m sample means\n",
    "\n",
    "    sample_means = np.empty(m)\n",
    "    for j in range(m):\n",
    "      sample_means[j] = draw_means(X_distribution, n) \n",
    "\n",
    "    # Generate a histogram\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.hist(sample_means, bins=30, alpha=0.5, density=True)\n",
    "    μ = X_distribution.mean()  # Get the population mean\n",
    "    σ = X_distribution.std()    # and the standard deviation\n",
    "    ax.axvline(x=μ, ls=\"--\", c=\"k\", label=fr\"$\\mu = {μ}$\")\n",
    "     \n",
    "    ax.set_xlim(μ - σ, μ + σ)\n",
    "    ax.set_xlabel(r'$\\bar X_n$', size=12)\n",
    "    ax.set_ylabel('density', size=12)\n",
    "    ax.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c261309",
   "metadata": {},
   "source": [
    "Now we call the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "715b8551",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# pick a distribution to draw each $X_i$ from\n",
    "X_distribution = st.norm(loc=5, scale=2) \n",
    "# Call the function\n",
    "generate_histogram(X_distribution, n=1_000, m=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f094d2b2",
   "metadata": {},
   "source": [
    "We can see that the distribution of $ \\bar X $ is clustered around $ \\mathbb E X $\n",
    "as expected.\n",
    "\n",
    "Let’s vary `n` to see how the distribution of the sample mean changes.\n",
    "\n",
    "We will use a [violin plot](https://intro.quantecon.org/prob_dist.html#violin-plots) to show the different distributions.\n",
    "\n",
    "Each distribution in the violin plot represents the distribution of $ X_n $ for some $ n $, calculated by simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e1dbd72",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def means_violin_plot(distribution,  \n",
    "                      ns = [1_000, 10_000, 100_000],\n",
    "                      m = 10_000):\n",
    "\n",
    "    data = []\n",
    "    for n in ns:\n",
    "        sample_means = [draw_means(distribution, n) for i in range(m)]\n",
    "        data.append(sample_means)\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "\n",
    "    ax.violinplot(data)\n",
    "    μ = distribution.mean()\n",
    "    ax.axhline(y=μ, ls=\"--\", c=\"k\", label=fr\"$\\mu = {μ}$\")\n",
    "\n",
    "    labels=[fr'$n = {n}$' for n in ns]\n",
    "\n",
    "    ax.set_xticks(np.arange(1, len(labels) + 1), labels=labels)\n",
    "    ax.set_xlim(0.25, len(labels) + 0.75)\n",
    "\n",
    "\n",
    "    plt.subplots_adjust(bottom=0.15, wspace=0.05)\n",
    "\n",
    "    ax.set_ylabel('density', size=12)\n",
    "    ax.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b20b9f9c",
   "metadata": {},
   "source": [
    "Let’s try with a normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26d73777",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "means_violin_plot(st.norm(loc=5, scale=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77b1e3ff",
   "metadata": {},
   "source": [
    "As $ n $ gets large, more probability mass clusters around the population mean $ \\mu $.\n",
    "\n",
    "Now let’s try with a Beta distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43fd5e06",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "means_violin_plot(st.beta(6, 6))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da8f5981",
   "metadata": {},
   "source": [
    "We get a similar result."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a7a2d34",
   "metadata": {},
   "source": [
    "## Breaking the LLN\n",
    "\n",
    "We have to pay attention to the assumptions in the statement of the LLN.\n",
    "\n",
    "If these assumptions do not hold, then the LLN might fail."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcea7dd1",
   "metadata": {},
   "source": [
    "### Infinite first moment\n",
    "\n",
    "As indicated by the theorem, the LLN can break when $ \\mathbb E |X| $ is not finite.\n",
    "\n",
    "We can demonstrate this using the [Cauchy distribution](https://en.wikipedia.org/wiki/Cauchy_distribution).\n",
    "\n",
    "The Cauchy distribution has the following property:\n",
    "\n",
    "If $ X_1, \\ldots, X_n $ are IID and Cauchy, then so is $ \\bar X_n $.\n",
    "\n",
    "This means that the distribution of $ \\bar X_n $ does not eventually concentrate on a single number.\n",
    "\n",
    "Hence the LLN does not hold.\n",
    "\n",
    "The LLN fails to hold here because the assumption $ \\mathbb E|X| < \\infty $ is violated by the Cauchy distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b4e35bb",
   "metadata": {},
   "source": [
    "### Failure of the IID condition\n",
    "\n",
    "The LLN can also fail to hold when the IID assumption is violated."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "105e9b7b",
   "metadata": {},
   "source": [
    "### \n",
    "\n",
    "$$\n",
    "X_0 \\sim N(0,1)\n",
    "    \\quad \\text{and} \\quad\n",
    "    X_i = X_{i-1} \\quad \\text{for} \\quad i = 1, ..., n\n",
    "$$\n",
    "\n",
    "In this case,\n",
    "\n",
    "$$\n",
    "\\bar X_n = \\frac{1}{n} \\sum_{i=1}^n X_i = X_0 \\sim N(0,1)\n",
    "$$\n",
    "\n",
    "Therefore, the distribution of $ \\bar X_n $ is $ N(0,1) $ for all $ n $!\n",
    "\n",
    "Does this contradict the LLN, which says that the distribution of $ \\bar X_n $\n",
    "collapses to the single point $ \\mu $?\n",
    "\n",
    "No, the LLN is correct — the issue is that its assumptions are not\n",
    "satisfied.\n",
    "\n",
    "In particular, the sequence $ X_1, \\ldots, X_n $ is not independent.\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">Although in this case the violation of IID breaks the LLN, there *are* situations\n",
    "where IID fails but the LLN still holds.\n",
    "\n",
    "We will show an example in the [exercise](#lln_ex3)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81da6ca6",
   "metadata": {},
   "source": [
    "## Central limit theorem\n",
    "\n",
    "\n",
    "<a id='index-2'></a>\n",
    "Next, we turn to the central limit theorem (CLT), which tells us about the\n",
    "distribution of the deviation between sample averages and population means."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa726a17",
   "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 IID setting, it tells us the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b029f83",
   "metadata": {},
   "source": [
    "### \n",
    "\n",
    "If $ X_1, \\ldots, X_n $ is IID with common mean $ \\mu $ and common variance\n",
    "$ \\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{20.3}\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 $.\n",
    "\n",
    "The striking implication of the CLT is that for any distribution with\n",
    "finite [second moment](https://en.wikipedia.org/wiki/Moment_%28mathematics%29), the simple operation of adding independent\n",
    "copies always leads to a Gaussian(Normal) curve."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79185483",
   "metadata": {},
   "source": [
    "### Simulation 1\n",
    "\n",
    "Since the CLT seems almost magical, running simulations that verify its implications is one good way to build understanding.\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.)\n",
    "\n",
    "\n",
    "<a id='sim-one'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d198d46c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Set parameters\n",
    "n = 250         # Choice of n\n",
    "k = 1_000_000        # Number of draws of Y_n\n",
    "distribution = st.expon(2) # Exponential distribution, λ = 1/2\n",
    "μ, σ = 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 * σ, 3 * σ\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.hist(Y, bins=60, alpha=0.4, density=True)\n",
    "xgrid = np.linspace(xmin, xmax, 200)\n",
    "ax.plot(xgrid, st.norm.pdf(xgrid, scale=σ), \n",
    "        'k-', lw=2, label=r'$N(0, \\sigma^2)$')\n",
    "ax.set_xlabel(r\"$Y_n$\", size=12)\n",
    "ax.set_ylabel(r\"$density$\", size=12)\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be6be5ec",
   "metadata": {},
   "source": [
    "(Notice the absence of for loops — every operation is vectorized, meaning that the major calculations are all shifted to fast C code.)\n",
    "\n",
    "The fit to the normal density is already tight and can be further improved by increasing `n`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7f94868",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4492736",
   "metadata": {},
   "source": [
    "## Exercise 20.1\n",
    "\n",
    "Repeat the simulation [above](#sim-one) with the [Beta distribution](https://en.wikipedia.org/wiki/Beta_distribution).\n",
    "\n",
    "You can choose any $ \\alpha > 0 $ and $ \\beta > 0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "544adfb8",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 20.1](https://intro.quantecon.org/#lln_ex1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31ef1c1a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Set parameters\n",
    "n = 250         # Choice of n\n",
    "k = 1_000_000        # Number of draws of Y_n\n",
    "distribution = st.beta(2,2) # We chose Beta(2, 2) as an example\n",
    "μ, σ = 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 * σ, 3 * σ\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.hist(Y, bins=60, alpha=0.4, density=True)\n",
    "ax.set_xlabel(r\"$Y_n$\", size=12)\n",
    "ax.set_ylabel(r\"$density$\", size=12)\n",
    "xgrid = np.linspace(xmin, xmax, 200)\n",
    "ax.plot(xgrid, st.norm.pdf(xgrid, scale=σ), 'k-', lw=2, label=r'$N(0, \\sigma^2)$')\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "324bb512",
   "metadata": {},
   "source": [
    "## Exercise 20.2\n",
    "\n",
    "At the start of this lecture we discussed Bernoulli random variables.\n",
    "\n",
    "NumPy doesn’t provide a `bernoulli` function that we can sample from.\n",
    "\n",
    "However, we can generate a draw of Bernoulli $ X $ using NumPy via"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a4f5f87",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "U = np.random.rand()\n",
    "X = 1 if U < p else 0\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e212a177",
   "metadata": {},
   "source": [
    "Explain why this provides a random variable $ X $ with the right distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1df1d29",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 20.2](https://intro.quantecon.org/#lln_ex2)\n",
    "\n",
    "We can write $ X $ as $ X = \\mathbf 1\\{U < p\\} $ where $ \\mathbf 1 $ is the\n",
    "[indicator function](https://en.wikipedia.org/wiki/Indicator_function) (i.e.,\n",
    "1 if the statement is true and zero otherwise).\n",
    "\n",
    "Here we generated a uniform draw $ U $ on $ [0,1] $ and then used the fact that\n",
    "\n",
    "$$\n",
    "\\mathbb P\\{0 \\leq U < p\\} = p - 0 = p\n",
    "$$\n",
    "\n",
    "This means that $ X = \\mathbf 1\\{U < p\\} $ has the right distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b5d05ca",
   "metadata": {},
   "source": [
    "## Exercise 20.3\n",
    "\n",
    "We mentioned above that LLN can still hold sometimes when IID is violated.\n",
    "\n",
    "Let’s investigate this claim further.\n",
    "\n",
    "Consider the AR(1) process\n",
    "\n",
    "$$\n",
    "X_{t+1} = \\alpha + \\beta X_t + \\sigma \\epsilon _{t+1}\n",
    "$$\n",
    "\n",
    "where $ \\alpha, \\beta, \\sigma $ are constants and $ \\epsilon_1, \\epsilon_2,\n",
    "\\ldots $ are IID and standard normal.\n",
    "\n",
    "Suppose that\n",
    "\n",
    "$$\n",
    "X_0 \\sim N \\left(\\frac{\\alpha}{1-\\beta}, \\frac{\\sigma^2}{1-\\beta^2}\\right)\n",
    "$$\n",
    "\n",
    "This process violates the independence assumption of the LLN\n",
    "(since $ X_{t+1} $ depends on the value of $ X_t $).\n",
    "\n",
    "However, the next exercise teaches us that LLN type convergence of the sample\n",
    "mean to the population mean still occurs.\n",
    "\n",
    "1. Prove that the sequence $ X_1, X_2, \\ldots $ is identically distributed.  \n",
    "1. Show that LLN convergence holds using simulations with $ \\alpha = 0.8 $, $ \\beta = 0.2 $.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acac7656",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 20.3](https://intro.quantecon.org/#lln_ex3)\n",
    "\n",
    "**Q1 Solution**\n",
    "\n",
    "Regarding part 1, we claim that $ X_t $ has the same distribution as $ X_0 $ for\n",
    "all $ t $.\n",
    "\n",
    "To construct a proof, we suppose that the claim is true for $ X_t $.\n",
    "\n",
    "Now we claim it is also true for $ X_{t+1} $.\n",
    "\n",
    "Observe that we have the correct mean:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\mathbb E X_{t+1} &= \\alpha + \\beta \\mathbb E X_t \\\\\n",
    "    &= \\alpha + \\beta \\frac{\\alpha}{1-\\beta} \\\\\n",
    "    &= \\frac{\\alpha}{1-\\beta}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We also have the correct variance:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\mathrm{Var}(X_{t+1}) &= \\beta^2 \\mathrm{Var}(X_{t}) + \\sigma^2\\\\\n",
    "    &= \\frac{\\beta^2\\sigma^2}{1-\\beta^2} + \\sigma^2 \\\\\n",
    "    &= \\frac{\\sigma^2}{1-\\beta^2}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Finally, since both $ X_t $ and $ \\epsilon_0 $ are normally distributed and\n",
    "independent from each other, any linear combination of these two variables is\n",
    "also normally distributed.\n",
    "\n",
    "We have now shown that\n",
    "\n",
    "$$\n",
    "X_{t+1} \\sim \n",
    "    N \\left(\\frac{\\alpha}{1-\\beta}, \\frac{\\sigma^2}{1-\\beta^2}\\right)\n",
    "$$\n",
    "\n",
    "We can conclude this AR(1) process violates the independence assumption but is\n",
    "identically distributed.\n",
    "\n",
    "**Q2 Solution**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79e0ea01",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "σ = 10\n",
    "α = 0.8\n",
    "β = 0.2\n",
    "n = 100_000\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "x = np.ones(n)\n",
    "x[0] = st.norm.rvs(α/(1-β), α**2/(1-β**2))\n",
    "ϵ = st.norm.rvs(size=n+1)\n",
    "means = np.ones(n)\n",
    "means[0] = x[0]\n",
    "for t in range(n-1):\n",
    "    x[t+1] = α + β * x[t] + σ * ϵ[t+1]\n",
    "    means[t+1] = np.mean(x[:t+1])\n",
    "\n",
    "\n",
    "ax.scatter(range(100, n), means[100:n], s=10, alpha=0.5)\n",
    "\n",
    "ax.set_xlabel(r\"$n$\", size=12)\n",
    "ax.set_ylabel(r\"$\\bar X_n$\", size=12)\n",
    "yabs_max = max(ax.get_ylim(), key=abs)\n",
    "ax.axhline(y=α/(1-β), ls=\"--\", lw=3, \n",
    "           label=r\"$\\mu = \\frac{\\alpha}{1-\\beta}$\", \n",
    "           color = 'black')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e082906e",
   "metadata": {},
   "source": [
    "We see the convergence of $ \\bar x $ around $ \\mu $ even when the independence assumption is violated."
   ]
  }
 ],
 "metadata": {
  "date": 1750037600.481144,
  "filename": "lln_clt.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "LLN and CLT"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}