diff --git a/lectures/_static/quant-econ.bib b/lectures/_static/quant-econ.bib index 08a19b216..7e02441fe 100644 --- a/lectures/_static/quant-econ.bib +++ b/lectures/_static/quant-econ.bib @@ -2,6 +2,36 @@ QuantEcon Bibliography File used in conjuction with sphinxcontrib-bibtex package Note: Extended Information (like abstracts, doi, url's etc.) can be found in quant-econ-extendedinfo.bib file in _static/ ### + +@article{shannon1948mathematical, + title={A mathematical theory of communication}, + author={Shannon, Claude E}, + journal={The Bell system technical journal}, + volume={27}, + number={3}, + pages={379--423}, + year={1948}, + publisher={Nokia Bell Labs} +} + +@article{kullback1951information, + title={On Information and Sufficiency}, + author={Kullback, Solomon and Leibler, Richard A}, + journal={The Annals of Mathematical Statistics}, + volume={22}, + number={1}, + pages={79--86}, + year={1951}, + publisher={JSTOR} +} + +@book{kullback1997information, + title={Information theory and statistics}, + author={Kullback, Solomon}, + year={1997}, + publisher={Courier Corporation} +} + @book{friedman1953essays, title={Essays in positive economics}, author={Friedman, Milton}, diff --git a/lectures/_toc.yml b/lectures/_toc.yml index 7fb62a7e8..0103e311b 100644 --- a/lectures/_toc.yml +++ b/lectures/_toc.yml @@ -33,7 +33,10 @@ parts: - caption: Statistics and Information numbered: true chapters: + - file: divergence_measures - file: likelihood_ratio_process + - file: likelihood_ratio_process_2 + - file: likelihood_var - file: imp_sample - file: wald_friedman - file: wald_friedman_2 diff --git a/lectures/divergence_measures.md b/lectures/divergence_measures.md new file mode 100644 index 000000000..b95d30d4e --- /dev/null +++ b/lectures/divergence_measures.md @@ -0,0 +1,569 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.17.1 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +(divergence_measures)= +```{raw} jupyter +
+ + QuantEcon + +
+``` + +# Statistical Divergence Measures + +```{contents} Contents +:depth: 2 +``` + +## Overview + +A statistical divergence quantifies discrepancies between two distinct + probability distributions that can be challenging to distinguish for the following reason: + + * every event that has positive probability under one of the distributions also has positive probability under the other distribution + + * this means that there is no "smoking gun" event whose occurrence tells a statistician that one of the probability distributions surely governs the data + +A statistical divergence is a **function** that maps two probability distributions into a nonnegative real number. + +Statistical divergence functions play important roles in statistics, information theory, and what many people now call "machine learning". + +This lecture describes three divergence measures: + +* **Kullback–Leibler (KL) divergence** +* **Jensen–Shannon (JS) divergence** +* **Chernoff entropy** + +These will appear in several quantecon lectures. + +Let's start by importing the necessary Python tools. + +```{code-cell} ipython3 +import matplotlib.pyplot as plt +import numpy as np +from numba import vectorize, jit +from math import gamma +from scipy.integrate import quad +from scipy.optimize import minimize_scalar +import pandas as pd +from IPython.display import display, Math +``` + + + + + +## Primer on entropy, cross-entropy, KL divergence + +Before diving in, we'll introduce some useful concepts in a simple setting. + +We'll temporarily assume that $f$ and $g$ are two probability mass functions for discrete random variables +on state space $I = \{1, 2, \ldots, n\}$ that satisfy $f_i \geq 0, \sum_{i} f_i =1, g_i \geq 0, \sum_{i} g_i =1$. + +We follow some statisticians and information theorists who define the **surprise** or **surprisal** +associated with having observed a single draw $x = i$ from distribution $f$ as + +$$ +\log\left(\frac{1}{f_i}\right) +$$ + +They then define the **information** that you can anticipate to gather from observing a single realization +as the expected surprisal + +$$ +H(f) = \sum_i f_i \log\left(\frac{1}{f_i}\right). +$$ + +Claude Shannon {cite}`shannon1948mathematical` called $H(f)$ the **entropy** of distribution $f$. + + +```{note} +By maximizing $H(f)$ with respect to $\{f_1, f_2, \ldots, f_n\}$ subject to $\sum_i f_i = 1$, we can verify that the distribution +that maximizes entropy is the uniform distribution +$ +f_i = \frac{1}{n} . +$ +Entropy $H(f)$ for the uniform distribution evidently equals $- \log(n)$. +``` + + + +Kullback and Leibler {cite}`kullback1951information` define the amount of information that a single draw of $x$ provides for distinguishing $f$ from $g$ as the log likelihood ratio + +$$ +\log \frac{f(x)}{g(x)} +$$ + + + + + +The following two concepts are widely used to compare two distributions $f$ and $g$. + + + +**Cross-Entropy:** + +\begin{equation} +H(f,g) = -\sum_{i} f_i \log g_i +\end{equation} + + + +**Kullback-Leibler (KL) Divergence:** +\begin{equation} +D_{KL}(f \parallel g) = \sum_{i} f_i \log\left[\frac{f_i}{g_i}\right] +\end{equation} + +These concepts are related by the following equality. + +$$ +D_{KL}(f \parallel g) = H(f,g) - H(f) +$$ (eq:KLcross) + +To prove {eq}`eq:KLcross`, note that + + +\begin{align} +D_{KL}(f \parallel g) &= \sum_{i} f_i \log\left[\frac{f_i}{g_i}\right] \\ +&= \sum_{i} f_i \left[\log f_i - \log g_i\right] \\ +&= \sum_{i} f_i \log f_i - \sum_{i} f_i \log g_i \\ +&= -H(f) + H(f,g) \\ +&= H(f,g) - H(f) +\end{align} + +Remember that $H(f)$ is the anticipated surprisal from drawing $x$ from $f$. + +Then the above equation tells us that the KL divergence is an anticipated "excess surprise" that comes from anticipating that $x$ is drawn from $f$ when it is +actually drawn from $g$. + + +## Two Beta distributions: running example + +We'll use Beta distributions extensively to illustrate concepts. + +The Beta distribution is particularly convenient as it's defined on $[0,1]$ and exhibits diverse shapes by appropriately choosing its two parameters. + +The density of a Beta distribution with parameters $a$ and $b$ is given by + +$$ +f(z; a, b) = \frac{\Gamma(a+b) z^{a-1} (1-z)^{b-1}}{\Gamma(a) \Gamma(b)} +\quad \text{where} \quad +\Gamma(p) := \int_{0}^{\infty} x^{p-1} e^{-x} dx +$$ + +Let's define parameters and density functions in Python + +```{code-cell} ipython3 +# Parameters in the two Beta distributions +F_a, F_b = 1, 1 +G_a, G_b = 3, 1.2 + +@vectorize +def p(x, a, b): + r = gamma(a + b) / (gamma(a) * gamma(b)) + return r * x** (a-1) * (1 - x) ** (b-1) + +# The two density functions +f = jit(lambda x: p(x, F_a, F_b)) +g = jit(lambda x: p(x, G_a, G_b)) + +# Plot the distributions +x_range = np.linspace(0.001, 0.999, 1000) +f_vals = [f(x) for x in x_range] +g_vals = [g(x) for x in x_range] + +plt.figure(figsize=(10, 6)) +plt.plot(x_range, f_vals, 'b-', linewidth=2, label=r'$f(x) \sim \text{Beta}(1,1)$') +plt.plot(x_range, g_vals, 'r-', linewidth=2, label=r'$g(x) \sim \text{Beta}(3,1.2)$') + +# Fill overlap region +overlap = np.minimum(f_vals, g_vals) +plt.fill_between(x_range, 0, overlap, alpha=0.3, color='purple', label='overlap') + +plt.xlabel('x') +plt.ylabel('density') +plt.legend() +plt.show() +``` + + + +(rel_entropy)= +## Kullback–Leibler divergence + +Our first divergence function is the **Kullback–Leibler (KL) divergence**. + +For probability densities (or pmfs) $f$ and $g$ it is defined by + +$$ +D_{KL}(f\|g) = KL(f, g) = \int f(x) \log \frac{f(x)}{g(x)} \, dx. +$$ + +We can interpret $D_{KL}(f\|g)$ as the expected excess log loss (expected excess surprisal) incurred when we use $g$ while the data are generated by $f$. + +It has several important properties: + +- Non-negativity (Gibbs' inequality): $D_{KL}(f\|g) \ge 0$ with equality if and only if $f=g$ almost everywhere. +- Asymmetry: $D_{KL}(f\|g) \neq D_{KL}(g\|f)$ in general (hence it is not a metric) +- Information decomposition: + $D_{KL}(f\|g) = H(f,g) - H(f)$, where $H(f,g)$ is the cross entropy and $H(f)$ is the Shannon entropy of $f$. +- Chain rule: For joint distributions $f(x, y)$ and $g(x, y)$, + $D_{KL}(f(x,y)\|g(x,y)) = D_{KL}(f(x)\|g(x)) + E_{f}\left[D_{KL}(f(y|x)\|g(y|x))\right]$ + +KL divergence plays a central role in statistical inference, including model selection and hypothesis testing. + +{doc}`likelihood_ratio_process` describes a link between KL divergence and the expected log likelihood ratio, +and the lecture {doc}`wald_friedman` connects it to the test performance of the sequential probability ratio test. + +Let's compute the KL divergence between our example distributions $f$ and $g$. + +```{code-cell} ipython3 +def compute_KL(f, g): + """ + Compute KL divergence KL(f, g) via numerical integration + """ + def integrand(w): + fw = f(w) + gw = g(w) + return fw * np.log(fw / gw) + val, _ = quad(integrand, 1e-5, 1-1e-5) + return val + +# Compute KL divergences between our example distributions +kl_fg = compute_KL(f, g) +kl_gf = compute_KL(g, f) + +print(f"KL(f, g) = {kl_fg:.4f}") +print(f"KL(g, f) = {kl_gf:.4f}") +``` + +The asymmetry of KL divergence has important practical implications. + +$D_{KL}(f\|g)$ penalizes regions where $f > 0$ but $g$ is close to zero, reflecting the cost of using $g$ to model $f$ and vice versa. + +## Jensen-Shannon divergence + +Sometimes we want a symmetric measure of divergence that captures the difference between two distributions without favoring one over the other. + +This often arises in applications like clustering, where we want to compare distributions without assuming one is the true model. + +The **Jensen-Shannon (JS) divergence** symmetrizes KL divergence by comparing both distributions to their mixture: + +$$ +JS(f,g) = \frac{1}{2} D_{KL}(f\|m) + \frac{1}{2} D_{KL}(g\|m), \quad m = \frac{1}{2}(f+g). +$$ + +where $m$ is a mixture distribution that averages $f$ and $g$ + +Let's also visualize the mixture distribution $m$: + +```{code-cell} ipython3 +def m(x): + return 0.5 * (f(x) + g(x)) + +m_vals = [m(x) for x in x_range] + +plt.figure(figsize=(10, 6)) +plt.plot(x_range, f_vals, 'b-', linewidth=2, label=r'$f(x)$') +plt.plot(x_range, g_vals, 'r-', linewidth=2, label=r'$g(x)$') +plt.plot(x_range, m_vals, 'g--', linewidth=2, label=r'$m(x) = \frac{1}{2}(f(x) + g(x))$') + +plt.xlabel('x') +plt.ylabel('density') +plt.legend() +plt.show() +``` + +The JS divergence has several useful properties: + +- Symmetry: $JS(f,g)=JS(g,f)$. +- Boundedness: $0 \le JS(f,g) \le \log 2$. +- Its square root $\sqrt{JS}$ is a metric (Jensen–Shannon distance) on the space of probability distributions. +- JS divergence equals the mutual information between a binary random variable $Z \sim \text{Bernoulli}(1/2)$ indicating the source and a sample $X$ drawn from $f$ if $Z=0$ or from $g$ if $Z=1$. + +The Jensen–Shannon divergence plays a key role in the optimization of certain +generative models, as it is bounded, symmetric, and smoother than KL divergence, +often providing more stable gradients for training. + +Let's compute the JS divergence between our example distributions $f$ and $g$ + +```{code-cell} ipython3 +def compute_JS(f, g): + """Compute Jensen-Shannon divergence.""" + def m(w): + return 0.5 * (f(w) + g(w)) + js_div = 0.5 * compute_KL(f, m) + 0.5 * compute_KL(g, m) + return js_div + +js_div = compute_JS(f, g) +print(f"Jensen-Shannon divergence JS(f,g) = {js_div:.4f}") +``` + +We can easily generalize to more than two distributions using the generalized Jensen-Shannon divergence with weights $\alpha = (\alpha_i)_{i=1}^{n}$: + +$$ +JS_\alpha(f_1, \ldots, f_n) = +H\left(\sum_{i=1}^n \alpha_i f_i\right) - \sum_{i=1}^n \alpha_i H(f_i) +$$ + +where: +- $\alpha_i \geq 0$ and $\sum_{i=1}^n \alpha_i = 1$, and +- $H(f) = -\int f(x) \log f(x) dx$ is the **Shannon entropy** of distribution $f$ + +## Chernoff entropy + +Chernoff entropy originates from early applications of the [theory of large deviations](https://en.wikipedia.org/wiki/Large_deviations_theory), which refines central limit approximations by providing exponential decay rates for rare events. + + +For densities $f$ and $g$ the Chernoff entropy is + +$$ +C(f,g) = - \log \min_{\phi \in (0,1)} \int f^{\phi}(x) g^{1-\phi}(x) \, dx. +$$ + +Remarks: + +- The inner integral is the **Chernoff coefficient**. +- At $\phi=1/2$ it becomes the **Bhattacharyya coefficient** $\int \sqrt{f g}$. +- In binary hypothesis testing with $T$ iid observations, the optimal error probability decays as $e^{-C(f,g) T}$. + +We will see an example of the third point in the lecture {doc}`likelihood_ratio_process`, +where we study the Chernoff entropy in the context of model selection. + +Let's compute the Chernoff entropy between our example distributions $f$ and $g$. + +```{code-cell} ipython3 +def chernoff_integrand(ϕ, f, g): + """Integral entering Chernoff entropy for a given ϕ.""" + def integrand(w): + return f(w)**ϕ * g(w)**(1-ϕ) + result, _ = quad(integrand, 1e-5, 1-1e-5) + return result + +def compute_chernoff_entropy(f, g): + """Compute Chernoff entropy C(f,g).""" + def objective(ϕ): + return chernoff_integrand(ϕ, f, g) + result = minimize_scalar(objective, bounds=(1e-5, 1-1e-5), method='bounded') + min_value = result.fun + ϕ_optimal = result.x + chernoff_entropy = -np.log(min_value) + return chernoff_entropy, ϕ_optimal + +C_fg, ϕ_optimal = compute_chernoff_entropy(f, g) +print(f"Chernoff entropy C(f,g) = {C_fg:.4f}") +print(f"Optimal ϕ = {ϕ_optimal:.4f}") +``` + +## Comparing divergence measures + +We now compare these measures across several pairs of Beta distributions + +```{code-cell} ipython3 +:tags: [hide-input] + +distribution_pairs = [ + # (f_params, g_params) + ((1, 1), (0.1, 0.2)), + ((1, 1), (0.3, 0.3)), + ((1, 1), (0.3, 0.4)), + ((1, 1), (0.5, 0.5)), + ((1, 1), (0.7, 0.6)), + ((1, 1), (0.9, 0.8)), + ((1, 1), (1.1, 1.05)), + ((1, 1), (1.2, 1.1)), + ((1, 1), (1.5, 1.2)), + ((1, 1), (2, 1.5)), + ((1, 1), (2.5, 1.8)), + ((1, 1), (3, 1.2)), + ((1, 1), (4, 1)), + ((1, 1), (5, 1)) +] + +# Create comparison table +results = [] +for i, ((f_a, f_b), (g_a, g_b)) in enumerate(distribution_pairs): + f = jit(lambda x, a=f_a, b=f_b: p(x, a, b)) + g = jit(lambda x, a=g_a, b=g_b: p(x, a, b)) + kl_fg = compute_KL(f, g) + kl_gf = compute_KL(g, f) + js_div = compute_JS(f, g) + chernoff_ent, _ = compute_chernoff_entropy(f, g) + results.append({ + 'Pair (f, g)': f"\\text{{Beta}}({f_a},{f_b}), \\text{{Beta}}({g_a},{g_b})", + 'KL(f, g)': f"{kl_fg:.4f}", + 'KL(g, f)': f"{kl_gf:.4f}", + 'JS': f"{js_div:.4f}", + 'C': f"{chernoff_ent:.4f}" + }) + +df = pd.DataFrame(results) +# Sort by JS divergence +df['JS_numeric'] = df['JS'].astype(float) +df = df.sort_values('JS_numeric').drop('JS_numeric', axis=1) + +columns = ' & '.join([f'\\text{{{col}}}' for col in df.columns]) +rows = ' \\\\\n'.join( + [' & '.join([f'{val}' for val in row]) + for row in df.values]) + +latex_code = rf""" +\begin{{array}}{{lcccc}} +{columns} \\ +\hline +{rows} +\end{{array}} +""" + +display(Math(latex_code)) +``` + +We can clearly see co-movement across the divergence measures as we vary the parameters of the Beta distributions. + +Next we visualize relationships among KL, JS, and Chernoff entropy. + +```{code-cell} ipython3 +kl_fg_values = [float(result['KL(f, g)']) for result in results] +js_values = [float(result['JS']) for result in results] +chernoff_values = [float(result['C']) for result in results] + +fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + +axes[0].scatter(kl_fg_values, js_values, alpha=0.7, s=60) +axes[0].set_xlabel('KL divergence KL(f, g)') +axes[0].set_ylabel('JS divergence') +axes[0].set_title('JS divergence vs KL divergence') + +axes[1].scatter(js_values, chernoff_values, alpha=0.7, s=60) +axes[1].set_xlabel('JS divergence') +axes[1].set_ylabel('Chernoff entropy') +axes[1].set_title('Chernoff entropy vs JS divergence') + +plt.tight_layout() +plt.show() +``` + +We now generate plots illustrating how overlap visually diminishes as divergence measures increase. + + +```{code-cell} ipython3 +param_grid = [ + ((1, 1), (1, 1)), + ((1, 1), (1.5, 1.2)), + ((1, 1), (2, 1.5)), + ((1, 1), (3, 1.2)), + ((1, 1), (5, 1)), + ((1, 1), (0.3, 0.3)) +] +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +def plot_dist_diff(para_grid): + """Plot overlap of selected Beta distribution pairs.""" + + fig, axes = plt.subplots(3, 2, figsize=(15, 12)) + divergence_data = [] + for i, ((f_a, f_b), (g_a, g_b)) in enumerate(param_grid): + row, col = divmod(i, 2) + f = jit(lambda x, a=f_a, b=f_b: p(x, a, b)) + g = jit(lambda x, a=g_a, b=g_b: p(x, a, b)) + kl_fg = compute_KL(f, g) + js_div = compute_JS(f, g) + chernoff_ent, _ = compute_chernoff_entropy(f, g) + divergence_data.append({ + 'f_params': (f_a, f_b), + 'g_params': (g_a, g_b), + 'kl_fg': kl_fg, + 'js_div': js_div, + 'chernoff': chernoff_ent + }) + x_range = np.linspace(0, 1, 200) + f_vals = [f(x) for x in x_range] + g_vals = [g(x) for x in x_range] + axes[row, col].plot(x_range, f_vals, 'b-', + linewidth=2, label=f'f ~ Beta({f_a},{f_b})') + axes[row, col].plot(x_range, g_vals, 'r-', + linewidth=2, label=f'g ~ Beta({g_a},{g_b})') + overlap = np.minimum(f_vals, g_vals) + axes[row, col].fill_between(x_range, 0, + overlap, alpha=0.3, color='purple', label='overlap') + axes[row, col].set_title( + f'KL(f,g)={kl_fg:.3f}, JS={js_div:.3f}, C={chernoff_ent:.3f}', + fontsize=12) + axes[row, col].legend(fontsize=12) + plt.tight_layout() + plt.show() + return divergence_data + +divergence_data = plot_dist_diff(param_grid) +``` + + + +## KL divergence and maximum-likelihood estimation + + +Given a sample of $n$ observations $X = \{x_1, x_2, \ldots, x_n\}$, the **empirical distribution** is + +$$p_e(x) = \frac{1}{n} \sum_{i=1}^n \delta(x - x_i)$$ + +where $\delta(x - x_i)$ is the Dirac delta function centered at $x_i$: + +$$ +\delta(x - x_i) = \begin{cases} ++\infty & \text{if } x = x_i \\ +0 & \text{if } x \neq x_i +\end{cases} +$$ + +- **Discrete probability measure**: Assigns probability $\frac{1}{n}$ to each observed data point +- **Empirical expectation**: $\langle X \rangle_{p_e} = \frac{1}{n} \sum_{i=1}^n x_i = \bar{\mu}$ +- **Support**: Only on the observed data points $\{x_1, x_2, \ldots, x_n\}$ + + +The KL divergence from the empirical distribution $p_e$ to a parametric model $p_\theta(x)$ is: + +$$D_{KL}(p_e \parallel p_\theta) = \int p_e(x) \log \frac{p_e(x)}{p_\theta(x)} dx$$ + +Using the mathematics of the Dirac delta function, it follows that + +$$D_{KL}(p_e \parallel p_\theta) = \sum_{i=1}^n \frac{1}{n} \log \frac{\left(\frac{1}{n}\right)}{p_\theta(x_i)}$$ + +$$= \frac{1}{n} \sum_{i=1}^n \log \frac{1}{n} - \frac{1}{n} \sum_{i=1}^n \log p_\theta(x_i)$$ + +$$= -\log n - \frac{1}{n} \sum_{i=1}^n \log p_\theta(x_i)$$ + +Since the log-likelihood function for parameter $\theta$ is: + +$$ +\ell(\theta; X) = \sum_{i=1}^n \log p_\theta(x_i) , +$$ + +it follows that maximum likelihood chooses parameters to minimize + +$$ D_{KL}(p_e \parallel p_\theta) $$ + + +Thus, MLE is equivalent to minimizing the KL divergence from the empirical distribution to the statistical model $p_\theta$. + +## Related lectures + +This lecture has introduced tools that we'll encounter elsewhere. + +- Other quantecon lectures that apply connections between divergence measures and statistical inference include {doc}`likelihood_ratio_process`, {doc}`wald_friedman`, and {doc}`mix_model`. + +- Statistical divergence functions also take center stage in {doc}`likelihood_ratio_process_2` that studies Lawrence Blume and David Easley's model of heterogeneous beliefs and financial markets. diff --git a/lectures/likelihood_bayes.md b/lectures/likelihood_bayes.md index 2787dbe1b..0d2a4fe1b 100644 --- a/lectures/likelihood_bayes.md +++ b/lectures/likelihood_bayes.md @@ -4,7 +4,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.2 + jupytext_version: 1.17.1 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -194,8 +194,6 @@ l_arr_f = simulate(F_a, F_b, N=50000) l_seq_f = np.cumprod(l_arr_f, axis=1) ``` - - ## Likelihood Ratio Processes and Bayes’ Law Let $\pi_0 \in [0,1]$ be a Bayesian statistician's prior probability that nature generates $w^t$ as a sequence of i.i.d. draws from @@ -504,7 +502,7 @@ We will first use this sequence to study how $\pi_t$ behaves. ```{note} Later, we can use it to study how a statistician who knows that nature generates data from an $x$-mixture of $f$ and $g$ could construct maximum likelihood or Bayesian estimators of $x$ along with the free parameters of $f$ and $g$. -``` +``` ```{code-cell} ipython3 x_true = 0.5 @@ -1027,4 +1025,4 @@ The conditional variance is nearly zero only when the agent is almost sure that ## Related Lectures This lecture has been devoted to building some useful infrastructure that will help us understand inferences that are the foundations of -results described in {doc}`this lecture ` and {doc}`this lecture ` and {doc}`this lecture `. \ No newline at end of file +results described in {doc}`this lecture ` and {doc}`this lecture ` and {doc}`this lecture `. diff --git a/lectures/likelihood_ratio_process.md b/lectures/likelihood_ratio_process.md index c9a7357bb..5a7213702 100644 --- a/lectures/likelihood_ratio_process.md +++ b/lectures/likelihood_ratio_process.md @@ -38,7 +38,6 @@ Among the things that we'll learn are * How a likelihood ratio process is a key ingredient in frequentist hypothesis testing * How a **receiver operator characteristic curve** summarizes information about a false alarm probability and power in frequentist hypothesis testing * How a statistician can combine frequentist probabilities of type I and type II errors to form posterior probabilities of mistakes in a model selection or in an individual-classification problem -* How likelihood ratios helped Lawrence Blume and David Easley formulate an answer to ''If you're so smart, why aren't you rich?'' {cite}`blume2006if` * How to use a Kullback-Leibler divergence to quantify the difference between two probability distributions with the same support * How during World War II the United States Navy devised a decision rule for doing quality control on lots of ammunition, a topic that sets the stage for {doc}`this lecture ` * A peculiar property of likelihood ratio processes @@ -127,36 +126,71 @@ ratio process by generating a sequence $w^t$ from one of the two probability distributions, for example, a sequence of IID draws from $g$. ```{code-cell} ipython3 -# Parameters in the two Beta distributions. +# Parameters for the two Beta distributions F_a, F_b = 1, 1 G_a, G_b = 3, 1.2 @vectorize def p(x, a, b): + """Beta distribution density function.""" r = gamma(a + b) / (gamma(a) * gamma(b)) return r * x** (a-1) * (1 - x) ** (b-1) -# The two density functions. f = jit(lambda x: p(x, F_a, F_b)) g = jit(lambda x: p(x, G_a, G_b)) -``` -```{code-cell} ipython3 -@jit -def simulate(a, b, T=50, N=500): - ''' - Generate N sets of T observations of the likelihood ratio, - return as N x T matrix. - ''' +def create_beta_density(a, b): + """Create a beta density function with specified parameters.""" + return jit(lambda x: p(x, a, b)) - l_arr = np.empty((N, T)) +def likelihood_ratio(w, f_func, g_func): + """Compute likelihood ratio for observation(s) w.""" + return f_func(w) / g_func(w) +@jit +def simulate_likelihood_ratios(a, b, f_func, g_func, T=50, N=500): + """ + Generate N sets of T observations of the likelihood ratio. + """ + l_arr = np.empty((N, T)) for i in range(N): for j in range(T): w = np.random.beta(a, b) - l_arr[i, j] = f(w) / g(w) - + l_arr[i, j] = f_func(w) / g_func(w) return l_arr + +def simulate_sequences(distribution, f_func, g_func, + F_params=(1, 1), G_params=(3, 1.2), T=50, N=500): + """ + Generate N sequences of T observations from specified distribution. + """ + if distribution == 'f': + a, b = F_params + elif distribution == 'g': + a, b = G_params + else: + raise ValueError("distribution must be 'f' or 'g'") + + l_arr = simulate_likelihood_ratios(a, b, f_func, g_func, T, N) + l_seq = np.cumprod(l_arr, axis=1) + return l_arr, l_seq + +def plot_likelihood_paths(l_seq, title="Likelihood ratio paths", + ylim=None, n_paths=None): + """Plot likelihood ratio paths.""" + N, T = l_seq.shape + n_show = n_paths or min(N, 100) + + plt.figure(figsize=(10, 6)) + for i in range(n_show): + plt.plot(range(T), l_seq[i, :], color='b', lw=0.8, alpha=0.5) + + if ylim: + plt.ylim(ylim) + plt.title(title) + plt.xlabel('t') + plt.ylabel('$L(w^t)$') + plt.show() ``` (nature_likeli)= @@ -166,19 +200,11 @@ We first simulate the likelihood ratio process when nature permanently draws from $g$. ```{code-cell} ipython3 -l_arr_g = simulate(G_a, G_b) -l_seq_g = np.cumprod(l_arr_g, axis=1) -``` - -```{code-cell} ipython3 -N, T = l_arr_g.shape - -for i in range(N): - - plt.plot(range(T), l_seq_g[i, :], color='b', lw=0.8, alpha=0.5) - -plt.ylim([0, 3]) -plt.title("$L(w^{t})$ paths"); +# Simulate when nature draws from g +l_arr_g, l_seq_g = simulate_sequences('g', f, g, (F_a, F_b), (G_a, G_b)) +plot_likelihood_paths(l_seq_g, + title="$L(w^{t})$ paths when nature draws from g", + ylim=[0, 3]) ``` Evidently, as sample length $T$ grows, most probability mass @@ -189,6 +215,7 @@ paths $L\left(w^{t}\right)$ that fall in the interval $\left[0, 0.01\right]$. ```{code-cell} ipython3 +N, T = l_arr_g.shape plt.plot(range(T), np.sum(l_seq_g <= 0.01, axis=0) / N) plt.show() ``` @@ -254,8 +281,8 @@ calculate the unconditional mean of $L\left(w^t\right)$ by averaging across these many paths at each $t$. ```{code-cell} ipython3 -l_arr_g = simulate(G_a, G_b, N=50000) -l_seq_g = np.cumprod(l_arr_g, axis=1) +l_arr_g, l_seq_g = simulate_sequences('g', + f, g, (F_a, F_b), (G_a, G_b), N=50000) ``` It would be useful to use simulations to verify that unconditional means @@ -301,8 +328,9 @@ Simulations below confirm this conclusion. Please note the scale of the $y$ axis. ```{code-cell} ipython3 -l_arr_f = simulate(F_a, F_b, N=50000) -l_seq_f = np.cumprod(l_arr_f, axis=1) +# Simulate when nature draws from f +l_arr_f, l_seq_f = simulate_sequences('f', f, g, + (F_a, F_b), (G_a, G_b), N=50000) ``` ```{code-cell} ipython3 @@ -447,29 +475,35 @@ is what makes it possible eventually to distinguish $q=f$ from $q=g$. ```{code-cell} ipython3 -fig, axs = plt.subplots(2, 2, figsize=(12, 8)) -fig.suptitle('distribution of $log(L(w^t))$ under f or under g', fontsize=15) - -for i, t in enumerate([1, 7, 14, 21]): - nr = i // 2 - nc = i % 2 - - axs[nr, nc].axvline(np.log(c), color="k", ls="--") - - hist_f, x_f = np.histogram(np.log(l_seq_f[:, t]), 200, density=True) - hist_g, x_g = np.histogram(np.log(l_seq_g[:, t]), 200, density=True) - - axs[nr, nc].plot(x_f[1:], hist_f, label="dist under f") - axs[nr, nc].plot(x_g[1:], hist_g, label="dist under g") - - for i, (x, hist, label) in enumerate(zip([x_f, x_g], [hist_f, hist_g], ["Type I error", "Type II error"])): - ind = x[1:] <= np.log(c) if i == 0 else x[1:] > np.log(c) - axs[nr, nc].fill_between(x[1:][ind], hist[ind], alpha=0.5, label=label) - - axs[nr, nc].legend() - axs[nr, nc].set_title(f"t={t}") +def plot_log_histograms(l_seq_f, l_seq_g, c=1, time_points=[1, 7, 14, 21]): + """Plot log likelihood ratio histograms.""" + fig, axs = plt.subplots(2, 2, figsize=(12, 8)) + + for i, t in enumerate(time_points): + nr, nc = i // 2, i % 2 + + axs[nr, nc].axvline(np.log(c), color="k", ls="--") + + hist_f, x_f = np.histogram(np.log(l_seq_f[:, t]), 200, density=True) + hist_g, x_g = np.histogram(np.log(l_seq_g[:, t]), 200, density=True) + + axs[nr, nc].plot(x_f[1:], hist_f, label="dist under f") + axs[nr, nc].plot(x_g[1:], hist_g, label="dist under g") + + # Fill error regions + for j, (x, hist, label) in enumerate( + zip([x_f, x_g], [hist_f, hist_g], + ["Type I error", "Type II error"])): + ind = x[1:] <= np.log(c) if j == 0 else x[1:] > np.log(c) + axs[nr, nc].fill_between(x[1:][ind], hist[ind], + alpha=0.5, label=label) + + axs[nr, nc].legend() + axs[nr, nc].set_title(f"t={t}") + + plt.show() -plt.show() +plot_log_histograms(l_seq_f, l_seq_g, c=c) ``` In the above graphs, @@ -485,19 +519,43 @@ $t$ * the probability of a false alarm monotonically decreases with increases in $t$. ```{code-cell} ipython3 -PD = np.empty(T) -PFA = np.empty(T) -for t in range(T): - PD[t] = np.sum(l_seq_g[:, t] < c) / N - PFA[t] = np.sum(l_seq_f[:, t] < c) / N +def compute_error_probabilities(l_seq_f, l_seq_g, c=1): + """ + Compute Type I and Type II error probabilities. + """ + N, T = l_seq_f.shape + + # Type I error (false alarm) - reject H0 when true + PFA = np.array([np.sum(l_seq_f[:, t] < c) / N for t in range(T)]) + + # Type II error - accept H0 when false + beta = np.array([np.sum(l_seq_g[:, t] >= c) / N for t in range(T)]) + + # Probability of detection (power) + PD = np.array([np.sum(l_seq_g[:, t] < c) / N for t in range(T)]) + + return { + 'alpha': PFA, + 'beta': beta, + 'PD': PD, + 'PFA': PFA + } -plt.plot(range(T), PD, label="Probability of detection") -plt.plot(range(T), PFA, label="Probability of false alarm") -plt.xlabel("t") -plt.title("$c=1$") -plt.legend() -plt.show() +def plot_error_probabilities(error_dict, T, c=1, title_suffix=""): + """Plot error probabilities over time.""" + plt.figure(figsize=(10, 6)) + plt.plot(range(T), error_dict['PD'], label="Probability of detection") + plt.plot(range(T), error_dict['PFA'], label="Probability of false alarm") + plt.xlabel("t") + plt.ylabel("Probability") + plt.title(f"Error Probabilities (c={c}){title_suffix}") + plt.legend() + plt.show() + +error_probs = compute_error_probabilities(l_seq_f, l_seq_g, c=c) +N, T = l_seq_f.shape +plot_error_probabilities(error_probs, T, c) ``` For a given sample size $t$, the threshold $c$ uniquely pins down probabilities @@ -513,24 +571,32 @@ Below, we plot receiver operating characteristic curves for different sample sizes $t$. ```{code-cell} ipython3 -PFA = np.arange(0, 100, 1) - -for t in range(1, 15, 4): - percentile = np.percentile(l_seq_f[:, t], PFA) - PD = [np.sum(l_seq_g[:, t] < p) / N for p in percentile] - - plt.plot(PFA / 100, PD, label=f"t={t}") +def plot_roc_curves(l_seq_f, l_seq_g, t_values=[1, 5, 9, 13], N=None): + """Plot ROC curves for different sample sizes.""" + if N is None: + N = l_seq_f.shape[0] + + PFA = np.arange(0, 100, 1) + + plt.figure(figsize=(10, 6)) + for t in t_values: + percentile = np.percentile(l_seq_f[:, t], PFA) + PD = [np.sum(l_seq_g[:, t] < p) / N for p in percentile] + plt.plot(PFA / 100, PD, label=f"t={t}") + + plt.scatter(0, 1, label="perfect detection") + plt.plot([0, 1], [0, 1], color='k', ls='--', label="random detection") + + plt.arrow(0.5, 0.5, -0.15, 0.15, head_width=0.03) + plt.text(0.35, 0.7, "better") + plt.xlabel("Probability of false alarm") + plt.ylabel("Probability of detection") + plt.legend() + plt.title("ROC Curve") + plt.show() -plt.scatter(0, 1, label="perfect detection") -plt.plot([0, 1], [0, 1], color='k', ls='--', label="random detection") -plt.arrow(0.5, 0.5, -0.15, 0.15, head_width=0.03) -plt.text(0.35, 0.7, "better") -plt.xlabel("Probability of false alarm") -plt.ylabel("Probability of detection") -plt.legend() -plt.title("Receiver Operating Characteristic Curve") -plt.show() +plot_roc_curves(l_seq_f, l_seq_g, t_values=range(1, 15, 4), N=N) ``` Notice that as $t$ increases, we are assured a larger probability @@ -589,9 +655,8 @@ control tests during World War II. A Navy Captain who had been ordered to perform tests of this kind had doubts about it that he presented to Milton Friedman, as we describe in {doc}`this lecture `. - -(rel_entropy)= -## Kullback–Leibler Divergence +(llr_h)= +### A third distribution $h$ Now let's consider a case in which neither $g$ nor $f$ generates the data. @@ -601,11 +666,7 @@ Instead, a third distribution $h$ does. Let's study how accumulated likelihood ratios $L$ behave when $h$ governs the data. -A key tool here is called **Kullback–Leibler divergence**. - -It is also called **relative entropy**. - -It measures how one probability distribution differs from another. +A key tool here is called **Kullback–Leibler divergence** we studied in {doc}`divergence_measures`. In our application, we want to measure how much $f$ or $g$ diverges from $h$ @@ -623,16 +684,13 @@ $$ $$ \begin{aligned} -K_{g} = D_{KL}\bigl(h\|g\bigr) = KL(h,g) +K_{g} = D_{KL}\bigl(h\|g\bigr) = KL(h, g) &= E_{h}\left[\log\frac{h(w)}{g(w)}\right] \\ &= \int \log\left(\frac{h(w)}{g(w)}\right)h(w)dw . \end{aligned} $$ -+++ - -Let's compute the Kullback–Leibler discrepancies by quadrature -integration. +Let's compute the Kullback–Leibler discrepancies using the same code in {doc}`divergence_measures`. ```{code-cell} ipython3 def compute_KL(f, g): @@ -642,19 +700,13 @@ def compute_KL(f, g): integrand = lambda w: f(w) * np.log(f(w) / g(w)) val, _ = quad(integrand, 1e-5, 1-1e-5) return val -``` - -Next we create a helper function to compute KL divergence with respect to a reference distribution $h$ -```{code-cell} ipython3 def compute_KL_h(h, f, g): """ - Compute KL divergence with reference distribution h + Compute KL divergences with respect to reference distribution h """ - Kf = compute_KL(h, f) Kg = compute_KL(h, g) - return Kf, Kg ``` @@ -666,13 +718,11 @@ There is a mathematical relationship between likelihood ratios and KL divergence When data is generated by distribution $h$, the expected log likelihood ratio is: $$ -\frac{1}{t} E_{h}\!\bigl[\log L_t\bigr] = KL(h, g) - KL(h, f) = K_g - K_f +\frac{1}{t} E_{h}\!\bigl[\log L_t\bigr] = K_g - K_f $$ (eq:kl_likelihood_link) where $L_t=\prod_{j=1}^{t}\frac{f(w_j)}{g(w_j)}$ is the likelihood ratio process. -(For the proof, see [this note](https://nowak.ece.wisc.edu/ece830/ece830_fall11_lecture7.pdf).) - Equation {eq}`eq:kl_likelihood_link` tells us that: - When $K_g < K_f$ (i.e., $g$ is closer to $h$ than $f$ is), the expected log likelihood ratio is negative, so $L\left(w^t\right) \rightarrow 0$. - When $K_g > K_f$ (i.e., $f$ is closer to $h$ than $g$ is), the expected log likelihood ratio is positive, so $L\left(w^t\right) \rightarrow + \infty$. @@ -774,892 +824,337 @@ Note that - In the first figure, $\log L(w^t)$ diverges to $\infty$ because $K_g > K_f$. - In the second figure, we still have $K_g > K_f$, but the difference is smaller, so $L(w^t)$ diverges to infinity at a slower pace. - In the last figure, $\log L(w^t)$ diverges to $-\infty$ because $K_g < K_f$. -- The black dotted line, $t \left(KL(h,g) - KL(h, f)\right)$, closely fits the paths verifying {eq}`eq:kl_likelihood_link`. +- The black dotted line, $t \left(D_{KL}(h\|g) - D_{KL}(h\|f)\right)$, closely fits the paths verifying {eq}`eq:kl_likelihood_link`. These observations align with the theory. -In the [next section](hetero_agent), we will see an application of these ideas. - - -(hetero_agent)= -## Heterogeneous Beliefs and Financial Markets - -A likelihood ratio process lies behind Lawrence Blume and David Easley's answer to their question -''If you're so smart, why aren't you rich?'' {cite}`blume2006if`. - -Blume and Easley constructed formal models to study how differences of opinions about probabilities governing risky income processes would influence outcomes and be reflected in prices of stocks, bonds, and insurance policies that individuals use to share and hedge risks. - -```{note} -{cite}`alchian1950uncertainty` and {cite}`friedman1953essays` can conjectured that, by rewarding traders with more realistic probability models, competitive markets in financial securities put wealth in the hands of better informed traders and help -make prices of risky assets reflect realistic probability assessments. -``` - +In {doc}`likelihood_ratio_process_2`, we will see an application of these ideas. -Here we'll provide an example that illustrates basic components of Blume and Easley's analysis. +## Hypothesis testing and classification -We'll focus only on their analysis of an environment with complete markets in which trades in all conceivable risky securities are possible. - -We'll study two alternative arrangements: - -* perfect socialism in which individuals surrender their endowments of consumption goods each period to a central planner who then dictatorially allocates those goods -* a decentralized system of competitive markets in which selfish price-taking individuals voluntarily trade with each other in competitive markets - -The fundamental theorems of welfare economics will apply and assure us that these two arrangements end up producing exactly the same allocation of consumption goods to individuals **provided** that the social planner assigns an appropriate set of **Pareto weights**. - -```{note} -You can learn about how the two welfare theorems are applied in modern macroeconomic models in {doc}`this lecture on a planning problem ` and {doc}`this lecture on a related competitive equilibrium `. -``` +This section discusses another application of likelihood ratio processes. -### The setting +We describe how a statistician can combine frequentist probabilities of type I and type II errors in order to -Let the random variable $s_t \in (0,1)$ at time $t =0, 1, 2, \ldots$ be distributed according to the same Beta distribution with parameters -$\theta = \{\theta_1, \theta_2\}$. +* compute an anticipated frequency of selecting a wrong model based on a sample length $T$ +* compute an anticipated error rate in a classification problem -We'll denote this probability density as +We consider a situation in which nature generates data by mixing known densities $f$ and $g$ with known mixing +parameter $\pi_{-1} \in (0,1)$ so that the random variable $w$ is drawn from the density $$ -\pi(s_t|\theta) +h (w) = \pi_{-1} f(w) + (1-\pi_{-1}) g(w) $$ -Below, we'll often just write $\pi(s_t)$ instead of $\pi(s_t|\theta)$ to save space. +We assume that the statistician knows the densities $f$ and $g$ and also the mixing parameter $\pi_{-1}$. -Let $s_t \equiv y_t^1$ be the endowment of a nonstorable consumption good that a person we'll call "agent 1" receives at time $t$. +Below, we'll set $\pi_{-1} = .5$, although much of the analysis would follow through with other settings of $\pi_{-1} \in (0,1)$. -Let a history $s^t = [s_t, s_{t-1}, \ldots, s_0]$ be a sequence of i.i.d. random variables with joint distribution +We assume that $f$ and $g$ both put positive probabilities on the same intervals of possible realizations of the random variable $W$. -$$ -\pi_t(s^t) = \pi(s_t) \pi(s_{t-1}) \cdots \pi(s_0) -$$ + -So in our example, the history $s^t$ is a comprehensive record of agent $1$'s endowments of the consumption good from time $0$ up to time $t$. +In the simulations below, we specify that $f$ is a $\text{Beta}(1, 1)$ distribution and that $g$ is $\text{Beta}(3, 1.2)$ distribution. -If agent $1$ were to live on an island by himself, agent $1$'s consumption $c^1(s_t)$ at time $t$ is +We consider two alternative timing protocols. -$$c^1(s_t) = y_t^1 = s_t. $$ + * Timing protocol 1 is for the model selection problem + * Timing protocol 2 is for the individual classification problem -But in our model, agent 1 is not alone. +**Timing Protocol 1:** Nature flips a coin only **once** at time $t=-1$ and with probability $\pi_{-1}$ generates a sequence $\{w_t\}_{t=1}^T$ +of IID draws from $f$ and with probability $1-\pi_{-1}$ generates a sequence $\{w_t\}_{t=1}^T$ +of IID draws from $g$. -### Nature and agents' beliefs +**Timing Protocol 2.** Nature flips a coin **often**. At each time $t \geq 0$, nature flips a coin and with probability $\pi_{-1}$ draws $w_t$ from $f$ and with probability $1-\pi_{-1}$ draws $w_t$ from $g$. -Nature draws i.i.d. sequences $\{s_t\}_{t=0}^\infty$ from $\pi_t(s^t)$. +Here is Python code that we'll use to implement timing protocol 1 and 2 -* so $\pi$ without a superscript is nature's model -* but in addition to nature, there are other entities inside our model -- artificial people that we call "agents" -* each agent has a sequence of probability distributions over $s^t$ for $t=0, \ldots$ -* agent $i$ thinks that nature draws i.i.d. sequences $\{s_t\}_{t=0}^\infty$ from $\{\pi_t^i(s^t)\}_{t=0}^\infty$ - * agent $i$ is mistaken unless $\pi_t^i(s^t) = \pi_t(s^t)$ +```{code-cell} ipython3 +def protocol_1(π_minus_1, T, N=1000, F_params=(1, 1), G_params=(3, 1.2)): + """ + Simulate Protocol 1: Nature decides once at t=-1 which model to use. + """ + F_a, F_b = F_params + G_a, G_b = G_params + + # Single coin flip for the true model + true_models_F = np.random.rand(N) < π_minus_1 + sequences = np.empty((N, T)) + + n_f = np.sum(true_models_F) + n_g = N - n_f + + if n_f > 0: + sequences[true_models_F, :] = np.random.beta(F_a, F_b, (n_f, T)) + if n_g > 0: + sequences[~true_models_F, :] = np.random.beta(G_a, G_b, (n_g, T)) + + return sequences, true_models_F -```{note} -A **rational expectations** model would set $\pi_t^i(s^t) = \pi_t(s^t)$ for all agents $i$. +def protocol_2(π_minus_1, T, N=1000, F_params=(1, 1), G_params=(3, 1.2)): + """ + Simulate Protocol 2: Nature decides at each time step which model to use. + """ + F_a, F_b = F_params + G_a, G_b = G_params + + # Coin flips for each time step + true_models_F = np.random.rand(N, T) < π_minus_1 + sequences = np.empty((N, T)) + + n_f = np.sum(true_models_F) + n_g = N * T - n_f + + if n_f > 0: + sequences[true_models_F] = np.random.beta(F_a, F_b, n_f) + if n_g > 0: + sequences[~true_models_F] = np.random.beta(G_a, G_b, n_g) + + return sequences, true_models_F ``` -There are two agents named $i=1$ and $i=2$. - -At time $t$, agent $1$ receives an endowment - -$$ -y_t^1 = s_t -$$ - -of a nonstorable consumption good, while agent $2$ receives an endowment of - -$$ -y_t^2 = 1 - s_t -$$ +**Remark:** Under timing protocol 2, the $\{w_t\}_{t=1}^T$ is a sequence of IID draws from $h(w)$. Under timing protocol 1, the $\{w_t\}_{t=1}^T$ is +not IID. It is **conditionally IID** -- meaning that with probability $\pi_{-1}$ it is a sequence of IID draws from $f(w)$ and with probability $1-\pi_{-1}$ it is a sequence of IID draws from $g(w)$. For more about this, see {doc}`this lecture about exchangeability `. -The aggregate endowment of the consumption good is +We again deploy a **likelihood ratio process** with time $t$ component being the likelihood ratio $$ -y_t^1 + y_t^2 = 1 +\ell (w_t)=\frac{f\left(w_t\right)}{g\left(w_t\right)},\quad t\geq1. $$ -at each date $t \geq 0$. - -At date $t$ agent $i$ consumes $c_t^i(s^t)$ of the good. - -A (non wasteful) feasible allocation of the aggregate endowment of $1$ each period satisfies +The **likelihood ratio process** for sequence $\left\{ w_{t}\right\} _{t=1}^{\infty}$ is $$ -c_t^1 + c_t^2 = 1 . +L\left(w^{t}\right)=\prod_{i=1}^{t} \ell (w_i), $$ -### A social risk-sharing arrangement - -In order to share risks, a benevolent social planner will dictate a history-dependent consumption allocation in the form of a sequence of functions - -$$ -c_t^i = c_t^i(s^t) -$$ +For shorthand we'll write $L_t = L(w^t)$. -that satisfy +### Model selection mistake probability -$$ -c_t^1(s^t) + c_t^2(s^t) = 1 -$$ (eq:feasibility) +We first study a problem that assumes timing protocol 1. -for all $s^t$ for all $t \geq 0$. +Consider a decision maker who wants to know whether model $f$ or model $g$ governs a data set of length $T$ observations. -To design a socially optimal allocation, the social planner wants to know what agent $1$ believes about the endowment sequence and how they feel about bearing risks. +The decision makers has observed a sequence $\{w_t\}_{t=1}^T$. -As for the endowment sequences, agent $i$ believes that nature draws i.i.d. sequences from joint densities +On the basis of that observed sequence, a likelihood ratio test selects model $f$ when + $L_T \geq 1 $ and model $g$ when $L_T < 1$. + +When model $f$ generates the data, the probability that the likelihood ratio test selects the wrong model is -$$ -\pi_t^i(s^t) = \pi(s_t)^i \pi^i(s_{t-1}) \cdots \pi^i(s_0) $$ - -As for attitudes toward bearing risks, agent $i$ has a one-period utility function - -$$ -u(c_t^i) = \ln (c_t^i) -$$ - -with marginal utility of consumption in period $i$ - -$$ -u'(c_t^i) = \frac{1}{c_t^i} +p_f = {\rm Prob}\left(L_T < 1\Big| f\right) = \alpha_T . $$ -Putting its beliefs about its random endowment sequence and its attitudes toward bearing risks together, agent $i$ has intertemporal utility function +When model $g$ generates the data, the probability that the likelihood ratio test selects the wrong model is +$$ +p_g = {\rm Prob}\left(L_T \geq 1 \Big|g \right) = \beta_T. $$ -V^i = \sum_{t=0}^{\infty} \sum_{s^t} \delta^t u(c_t^i(s^t)) \pi_t^i(s^t) , -$$ (eq:objectiveagenti) -where $\delta \in (0,1)$ is an intertemporal discount factor, and $u(\cdot)$ is a strictly increasing, concave one-period utility function. +We can construct a probability that the likelihood ratio selects the wrong model by assigning a Bayesian prior probability of $\pi_{-1} = .5$ that nature selects model $f$ then averaging $p_f$ and $p_g$ to form the Bayesian posterior probability of a detection error equal to +$$ +p(\textrm{wrong decision}) = {1 \over 2} (\alpha_T + \beta_T) . +$$ (eq:detectionerrorprob) -### The social planner's allocation problem +Now let's simulate timing protocol 1 and 2 and compute the error probabilities -The benevolent dictator has all the information it requires to choose a consumption allocation that maximizes the social welfare criterion +```{code-cell} ipython3 -$$ -W = \lambda V^1 + (1-\lambda) V^2 -$$ (eq:welfareW) +def compute_protocol_1_errors(π_minus_1, T_max, N_simulations, f_func, g_func, + F_params=(1, 1), G_params=(3, 1.2)): + """ + Compute error probabilities for Protocol 1. + """ + sequences, true_models = protocol_1( + π_minus_1, T_max, N_simulations, F_params, G_params) + l_ratios, L_cumulative = compute_likelihood_ratios(sequences, + f_func, g_func) + + T_range = np.arange(1, T_max + 1) + + mask_f = true_models + mask_g = ~true_models + + L_f = L_cumulative[mask_f, :] + L_g = L_cumulative[mask_g, :] + + α_T = np.mean(L_f < 1, axis=0) + β_T = np.mean(L_g >= 1, axis=0) + error_prob = 0.5 * (α_T + β_T) + + return { + 'T_range': T_range, + 'alpha': α_T, + 'beta': β_T, + 'error_prob': error_prob, + 'L_cumulative': L_cumulative, + 'true_models': true_models + } -where $\lambda \in [0,1]$ is a Pareto weight tells how much the planner likes agent $1$ and $1 - \lambda$ is a Pareto weight that tells how much the social planner likes agent $2$. +def compute_protocol_2_errors(π_minus_1, T_max, N_simulations, f_func, g_func, + F_params=(1, 1), G_params=(3, 1.2)): + """ + Compute error probabilities for Protocol 2. + """ + sequences, true_models = protocol_2(π_minus_1, + T_max, N_simulations, F_params, G_params) + l_ratios, _ = compute_likelihood_ratios(sequences, f_func, g_func) + + T_range = np.arange(1, T_max + 1) + + accuracy = np.empty(T_max) + for t in range(T_max): + predictions = (l_ratios[:, t] >= 1) + actual = true_models[:, t] + accuracy[t] = np.mean(predictions == actual) + + return { + 'T_range': T_range, + 'accuracy': accuracy, + 'l_ratios': l_ratios, + 'true_models': true_models + } +``` -Setting $\lambda = .5$ expresses ''egalitarian'' social preferences. +The following code visualizes the error probabilities for timing protocol 1 and 2 -Notice how social welfare criterion {eq}`eq:welfareW` takes into account both agents' preferences as represented by formula {eq}`eq:objectiveagenti`. +```{code-cell} ipython3 +:tags: [hide-input] -This means that the social planner knows and respects +def analyze_protocol_1(π_minus_1, T_max, N_simulations, f_func, g_func, + F_params=(1, 1), G_params=(3, 1.2)): + """Analyze Protocol 1""" + result = compute_protocol_1_errors(π_minus_1, T_max, N_simulations, + f_func, g_func, F_params, G_params) + + # Plot results + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) + + ax1.plot(result['T_range'], result['alpha'], 'b-', + label=r'$\alpha_T$', linewidth=2) + ax1.plot(result['T_range'], result['beta'], 'r-', + label=r'$\beta_T$', linewidth=2) + ax1.set_xlabel('$T$') + ax1.set_ylabel('error probability') + ax1.legend() + + ax2.plot(result['T_range'], result['error_prob'], 'g-', + label=r'$\frac{1}{2}(\alpha_T+\beta_T)$', linewidth=2) + ax2.set_xlabel('$T$') + ax2.set_ylabel('error probability') + ax2.legend() + + plt.tight_layout() + plt.show() + + # Print summary + print(f"At T={T_max}:") + print(f"α_{T_max} = {result['alpha'][-1]:.4f}") + print(f"β_{T_max} = {result['beta'][-1]:.4f}") + print(f"Model selection error probability = {result['error_prob'][-1]:.4f}") + + return result -* each agent's one period utility function $u(\cdot) = \ln(\cdot)$ -* each agent $i$'s probability model $\{\pi_t^i(s^t)\}_{t=0}^\infty$ +def analyze_protocol_2(π_minus_1, T_max, N_simulations, f_func, g_func, + theory_error=None, F_params=(1, 1), G_params=(3, 1.2)): + """Analyze Protocol 2.""" + result = compute_protocol_2_errors(π_minus_1, T_max, N_simulations, + f_func, g_func, F_params, G_params) + + # Plot results + plt.figure(figsize=(10, 6)) + plt.plot(result['T_range'], result['accuracy'], + 'b-', linewidth=2, label='empirical accuracy') + + if theory_error is not None: + plt.axhline(1 - theory_error, color='r', linestyle='--', + label=f'theoretical accuracy = {1 - theory_error:.4f}') + + plt.xlabel('$t$') + plt.ylabel('accuracy') + plt.legend() + plt.ylim(0.5, 1.0) + plt.show() + + return result -Consequently, we anticipate that these objects will appear in the social planner's rule for allocating the aggregate endowment each period. +def compare_protocols(result1, result2): + """Compare results from both protocols.""" + plt.figure(figsize=(10, 6)) + + plt.plot(result1['T_range'], result1['error_prob'], linewidth=2, + label='Protocol 1 (Model Selection)') + plt.plot(result2['T_range'], 1 - result2['accuracy'], + linestyle='--', linewidth=2, + label='Protocol 2 (classification)') + + plt.xlabel('$T$') + plt.ylabel('error probability') + plt.legend() + plt.show() +# Analyze Protocol 1 +π_minus_1 = 0.5 +T_max = 30 +N_simulations = 10_000 -First-order necessary conditions for maximizing welfare criterion {eq}`eq:welfareW` subject to the feasibility constraint {eq}`eq:feasibility` are +result_p1 = analyze_protocol_1(π_minus_1, T_max, N_simulations, + f, g, (F_a, F_b), (G_a, G_b)) +``` -$$\frac{\pi_t^2(s^t)}{\pi_t^1(s^t)} \frac{(1/c_t^2(s^t))}{(1/c_t^1(s^t))} = \frac{\lambda}{1 -\lambda}$$ +Notice how the model selection error probability approaches zero as $T$ grows. -which can be rearranged to become +### Classification +We now consider a problem that assumes timing protocol 2. +A decision maker wants to classify components of an observed sequence $\{w_t\}_{t=1}^T$ as having been drawn from either $f$ or $g$. +The decision maker uses the following classification rule: $$ -\frac{c_t^1(s^t)}{c_t^2(s^t)} = \frac{\lambda}{1- \lambda} l_t(s^t) -$$ (eq:allocationrule0) - +\begin{aligned} +w_t & \ {\rm is \ from \ } f \ {\rm if \ } l_t > 1 \\ +w_t & \ {\rm is \ from \ } g \ {\rm if \ } l_t \leq 1 . +\end{aligned} +$$ -where +Under this rule, the expected misclassification rate is -$$ l_t(s^t) = \frac{\pi_t^1(s^t)}{\pi_t^2(s^t)} $$ +$$ +p(\textrm{misclassification}) = {1 \over 2} (\tilde \alpha_t + \tilde \beta_t) +$$ (eq:classerrorprob) -is the likelihood ratio of agent 1's joint density to agent 2's joint density. +where $\tilde \alpha_t = {\rm Prob}(l_t < 1 \mid f)$ and $\tilde \beta_t = {\rm Prob}(l_t \geq 1 \mid g)$. -Using +Since for each $t$, the decision boundary is the same, the decision boundary can be computed as -$$c_t^1(s^t) + c_t^2(s^t) = 1$$ +```{code-cell} ipython3 +root = brentq(lambda w: f(w) / g(w) - 1, 0.001, 0.999) +``` -we can rewrite allocation rule {eq}`eq:allocationrule0` as +we can plot the distributions of $f$ and $g$ and the decision boundary +```{code-cell} ipython3 +:tags: [hide-input] +fig, ax = plt.subplots(figsize=(7, 6)) -$$\frac{c_t^1(s^t)}{1 - c_t^1(s^t)} = \frac{\lambda}{1-\lambda} l_t(s^t)$$ +w_range = np.linspace(1e-5, 1-1e-5, 1000) +f_values = [f(w) for w in w_range] +g_values = [g(w) for w in w_range] +ratio_values = [f(w)/g(w) for w in w_range] -or +ax.plot(w_range, f_values, 'b-', + label=r'$f(w) \sim Beta(1,1)$', linewidth=2) +ax.plot(w_range, g_values, 'r-', + label=r'$g(w) \sim Beta(3,1.2)$', linewidth=2) -$$c_t^1(s^t) = \frac{\lambda}{1-\lambda} l_t(s^t)(1 - c_t^1(s^t))$$ - -which implies that the social planner's allocation rule is - -$$ -c_t^1(s^t) = \frac{\lambda l_t(s^t)}{1-\lambda + \lambda l_t(s^t)} -$$ (eq:allocationrule1) - -If we define a temporary or **continuation Pareto weight** process as - -$$ -\lambda_t(s^t) = \frac{\lambda l_t(s^t)}{1-\lambda + \lambda l_t(s^t)}, -$$ - -then we can represent the social planner's allocation rule as - -$$ -c_t^1(s^t) = \lambda_t(s^t) . -$$ - - - - -### If you're so smart, $\ldots$ - - -Let's compute some values of limiting allocations {eq}`eq:allocationrule1` for some interesting possible limiting -values of the likelihood ratio process $l_t(s^t)$: - - $$l_\infty (s^\infty)= 1; \quad c_\infty^1 = \lambda$$ - - * In the above case, both agents are equally smart (or equally not smart) and the consumption allocation stays put at a $\lambda, 1 - \lambda$ split between the two agents. - -$$l_\infty (s^\infty) = 0; \quad c_\infty^1 = 0$$ - -* In the above case, agent 2 is ''smarter'' than agent 1, and agent 1's share of the aggregate endowment converges to zero. - - - - -$$l_\infty (s^\infty)= \infty; \quad c_\infty^1 = 1$$ - -* In the above case, agent 1 is smarter than agent 2, and agent 1's share of the aggregate endowment converges to 1. - -```{note} -These three cases are somehow telling us about how relative **wealths** of the agents evolve as time passes. -* when the two agents are equally smart and $\lambda \in (0,1)$, agent 1's wealth share stays at $\lambda$ perpetually. -* when agent 1 is smarter and $\lambda \in (0,1)$, agent 1 eventually "owns" the continuation entire continuation endowment and agent 2 eventually "owns" nothing. -* when agent 2 is smarter and $\lambda \in (0,1)$, agent 2 eventually "owns" the continuation entire continuation endowment and agent 1 eventually "owns" nothing. -Continuation wealths can be defined precisely after we introduce a competitive equilibrium **price** system below. -``` - - -Soon we'll do some simulations that will shed further light on possible outcomes. - -But before we do that, let's take a detour and study some "shadow prices" for the social planning problem that can readily be -converted to "equilibrium prices" for a competitive equilibrium. - -Doing this will allow us to connect our analysis with an argument of {cite}`alchian1950uncertainty` and {cite}`friedman1953essays` that competitive market processes can make prices of risky assets better reflect realistic probability assessments. - - - -### Competitive Equilibrium Prices - -Two fundamental welfare theorems for general equilibrium models lead us to anticipate that there is a connection between the allocation that solves the social planning problem we have been studying and the allocation in a **competitive equilibrium** with complete markets in history-contingent commodities. - -```{note} -For the two welfare theorems and their history, see . -Again, for applications to a classic macroeconomic growth model, see {doc}`this lecture on a planning problem ` and {doc}`this lecture on a related competitive equilibrium ` -``` - -Such a connection prevails for our model. - -We'll sketch it now. - -In a competitive equilibrium, there is no social planner that dictatorially collects everybody's endowments and then reallocates them. - -Instead, there is a comprehensive centralized market that meets at one point in time. - -There are **prices** at which price-taking agents can buy or sell whatever goods that they want. - -Trade is multilateral in the sense that that there is a "Walrasian auctioneer" who lives outside the model and whose job is to verify that -each agent's budget constraint is satisfied. - -That budget constraint involves the total value of the agent's endowment stream and the total value of its consumption stream. - -These values are computed at price vectors that the agents take as given -- they are "price-takers" who assume that they can buy or sell -whatever quantities that they want at those prices. - -Suppose that at time $-1$, before time $0$ starts, agent $i$ can purchase one unit $c_t(s^t)$ of consumption at time $t$ after history -$s^t$ at price $p_t(s^t)$. - -Notice that there is (very long) **vector** of prices. - - * there is one price $p_t(s^t)$ for each history $s^t$ at every date $t = 0, 1, \ldots, $. - * so there are as many prices as there are histories and dates. - -These prices determined at time $-1$ before the economy starts. - -The market meets once at time $-1$. - -At times $t =0, 1, 2, \ldots$ trades made at time $-1$ are executed. - - - -* in the background, there is an "enforcement" procedure that forces agents to carry out the exchanges or "deliveries" that they agreed to at time $-1$. - - - -We want to study how agents' beliefs influence equilibrium prices. - -Agent $i$ faces a **single** intertemporal budget constraint - -$$ -\sum_{t=0}^\infty\sum_{s^t} p_t(s^t) c_t^i (s^t) \leq \sum_{t=0}^\infty\sum_{s^t} p_t(s^t) y_t^i (s^t) -$$ (eq:budgetI) - -According to budget constraint {eq}`eq:budgetI`, trade is **multilateral** in the following sense - -* we can imagine that agent $i$ first sells his random endowment stream $\{y_t^i (s^t)\}$ and then uses the proceeds (i.e., his "wealth") to purchase a random consumption stream $\{c_t^i (s^t)\}$. - -Agent $i$ puts a Lagrange multiplier $\mu_i$ on {eq}`eq:budgetI` and once-and-for-all chooses a consumption plan $\{c^i_t(s^t)\}_{t=0}^\infty$ -to maximize criterion {eq}`eq:objectiveagenti` subject to budget constraint {eq}`eq:budgetI`. - -This means that the agent $i$ chooses many objects, namely, $c_t^i(s^t)$ for all $s^t$ for $t = 0, 1, 2, \ldots$. - - -For convenience, let's remind ourselves of criterion $V^i$ defined in {eq}`eq:objectiveagenti`: - -$$ -V^i = \sum_{t=0}^{\infty} \sum_{s^t} \delta^t u_t(c_t^i(s^t)) \pi_t^i(s^t) -$$ - -First-order necessary conditions for maximizing objective $V^i$ defined in {eq}`eq:objectiveagenti` with respect to $c_t^i(s^t)$ are - -$$ -\delta^t u'(c^i_t(s^t)) \pi_t^i(s^t) = \mu_i p_t(s^t) , -$$ - -which we can rearrange to obtain - -$$ -p_t(s^t) = \frac{ \delta^t \pi_t^i(s^t)}{\mu_i c^i_t(s^t)} -$$ (eq:priceequation1) - -for $i=1,2$. - -If we divide equation {eq}`eq:priceequation1` for agent $1$ by the appropriate version of equation {eq}`eq:priceequation1` for agent 2, use -$c^2_t(s^t) = 1 - c^1_t(s^t)$, and do some algebra, we'll obtain - -$$ -c_t^1(s^t) = \frac{\mu_1 l_t(s^t)}{\mu_2 + \mu_1 l_t(s^t)} . -$$ (eq:allocationce) - -We now engage in an extended "guess-and-verify" exercise that involves matching objects in our competitive equilibrium with objects in -our social planning problem. - -* we'll match consumption allocations in the planning problem with equilibrium consumption allocations in the competitive equilibrium -* we'll match "shadow" prices in the planning problem with competitive equilibrium prices. - -Notice that if we set $\mu_1 = \lambda$ and $\mu_2 = 1 -\lambda$, then formula {eq}`eq:allocationce` agrees with formula -{eq}`eq:allocationrule1`. - - * doing this amounts to choosing a **numeraire** or normalization for the price system $\{p_t(s^t)\}_{t=0}^\infty$ - -```{note} -For information about how a numeraire must be chosen to pin down the absolute price level in a model like ours that determines only -relative prices, see . -``` - -If we substitute formula {eq}`eq:allocationce` for $c_t^1(s^t)$ into formula {eq}`eq:priceequation1` and rearrange, we obtain - -$$ -p_t(s^t) = \frac{\delta^t}{\lambda(1-\lambda)} \pi_t^2(s^t) \bigl[1 - \lambda + \lambda l_t(s^t)\bigr] -$$ - -or - -$$ -p_t(s^t) = \frac{\delta^t}{\lambda(1-\lambda)} \bigl[(1 - \lambda) \pi_t^2(s^t) + \lambda \pi_t^1(s^t)\bigr] -$$ (eq:pformulafinal) - -According to formula {eq}`eq:pformulafinal`, we have the following possible limiting cases: - -* when $l_\infty = 0$, $c_\infty^1 = 0 $ and tails of competitive equilibrium prices reflect agent $2$'s probability model $\pi_t^2(s^t)$ according to $p_t(s^t) \propto \delta^t \pi_t^2(s^t) $ -* when $l_\infty = \infty$, $c_\infty^1 = 1 $ and tails of competitive equilibrium prices reflect agent $1$'s probability model $\pi_t^1(s^t)$ according to $p_t(s^t) \propto \delta^t \pi_t^1(s^t) $ -* for small $t$'s, competitive equilibrium prices reflect both agents' probability models. - -### Simulations - -Now let's implement some simulations when agent $1$ believes marginal density - -$$\pi^1(s_t) = f(s_t) $$ - -and agent $2$ believes marginal density - -$$ \pi^2(s_t) = g(s_t) $$ - -where $f$ and $g$ are Beta distributions like ones that we used in earlier sections of this lecture. - -Meanwhile, we'll assume that nature believes a marginal density - -$$ -\pi(s_t) = h(s_t) -$$ - -where $h(s_t)$ is perhaps a mixture of $f$ and $g$. - -Let's write a Python function that computes agent 1's consumption share - -```{code-cell} ipython3 -def simulate_blume_easley(sequences, f_belief=f, g_belief=g, λ=0.5): - """Simulate Blume-Easley model consumption shares.""" - l_ratios, l_cumulative = compute_likelihood_ratios(sequences, f_belief, g_belief) - c1_share = λ * l_cumulative / (1 - λ + λ * l_cumulative) - return l_cumulative, c1_share -``` - -Now let's use this function to generate sequences in which - -* nature draws from $f$ each period, or -* nature draws from $g$ each period, or -* or nature flips a fair coin each period to decide whether to draw from $f$ or $g$ - -```{code-cell} ipython3 -λ = 0.5 -T = 100 -N = 10000 - -# Nature follows f, g, or mixture -s_seq_f = np.random.beta(F_a, F_b, (N, T)) -s_seq_g = np.random.beta(G_a, G_b, (N, T)) - -h = jit(lambda x: 0.5 * f(x) + 0.5 * g(x)) -model_choices = np.random.rand(N, T) < 0.5 -s_seq_h = np.empty((N, T)) -s_seq_h[model_choices] = np.random.beta(F_a, F_b, size=model_choices.sum()) -s_seq_h[~model_choices] = np.random.beta(G_a, G_b, size=(~model_choices).sum()) - -l_cum_f, c1_f = simulate_blume_easley(s_seq_f) -l_cum_g, c1_g = simulate_blume_easley(s_seq_g) -l_cum_h, c1_h = simulate_blume_easley(s_seq_h) -``` - -Before looking at the figure below, have some fun by guessing whether agent 1 or agent 2 will have a larger and larger consumption share as time passes in our three cases. - -To make better guesses, let's visualize instances of the likelihood ratio processes in the three cases. - -```{code-cell} ipython3 -fig, axes = plt.subplots(2, 3, figsize=(15, 10)) - -titles = ["Nature = f", "Nature = g", "Nature = mixture"] -data_pairs = [(l_cum_f, c1_f), (l_cum_g, c1_g), (l_cum_h, c1_h)] - -for i, ((l_cum, c1), title) in enumerate(zip(data_pairs, titles)): - # Likelihood ratios - ax = axes[0, i] - for j in range(min(50, l_cum.shape[0])): - ax.plot(l_cum[j, :], alpha=0.3, color='blue') - ax.set_yscale('log') - ax.set_xlabel('time') - ax.set_ylabel('Likelihood ratio $l_t$') - ax.set_title(title) - ax.axhline(y=1, color='red', linestyle='--', alpha=0.5) - - # Consumption shares - ax = axes[1, i] - for j in range(min(50, c1.shape[0])): - ax.plot(c1[j, :], alpha=0.3, color='green') - ax.set_xlabel('time') - ax.set_ylabel("Agent 1's consumption share") - ax.set_ylim([0, 1]) - ax.axhline(y=λ, color='red', linestyle='--', alpha=0.5) - -plt.tight_layout() -plt.show() -``` - -In the left panel, nature chooses $f$. Agent 1's consumption reaches $1$ very quickly. - -In the middle panel, nature chooses $g$. Agent 1's consumption ratio tends to move towards $0$ but not as fast as in the first case. - -In the right panel, nature flips coins each period. We see a very similar pattern to the processes in the left panel. - -The figures in the top panel remind us of the discussion in [this section](KL_link). - -We invite readers to revisit [that section](rel_entropy) and try to infer the relationships among $KL(f, g)$, $KL(g, f)$, $KL(h, f)$, and $KL(h,g)$. - - -Let's compute values of KL divergence - -```{code-cell} ipython3 -shares = [np.mean(c1_f[:, -1]), np.mean(c1_g[:, -1]), np.mean(c1_h[:, -1])] -Kf_g, Kg_f = compute_KL(f, g), compute_KL(g, f) -Kf_h, Kg_h = compute_KL_h(h, f, g) - -print(f"Final shares: f={shares[0]:.3f}, g={shares[1]:.3f}, mix={shares[2]:.3f}") -print(f"KL divergences: \nKL(f,g)={Kf_g:.3f}, KL(g,f)={Kg_f:.3f}") -print(f"KL(h,f)={Kf_h:.3f}, KL(h,g)={Kg_h:.3f}") -``` - -We find that $KL(f,g) > KL(g,f)$ and $KL(h,g) > KL(h,f)$. - -The first inequality tells us that the average "surprise" from having belief $g$ when nature chooses $f$ is greater than the "surprise" from having belief $f$ when nature chooses $g$. - -This explains the difference between the first two panels we noted above. - -The second inequality tells us that agent 1's belief distribution $f$ is closer to nature's pick than agent 2's belief $g$. - -+++ - -To make this idea more concrete, let's compare two cases: - -- agent 1's belief distribution $f$ is close to agent 2's belief distribution $g$; -- agent 1's belief distribution $f$ is far from agent 2's belief distribution $g$. - - -We use the two distributions visualized below - -```{code-cell} ipython3 -def plot_distribution_overlap(ax, x_range, f_vals, g_vals, - f_label='f', g_label='g', - f_color='blue', g_color='red'): - """Plot two distributions with their overlap region.""" - ax.plot(x_range, f_vals, color=f_color, linewidth=2, label=f_label) - ax.plot(x_range, g_vals, color=g_color, linewidth=2, label=g_label) - - overlap = np.minimum(f_vals, g_vals) - ax.fill_between(x_range, 0, overlap, alpha=0.3, color='purple', label='Overlap') - ax.set_xlabel('x') - ax.set_ylabel('Density') - ax.legend() - -# Define close and far belief distributions -f_close = jit(lambda x: p(x, 1, 1)) -g_close = jit(lambda x: p(x, 1.1, 1.05)) - -f_far = jit(lambda x: p(x, 1, 1)) -g_far = jit(lambda x: p(x, 3, 1.2)) - -# Visualize the belief distributions -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) - -x_range = np.linspace(0.001, 0.999, 200) - -# Close beliefs -f_close_vals = [f_close(x) for x in x_range] -g_close_vals = [g_close(x) for x in x_range] -plot_distribution_overlap(ax1, x_range, f_close_vals, g_close_vals, - f_label='f (Beta(1, 1))', g_label='g (Beta(1.1, 1.05))') -ax1.set_title(f'Close Beliefs') - -# Far beliefs -f_far_vals = [f_far(x) for x in x_range] -g_far_vals = [g_far(x) for x in x_range] -plot_distribution_overlap(ax2, x_range, f_far_vals, g_far_vals, - f_label='f (Beta(1, 1))', g_label='g (Beta(3, 1.2))') -ax2.set_title(f'Far Beliefs') - -plt.tight_layout() -plt.show() -``` - -Let's draw the same consumption ratio plots as above for agent 1. - -We replace the simulation paths with median and percentiles to make the figure cleaner. - -Staring at the figure below, can we infer the relation between $KL(f,g)$ and $KL(g,f)$? - -From the right panel, can we infer the relation between $KL(h,g)$ and $KL(h,f)$? - -```{code-cell} ipython3 -fig, axes = plt.subplots(2, 3, figsize=(15, 10)) -nature_params = {'close': [(1, 1), (1.1, 1.05), (2, 1.5)], - 'far': [(1, 1), (3, 1.2), (2, 1.5)]} -nature_labels = ["Nature = f", "Nature = g", "Nature = h"] -colors = {'close': 'blue', 'far': 'red'} - -threshold = 1e-5 # "close to zero" cutoff - -for row, (f_belief, g_belief, label) in enumerate([ - (f_close, g_close, 'close'), - (f_far, g_far, 'far')]): - - for col, nature_label in enumerate(nature_labels): - params = nature_params[label][col] - s_seq = np.random.beta(params[0], params[1], (1000, 200)) - _, c1 = simulate_blume_easley(s_seq, f_belief, g_belief, λ) - - median_c1 = np.median(c1, axis=0) - p10, p90 = np.percentile(c1, [10, 90], axis=0) - - ax = axes[row, col] - color = colors[label] - ax.plot(median_c1, color=color, linewidth=2, label='Median') - ax.fill_between(range(len(median_c1)), p10, p90, alpha=0.3, color=color, label='10–90%') - ax.set_xlabel('time') - ax.set_ylabel("Agent 1's share") - ax.set_ylim([0, 1]) - ax.set_title(nature_label) - ax.axhline(y=λ, color='gray', linestyle='--', alpha=0.5) - below = np.where(median_c1 < threshold)[0] - above = np.where(median_c1 > 1-threshold)[0] - if below.size > 0: first_zero = (below[0], True) - elif above.size > 0: first_zero = (above[0], False) - else: first_zero = None - if first_zero is not None: - ax.axvline(x=first_zero[0], color='black', linestyle='--', - alpha=0.7, - label=fr'Median $\leq$ {threshold}' if first_zero[1] - else fr'Median $\geq$ 1-{threshold}') - ax.legend() - -plt.tight_layout() -plt.show() -``` - -Holding to our guesses, let's calculate the four values - -```{code-cell} ipython3 -# Close case -Kf_g, Kg_f = compute_KL(f_close, g_close), compute_KL(g_close, f_close) -Kf_h, Kg_h = compute_KL_h(h, f_close, g_close) - -print(f"KL divergences (close): \nKL(f,g)={Kf_g:.3f}, KL(g,f)={Kg_f:.3f}") -print(f"KL(h,f)={Kf_h:.3f}, KL(h,g)={Kg_h:.3f}") - -# Far case -Kf_g, Kg_f = compute_KL(f_far, g_far), compute_KL(g_far, f_far) -Kf_h, Kg_h = compute_KL_h(h, f_far, g_far) - -print(f"KL divergences (far): \nKL(f,g)={Kf_g:.3f}, KL(g,f)={Kg_f:.3f}") -print(f"KL(h,f)={Kf_h:.3f}, KL(h,g)={Kg_h:.3f}") -``` - -We find that in the first case, $KL(f,g) \approx KL(g,f)$ and both are relatively small, so although either agent 1 or agent 2 will eventually consume everything, convergence displaying in first two panels on the top is pretty slowly. - -In the first two panels at the bottom, we see convergence occurring faster (as indicated by the black dashed line) because the divergence gaps $KL(f, g)$ and $KL(g, f)$ are larger. - -Since $KL(f,g) > KL(g,f)$, we see faster convergence in the first panel at the bottom when nature chooses $f$ than in the second panel where nature chooses $g$. - -This ties in nicely with {eq}`eq:kl_likelihood_link`. - -## Hypothesis Testing and Classification - -This section discusses another application of likelihood ratio processes. - -We describe how a statistician can combine frequentist probabilities of type I and type II errors in order to - -* compute an anticipated frequency of selecting a wrong model based on a sample length $T$ -* compute an anticipated error rate in a classification problem - -We consider a situation in which nature generates data by mixing known densities $f$ and $g$ with known mixing -parameter $\pi_{-1} \in (0,1)$ so that the random variable $w$ is drawn from the density - -$$ -h (w) = \pi_{-1} f(w) + (1-\pi_{-1}) g(w) -$$ - -We assume that the statistician knows the densities $f$ and $g$ and also the mixing parameter $\pi_{-1}$. - -Below, we'll set $\pi_{-1} = .5$, although much of the analysis would follow through with other settings of $\pi_{-1} \in (0,1)$. - -We assume that $f$ and $g$ both put positive probabilities on the same intervals of possible realizations of the random variable $W$. - - - -In the simulations below, we specify that $f$ is a $\text{Beta}(1, 1)$ distribution and that $g$ is $\text{Beta}(3, 1.2)$ distribution. - -We consider two alternative timing protocols. - - * Timing protocol 1 is for the model selection problem - * Timing protocol 2 is for the individual classification problem - -**Timing Protocol 1:** Nature flips a coin only **once** at time $t=-1$ and with probability $\pi_{-1}$ generates a sequence $\{w_t\}_{t=1}^T$ -of IID draws from $f$ and with probability $1-\pi_{-1}$ generates a sequence $\{w_t\}_{t=1}^T$ -of IID draws from $g$. - -Let's write some Python code that implements timing protocol 1. - -```{code-cell} ipython3 -def protocol_1(π_minus_1, T, N=1000): - """ - Simulate Protocol 1: - Nature decides once at t=-1 which model to use. - """ - - # On-off coin flip for the true model - true_models_F = np.random.rand(N) < π_minus_1 - - sequences = np.empty((N, T)) - - n_f = np.sum(true_models_F) - n_g = N - n_f - if n_f > 0: - sequences[true_models_F, :] = np.random.beta(F_a, F_b, (n_f, T)) - if n_g > 0: - sequences[~true_models_F, :] = np.random.beta(G_a, G_b, (n_g, T)) - - return sequences, true_models_F -``` - -**Timing Protocol 2.** Nature flips a coin **often**. At each time $t \geq 0$, nature flips a coin and with probability $\pi_{-1}$ draws $w_t$ from $f$ and with probability $1-\pi_{-1}$ draws $w_t$ from $g$. - -Here is Python code that we'll use to implement timing protocol 2. - -```{code-cell} ipython3 -def protocol_2(π_minus_1, T, N=1000): - """ - Simulate Protocol 2: - Nature decides at each time step which model to use. - """ - - # Coin flips for each time t upto T - true_models_F = np.random.rand(N, T) < π_minus_1 - - sequences = np.empty((N, T)) - - n_f = np.sum(true_models_F) - n_g = N * T - n_f - if n_f > 0: - sequences[true_models_F] = np.random.beta(F_a, F_b, n_f) - if n_g > 0: - sequences[~true_models_F] = np.random.beta(G_a, G_b, n_g) - - return sequences, true_models_F -``` - -**Remark:** Under timing protocol 2, the $\{w_t\}_{t=1}^T$ is a sequence of IID draws from $h(w)$. Under timing protocol 1, the $\{w_t\}_{t=1}^T$ is -not IID. It is **conditionally IID** -- meaning that with probability $\pi_{-1}$ it is a sequence of IID draws from $f(w)$ and with probability $1-\pi_{-1}$ it is a sequence of IID draws from $g(w)$. For more about this, see {doc}`this lecture about exchangeability `. - -We again deploy a **likelihood ratio process** with time $t$ component being the likelihood ratio - -$$ -\ell (w_t)=\frac{f\left(w_t\right)}{g\left(w_t\right)},\quad t\geq1. -$$ - -The **likelihood ratio process** for sequence $\left\{ w_{t}\right\} _{t=1}^{\infty}$ is - -$$ -L\left(w^{t}\right)=\prod_{i=1}^{t} \ell (w_i), -$$ - -For shorthand we'll write $L_t = L(w^t)$. - -### Model Selection Mistake Probability - -We first study a problem that assumes timing protocol 1. - -Consider a decision maker who wants to know whether model $f$ or model $g$ governs a data set of length $T$ observations. - -The decision makers has observed a sequence $\{w_t\}_{t=1}^T$. - -On the basis of that observed sequence, a likelihood ratio test selects model $f$ when - $L_T \geq 1 $ and model $g$ when $L_T < 1$. - -When model $f$ generates the data, the probability that the likelihood ratio test selects the wrong model is - -$$ -p_f = {\rm Prob}\left(L_T < 1\Big| f\right) = \alpha_T . -$$ - -When model $g$ generates the data, the probability that the likelihood ratio test selects the wrong model is - -$$ -p_g = {\rm Prob}\left(L_T \geq 1 \Big|g \right) = \beta_T. -$$ - -We can construct a probability that the likelihood ratio selects the wrong model by assigning a Bayesian prior probability of $\pi_{-1} = .5$ that nature selects model $f$ then averaging $p_f$ and $p_g$ to form the Bayesian posterior probability of a detection error equal to - -$$ -p(\textrm{wrong decision}) = {1 \over 2} (\alpha_T + \beta_T) . -$$ (eq:detectionerrorprob) - -Now let's simulate timing protocol 1 and compute the error probabilities - -```{code-cell} ipython3 -# Set parameters -π_minus_1 = 0.5 -T_max = 30 -N_simulations = 10_000 - -sequences_p1, true_models_p1 = protocol_1( - π_minus_1, T_max, N_simulations) -l_ratios_p1, L_cumulative_p1 = compute_likelihood_ratios(sequences_p1, f, g) - -# Compute error probabilities for different sample sizes -T_range = np.arange(1, T_max + 1) - -# Boolean masks for true models -mask_f = true_models_p1 -mask_g = ~true_models_p1 - -# Select cumulative likelihoods for each model -L_f = L_cumulative_p1[mask_f, :] -L_g = L_cumulative_p1[mask_g, :] - -α_T = np.mean(L_f < 1, axis=0) -β_T = np.mean(L_g >= 1, axis=0) - -error_prob = 0.5 * (α_T + β_T) - -# Plot results -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) - -ax1.plot(T_range, α_T, 'b-', - label=r'$\alpha_T$', linewidth=2) -ax1.plot(T_range, β_T, 'r-', - label=r'$\beta_T$', linewidth=2) -ax1.set_xlabel('$T$') -ax1.set_ylabel('error probability') -ax1.legend() - -ax2.plot(T_range, error_prob, 'g-', - label=r'$\frac{1}{2}(\alpha_T+\beta_T)$', linewidth=2) -ax2.set_xlabel('$T$') -ax2.set_ylabel('error probability') -ax2.legend() - -plt.tight_layout() -plt.show() - -print(f"At T={T_max}:") -print(f"α_{T_max} = {α_T[-1]:.4f}") -print(f"β_{T_max} = {β_T[-1]:.4f}") -print(f"Model selection error probability = {error_prob[-1]:.4f}") -``` - -Notice how the model selection error probability approaches zero as $T$ grows. - -### Classification - -We now consider a problem that assumes timing protocol 2. - -A decision maker wants to classify components of an observed sequence $\{w_t\}_{t=1}^T$ as having been drawn from either $f$ or $g$. - -The decision maker uses the following classification rule: - -$$ -\begin{aligned} -w_t & \ {\rm is \ from \ f \ if \ } l_t > 1 \\ -w_t & \ {\rm is \ from \ g \ if \ } l_t \leq 1 . -\end{aligned} -$$ - -Under this rule, the expected misclassification rate is - -$$ -p(\textrm{misclassification}) = {1 \over 2} (\tilde \alpha_t + \tilde \beta_t) -$$ (eq:classerrorprob) - -where $\tilde \alpha_t = {\rm Prob}(l_t < 1 \mid f)$ and $\tilde \beta_t = {\rm Prob}(l_t \geq 1 \mid g)$. - -Since for each $t$, the decision boundary is the same, the decision boundary can be computed as - -```{code-cell} ipython3 -root = brentq(lambda w: f(w) / g(w) - 1, 0.001, 0.999) -``` - -we can plot the distributions of $f$ and $g$ and the decision boundary - -```{code-cell} ipython3 -:tags: [hide-input] - -fig, ax = plt.subplots(figsize=(7, 6)) - -w_range = np.linspace(1e-5, 1-1e-5, 1000) -f_values = [f(w) for w in w_range] -g_values = [g(w) for w in w_range] -ratio_values = [f(w)/g(w) for w in w_range] - -ax.plot(w_range, f_values, 'b-', - label=r'$f(w) \sim Beta(1,1)$', linewidth=2) -ax.plot(w_range, g_values, 'r-', - label=r'$g(w) \sim Beta(3,1.2)$', linewidth=2) - -type1_prob = 1 - beta_dist.cdf(root, F_a, F_b) -type2_prob = beta_dist.cdf(root, G_a, G_b) +type1_prob = 1 - beta_dist.cdf(root, F_a, F_b) +type2_prob = beta_dist.cdf(root, G_a, G_b) w_type1 = w_range[w_range >= root] f_type1 = [f(w) for w in w_type1] @@ -1718,41 +1213,16 @@ Now we simulate timing protocol 2 and compute the classification error probabili In the next cell, we also compare the theoretical classification accuracy to the empirical classification accuracy ```{code-cell} ipython3 -accuracy = np.empty(T_max) - -sequences_p2, true_sources_p2 = protocol_2( - π_minus_1, T_max, N_simulations) -l_ratios_p2, _ = compute_likelihood_ratios(sequences_p2, f, g) - -for t in range(T_max): - predictions = (l_ratios_p2[:, t] >= 1) - actual = true_sources_p2[:, t] - accuracy[t] = np.mean(predictions == actual) - -plt.figure(figsize=(10, 6)) -plt.plot(range(1, T_max + 1), accuracy, - 'b-', linewidth=2, label='empirical accuracy') -plt.axhline(1 - theory_error, color='r', linestyle='--', - label=f'theoretical accuracy = {1 - theory_error:.4f}') -plt.xlabel('$t$') -plt.ylabel('accuracy') -plt.legend() -plt.ylim(0.5, 1.0) -plt.show() +# Analyze Protocol 2 +result_p2 = analyze_protocol_2(π_minus_1, T_max, N_simulations, f, g, + theory_error, (F_a, F_b), (G_a, G_b)) ``` Let's watch decisions made by the two timing protocols as more and more observations accrue. ```{code-cell} ipython3 -fig, ax = plt.subplots(figsize=(7, 6)) - -ax.plot(T_range, error_prob, linewidth=2, - label='Protocol 1') -ax.plot(T_range, 1-accuracy, linestyle='--', linewidth=2, - label=f'Protocol 2') -ax.set_ylabel('error probability') -ax.legend() -plt.show() +# Compare both protocols +compare_protocols(result_p1, result_p2) ``` From the figure above, we can see: @@ -1765,279 +1235,15 @@ From the figure above, we can see: **Remark:** Think about how laws of large numbers are applied to compute error probabilities for the model selection problem and the classification problem. -## Special case: Markov chain models - -Consider two $n$-state irreducible and aperiodic Markov chain models on the same state space $\{1, 2, \ldots, n\}$ with positive transition matrices $P^{(f)}$, $P^{(g)}$ and initial distributions $\pi_0^{(f)}$, $\pi_0^{(g)}$. - -In this section, we assume nature chooses $f$. - -For a sample path $(x_0, x_1, \ldots, x_T)$, let $N_{ij}$ count transitions from state $i$ to $j$. - -The likelihood under model $m \in \{f, g\}$ is - -$$ -L_T^{(m)} = \pi_{0,x_0}^{(m)} \prod_{i=1}^n \prod_{j=1}^n \left(P_{ij}^{(m)}\right)^{N_{ij}} -$$ - -Hence, - -$$ -\log L_T^{(m)} =\log\pi_{0,x_0}^{(m)} +\sum_{i,j}N_{ij}\log P_{ij}^{(m)} -$$ - -The log-likelihood ratio is - -$$ -\log \frac{L_T^{(f)}}{L_T^{(g)}} = \log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} + \sum_{i,j}N_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} -$$ (eq:llr_markov) - -### KL divergence rate - -By the ergodic theorem for irreducible, aperiodic Markov chains, we have - -$$ -\frac{N_{ij}}{T} \xrightarrow{a.s.} \pi_i^{(f)}P_{ij}^{(f)} \quad \text{as } T \to \infty -$$ - -where $\boldsymbol{\pi}^{(f)}$ is the stationary distribution satisfying $\boldsymbol{\pi}^{(f)} = \boldsymbol{\pi}^{(f)} P^{(f)}$. - -Therefore, - -$$ -\frac{1}{T}\log \frac{L_T^{(f)}}{L_T^{(g)}} = \frac{1}{T}\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} + \frac{1}{T}\sum_{i,j}N_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} -$$ - -Taking the limit as $T \to \infty$, we have -- The first term: $\frac{1}{T}\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} \to 0$ -- The second term: $\frac{1}{T}\sum_{i,j}N_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} \xrightarrow{a.s.} \sum_{i,j}\pi_i^{(f)}P_{ij}^{(f)}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}}$ - -Define the **KL divergence rate** as - -$$ -h_{KL}(f, g) = \sum_{i=1}^n \pi_i^{(f)} \underbrace{\sum_{j=1}^n P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}}}_{=: KL(P_{i\cdot}^{(f)}, P_{i\cdot}^{(g)})} -$$ - -where $KL(P_{i\cdot}^{(f)}, P_{i\cdot}^{(g)})$ is the row-wise KL divergence. - - -By the ergodic theorem, we have - -$$ -\frac{1}{T}\log \frac{L_T^{(f)}}{L_T^{(g)}} \xrightarrow{a.s.} h_{KL}(f, g) \quad \text{as } T \to \infty -$$ - -Taking expectations and using the dominated convergence theorem, we can get - -$$ -\frac{1}{T}E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] \to h_{KL}(f, g) \quad \text{as } T \to \infty -$$ - -Here we invite readers to pause and compare this result with {eq}`eq:kl_likelihood_link`. - -Let's confirm this in the simulation below. - -### Simulations - -Let's implement simulations to illustrate these concepts with a three-state Markov chain. - -We start with writing out functions to compute the stationary distribution and the KL divergence rate for Markov chain models - -```{code-cell} ipython3 -def compute_stationary_dist(P): - """ - Compute stationary distribution of transition matrix P - """ - - eigenvalues, eigenvectors = np.linalg.eig(P.T) - idx = np.argmax(np.abs(eigenvalues)) - stationary = np.real(eigenvectors[:, idx]) - return stationary / stationary.sum() - -def markov_kl_divergence(P_f, P_g, pi_f): - """ - Compute KL divergence rate between two Markov chains - """ - if np.any((P_f > 0) & (P_g == 0)): - return np.inf - - valid_mask = (P_f > 0) & (P_g > 0) - - log_ratios = np.zeros_like(P_f) - log_ratios[valid_mask] = np.log(P_f[valid_mask] / P_g[valid_mask]) - - # Weight by stationary probabilities and sum - kl_rate = np.sum(pi_f[:, np.newaxis] * P_f * log_ratios) - - return kl_rate - -def simulate_markov_chain(P, pi_0, T, N_paths=1000): - """ - Simulate N_paths sample paths from a Markov chain - """ - mc = qe.MarkovChain(P, state_values=None) - - initial_states = np.random.choice(len(P), size=N_paths, p=pi_0) - - paths = np.zeros((N_paths, T+1), dtype=int) - - for i in range(N_paths): - path = mc.simulate(T+1, init=initial_states[i]) - paths[i, :] = path - - return paths - - -def compute_likelihood_ratio_markov(paths, P_f, P_g, π_0_f, π_0_g): - """ - Compute likelihood ratio process for Markov chain paths - """ - N_paths, T_plus_1 = paths.shape - T = T_plus_1 - 1 - L_ratios = np.ones((N_paths, T+1)) - - L_ratios[:, 0] = π_0_f[paths[:, 0]] / π_0_g[paths[:, 0]] - - for t in range(1, T+1): - prev_states = paths[:, t-1] - curr_states = paths[:, t] - - transition_ratios = P_f[prev_states, curr_states] \ - / P_g[prev_states, curr_states] - L_ratios[:, t] = L_ratios[:, t-1] * transition_ratios - - return L_ratios -``` - -Now let's create an example with two different 3-state Markov chains - -```{code-cell} ipython3 -P_f = np.array([[0.7, 0.2, 0.1], - [0.3, 0.5, 0.2], - [0.1, 0.3, 0.6]]) - -P_g = np.array([[0.5, 0.3, 0.2], - [0.2, 0.6, 0.2], - [0.2, 0.2, 0.6]]) - -# Compute stationary distributions -π_f = compute_stationary_dist(P_f) -π_g = compute_stationary_dist(P_g) - -print(f"Stationary distribution (f): {π_f}") -print(f"Stationary distribution (g): {π_g}") - -# Compute KL divergence rate -kl_rate_fg = markov_kl_divergence(P_f, P_g, π_f) -kl_rate_gf = markov_kl_divergence(P_g, P_f, π_g) - -print(f"\nKL divergence rate h(f, g): {kl_rate_fg:.4f}") -print(f"KL divergence rate h(g, f): {kl_rate_gf:.4f}") -``` - -We are now ready to simulate paths and visualize how likelihood ratios evolve. - -We'll verify $\frac{1}{T}E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = h_{KL}(f, g)$ starting from the stationary distribution by plotting both the empirical average and the line predicted by the theory - -```{code-cell} ipython3 -T = 500 -N_paths = 1000 -paths_from_f = simulate_markov_chain(P_f, π_f, T, N_paths) - -L_ratios_f = compute_likelihood_ratio_markov(paths_from_f, - P_f, P_g, - π_f, π_g) - -plt.figure(figsize=(10, 6)) - -# Plot individual paths -n_show = 50 -for i in range(n_show): - plt.plot(np.log(L_ratios_f[i, :]), - alpha=0.3, color='blue', lw=0.8) - -# Compute theoretical expectation -theory_line = kl_rate_fg * np.arange(T+1) -plt.plot(theory_line, 'k--', linewidth=2.5, - label=r'$T \times h_{KL}(f,g)$') - -# Compute empirical mean -avg_log_L = np.mean(np.log(L_ratios_f), axis=0) -plt.plot(avg_log_L, 'r-', linewidth=2.5, - label='empirical average', alpha=0.5) - -plt.axhline(y=0, color='gray', - linestyle='--', alpha=0.5) -plt.xlabel(r'$T$') -plt.ylabel(r'$\log L_T$') -plt.title('nature = $f$') -plt.legend() -plt.show() -``` - -Let's examine how the model selection error probability depends on sample size using the same simulation strategy in the previous section - -```{code-cell} ipython3 -def compute_selection_error(T_values, P_f, P_g, π_0_f, π_0_g, N_sim=1000): - """ - Compute model selection error probability for different sample sizes - """ - errors = [] - - for T in T_values: - # Simulate from both models - paths_f = simulate_markov_chain(P_f, π_0_f, T, N_sim//2) - paths_g = simulate_markov_chain(P_g, π_0_g, T, N_sim//2) - - # Compute likelihood ratios - L_f = compute_likelihood_ratio_markov(paths_f, - P_f, P_g, - π_0_f, π_0_g) - L_g = compute_likelihood_ratio_markov(paths_g, - P_f, P_g, - π_0_f, π_0_g) - - # Decision rule: choose f if L_T >= 1 - error_f = np.mean(L_f[:, -1] < 1) # Type I error - error_g = np.mean(L_g[:, -1] >= 1) # Type II error - - total_error = 0.5 * (error_f + error_g) - errors.append(total_error) - - return np.array(errors) - -# Compute error probabilities -T_values = np.arange(10, 201, 10) -errors = compute_selection_error(T_values, - P_f, P_g, - π_f, π_g) - -# Plot results -plt.figure(figsize=(10, 6)) -plt.plot(T_values, errors, linewidth=2) -plt.xlabel('$T$') -plt.ylabel('error probability') -plt.show() -``` - -## Measuring discrepancies between distributions +### Error probability and divergence measures A plausible guess is that the ability of a likelihood ratio to distinguish distributions $f$ and $g$ depends on how "different" they are. -But how should we measure discrepancies between distributions? - -We've already encountered one discrepancy measure -- the Kullback-Leibler (KL) divergence. +We have learnt some measures of "difference" between distributions in {doc}`divergence_measures`. -We now briefly explore two alternative discrepancy measures. - -### Chernoff entropy - -Chernoff entropy was motivated by an early application of the [theory of large deviations](https://en.wikipedia.org/wiki/Large_deviations_theory). - -```{note} -Large deviation theory provides refinements of the central limit theorem. -``` +Let's now study two more measures of "difference" between distributions that are useful in the context of model selection and classification. -The Chernoff entropy between probability densities $f$ and $g$ is defined as: +Recall that Chernoff entropy between probability densities $f$ and $g$ is defined as: $$ C(f,g) = - \log \min_{\phi \in (0,1)} \int f^\phi(x) g^{1-\phi}(x) dx @@ -2049,8 +1255,6 @@ $$ e^{-C(f,g)T} . $$ -Thus, Chernoff entropy is an upper bound on the exponential rate at which the selection error probability falls as sample size $T$ grows. - Let's compute Chernoff entropy numerically with some Python code ```{code-cell} ipython3 @@ -2060,7 +1264,7 @@ def chernoff_integrand(ϕ, f, g): """ def integrand(w): return f(w)**ϕ * g(w)**(1-ϕ) - + result, _ = quad(integrand, 1e-5, 1-1e-5) return result @@ -2073,7 +1277,6 @@ def compute_chernoff_entropy(f, g): # Find the minimum over ϕ in (0,1) result = minimize_scalar(objective, - # For numerical stability bounds=(1e-5, 1-1e-5), method='bounded') min_value = result.fun @@ -2081,7 +1284,6 @@ def compute_chernoff_entropy(f, g): chernoff_entropy = -np.log(min_value) return chernoff_entropy, ϕ_optimal - C_fg, ϕ_optimal = compute_chernoff_entropy(f, g) print(f"Chernoff entropy C(f,g) = {C_fg:.4f}") print(f"Optimal ϕ = {ϕ_optimal:.4f}") @@ -2098,7 +1300,7 @@ fig, ax = plt.subplots(figsize=(10, 6)) ax.semilogy(T_range, chernoff_bound, 'r-', linewidth=2, label=f'$e^{{-C(f,g)T}}$') -ax.semilogy(T_range, error_prob, 'b-', linewidth=2, +ax.semilogy(T_range, result_p1['error_prob'], 'b-', linewidth=2, label='Model selection error probability') ax.set_xlabel('T') @@ -2110,19 +1312,13 @@ plt.show() Evidently, $e^{-C(f,g)T}$ is an upper bound on the error rate. -### Jensen-Shannon divergence - -The [Jensen-Shannon divergence](https://en.wikipedia.org/wiki/Jensen%E2%80%93Shannon_divergence) is another divergence measure. - -For probability densities $f$ and $g$, the **Jensen-Shannon divergence** is defined as: - -$$ -D(f,g) = \frac{1}{2} KL(f, m) + \frac{1}{2} KL(g, m) -$$ (eq:compute_JS) +In `{doc}`divergence_measures`, we also studied **Jensen-Shannon divergence** as +a symmetric measure of distance between distributions. -where $m = \frac{1}{2}(f+g)$ is a mixture of $f$ and $g$. +We can use Jensen-Shannon divergence to measure the distance between distributions $f$ and $g$ and +compute how it covaries with the model selection error probability. -Below we compute Jensen-Shannon divergence numerically with some Python code +We also compute Jensen-Shannon divergence numerically with some Python code ```{code-cell} ipython3 def compute_JS(f, g): @@ -2136,24 +1332,19 @@ def compute_JS(f, g): return js_div ``` - -```{note} -We studied KL divergence in the [section above](rel_entropy) with respect to a reference distribution $h$. - -Recall that KL divergence $KL(f, g)$ measures expected excess surprisal from using misspecified model $g$ instead $f$ when $f$ is the true model. +Now let's return to our guess that the error probability at large sample sizes is related to the Chernoff entropy between two distributions. -Because in general $KL(f, g) \neq KL(g, f)$, KL divergence is not symmetric, but Jensen-Shannon divergence is symmetric. +We verify this by computing the correlation between the log of the error probability at $T=50$ under Timing Protocol 1 and the divergence measures. -(In fact, the square root of the Jensen-Shannon divergence is a metric referred to as the Jensen-Shannon distance.) +In the simulation below, nature draws $N / 2$ sequences from $g$ and $N/2$ sequences from $f$. -As {eq}`eq:compute_JS` shows, the Jensen-Shannon divergence computes average of the KL divergence of $f$ and $g$ with respect to a particular reference distribution $m$ defined below the equation. +```{note} +Nature does this rather than flipping a fair coin to decide whether to draw from $g$ or $f$ once and for all before each simulation of length $T$. ``` -Now let's create a comparison table showing KL divergence, Jensen-Shannon divergence, and Chernoff entropy for a set of pairs of Beta distributions. +We use the following pairs of Beta distributions for $f$ and $g$ as test cases ```{code-cell} ipython3 -:tags: [hide-input] - distribution_pairs = [ # (f_params, g_params) ((1, 1), (0.1, 0.2)), @@ -2167,169 +1358,13 @@ distribution_pairs = [ ((1, 1), (1.5, 1.2)), ((1, 1), (2, 1.5)), ((1, 1), (2.5, 1.8)), - ((1, 1), (3, 1.2)), - ((1, 1), (4, 1)), - ((1, 1), (5, 1)) -] - -# Create comparison table -results = [] -for i, ((f_a, f_b), (g_a, g_b)) in enumerate(distribution_pairs): - # Define the density functions - f = jit(lambda x, a=f_a, b=f_b: p(x, a, b)) - g = jit(lambda x, a=g_a, b=g_b: p(x, a, b)) - - # Compute measures - kl_fg = compute_KL(f, g) - kl_gf = compute_KL(g, f) - js_div = compute_JS(f, g) - chernoff_ent, _ = compute_chernoff_entropy(f, g) - - results.append({ - 'Pair (f, g)': f"\\text{{Beta}}({f_a},{f_b}), \\text{{Beta}}({g_a},{g_b})", - 'KL(f, g)': f"{kl_fg:.4f}", - 'KL(g, f)': f"{kl_gf:.4f}", - 'JS': f"{js_div:.4f}", - 'C': f"{chernoff_ent:.4f}" - }) - -df = pd.DataFrame(results) - -# Sort by JS divergence -df['JS_numeric'] = df['JS'].astype(float) -df = df.sort_values('JS_numeric').drop('JS_numeric', axis=1) - -# Generate LaTeX table manually -columns = ' & '.join([f'\\text{{{col}}}' for col in df.columns]) -rows = ' \\\\\n'.join( - [' & '.join([f'{val}' for val in row]) - for row in df.values]) - -latex_code = rf""" -\begin{{array}}{{lcccc}} -{columns} \\ -\hline -{rows} -\end{{array}} -""" - -display(Math(latex_code)) -``` - -The above table indicates how Jensen-Shannon divergence, and Chernoff entropy, and KL divergence covary as we alter $f$ and $g$. - -Let's also visualize how these diverge measures covary - -```{code-cell} ipython3 -kl_fg_values = [float(result['KL(f, g)']) for result in results] -js_values = [float(result['JS']) for result in results] -chernoff_values = [float(result['C']) for result in results] - -fig, axes = plt.subplots(1, 2, figsize=(12, 5)) - -# JS divergence and KL divergence -axes[0].scatter(kl_fg_values, js_values, alpha=0.7, s=60) -axes[0].set_xlabel('KL divergence KL(f, g)') -axes[0].set_ylabel('JS divergence') -axes[0].set_title('JS divergence and KL divergence') - -# Chernoff Entropy and JS divergence -axes[1].scatter(js_values, chernoff_values, alpha=0.7, s=60) -axes[1].set_xlabel('JS divergence') -axes[1].set_ylabel('Chernoff entropy') -axes[1].set_title('Chernoff entropy and JS divergence') - -plt.tight_layout() -plt.show() -``` - -To make the comparison more concrete, let's plot the distributions and the divergence measures for a few pairs of distributions. - -Note that the numbers on the title changes with the area of the overlaps of two distributions - -```{code-cell} ipython3 -:tags: [hide-input] - -def plot_dist_diff(): - """ - Plot overlap of two distributions and divergence measures - """ - - # Chose a subset of Beta distribution parameters - param_grid = [ - ((1, 1), (1, 1)), - ((1, 1), (1.5, 1.2)), - ((1, 1), (2, 1.5)), - ((1, 1), (3, 1.2)), - ((1, 1), (5, 1)), - ((1, 1), (0.3, 0.3)) - ] - - fig, axes = plt.subplots(3, 2, figsize=(15, 12)) - - divergence_data = [] - - for i, ((f_a, f_b), (g_a, g_b)) in enumerate(param_grid): - row = i // 2 - col = i % 2 - - # Create density functions - f = jit(lambda x, a=f_a, b=f_b: p(x, a, b)) - g = jit(lambda x, a=g_a, b=g_b: p(x, a, b)) - - # Compute divergence measures - kl_fg = compute_KL(f, g) - js_div = compute_JS(f, g) - chernoff_ent, _ = compute_chernoff_entropy(f, g) - - divergence_data.append({ - 'f_params': (f_a, f_b), - 'g_params': (g_a, g_b), - 'kl_fg': kl_fg, - 'js_div': js_div, - 'chernoff': chernoff_ent - }) - - # Plot distributions - x_range = np.linspace(0, 1, 200) - f_vals = [f(x) for x in x_range] - g_vals = [g(x) for x in x_range] - - axes[row, col].plot(x_range, f_vals, 'b-', linewidth=2, - label=f'f ~ Beta({f_a},{f_b})') - axes[row, col].plot(x_range, g_vals, 'r-', linewidth=2, - label=f'g ~ Beta({g_a},{g_b})') - - # Fill overlap region - overlap = np.minimum(f_vals, g_vals) - axes[row, col].fill_between(x_range, 0, overlap, alpha=0.3, - color='purple', label='overlap') - - # Add divergence information - axes[row, col].set_title( - f'KL(f, g)={kl_fg:.3f}, JS={js_div:.3f}, C={chernoff_ent:.3f}', - fontsize=12) - axes[row, col].legend(fontsize=14) - - plt.tight_layout() - plt.show() - - return divergence_data - -divergence_data = plot_dist_diff() + ((1, 1), (3, 1.2)), + ((1, 1), (4, 1)), + ((1, 1), (5, 1)) +] ``` -### Error probability and divergence measures - -Now let's return to our guess that the error probability at large sample sizes is related to the Chernoff entropy between two distributions. - -We verify this by computing the correlation between the log of the error probability at $T=50$ under Timing Protocol 1 and the divergence measures. - -In the simulation below, nature draws $N / 2$ sequences from $g$ and $N/2$ sequences from $f$. - -```{note} -Nature does this rather than flipping a fair coin to decide whether to draw from $g$ or $f$ once and for all before each simulation of length $T$. -``` +Now let's run the simmulation ```{code-cell} ipython3 # Parameters for simulation @@ -2431,14 +1466,258 @@ Evidently, Chernoff entropy and Jensen-Shannon entropy each covary tightly with We'll encounter related ideas in {doc}`wald_friedman` very soon. +(lrp_markov)= +## Markov chains + +Let's now look at a likelihood ratio process for a sequence of random variables that is not independently and identically distributed. + +Here we assume that the sequence is generated by a Markov chain on a finite state space. + +We consider two $n$-state irreducible and aperiodic Markov chain models on the same state space $\{1, 2, \ldots, n\}$ with positive transition matrices $P^{(f)}$, $P^{(g)}$ and initial distributions $\pi_0^{(f)}$, $\pi_0^{(g)}$. + +We assume that nature samples from chain $f$. + +For a sample path $(x_0, x_1, \ldots, x_T)$, let $N_{ij}$ count transitions from state $i$ to $j$. + +The likelihood process under model $m \in \{f, g\}$ is + +$$ +L_T^{(m)} = \pi_{0,x_0}^{(m)} \prod_{i=1}^n \prod_{j=1}^n \left(P_{ij}^{(m)}\right)^{N_{ij}} +$$ + +Hence, + +$$ +\log L_T^{(m)} =\log\pi_{0,x_0}^{(m)} +\sum_{i,j}N_{ij}\log P_{ij}^{(m)} +$$ + +The log-likelihood ratio is + +$$ +\log \frac{L_T^{(f)}}{L_T^{(g)}} = \log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} + \sum_{i,j}N_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} +$$ (eq:llr_markov) + +### KL divergence rate + +By the ergodic theorem for irreducible, aperiodic Markov chains, we have + +$$ +\frac{N_{ij}}{T} \xrightarrow{a.s.} \pi_i^{(f)}P_{ij}^{(f)} \quad \text{as } T \to \infty +$$ + +where $\boldsymbol{\pi}^{(f)}$ is the stationary distribution satisfying $\boldsymbol{\pi}^{(f)} = \boldsymbol{\pi}^{(f)} P^{(f)}$. + +Therefore, + +$$ +\frac{1}{T}\log \frac{L_T^{(f)}}{L_T^{(g)}} = \frac{1}{T}\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} + \frac{1}{T}\sum_{i,j}N_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} +$$ + +Taking the limit as $T \to \infty$, we have: +- The first term: $\frac{1}{T}\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} \to 0$ +- The second term: $\frac{1}{T}\sum_{i,j}N_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} \xrightarrow{a.s.} \sum_{i,j}\pi_i^{(f)}P_{ij}^{(f)}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}}$ + +Define the **KL divergence rate** as + +$$ +h_{KL}(f, g) = \sum_{i=1}^n \pi_i^{(f)} \underbrace{\sum_{j=1}^n P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}}}_{=: KL(P_{i\cdot}^{(f)}, P_{i\cdot}^{(g)})} +$$ + +where $KL(P_{i\cdot}^{(f)}, P_{i\cdot}^{(g)})$ is the row-wise KL divergence. + + +By the ergodic theorem, we have + +$$ +\frac{1}{T}\log \frac{L_T^{(f)}}{L_T^{(g)}} \xrightarrow{a.s.} h_{KL}(f, g) \quad \text{as } T \to \infty +$$ + +Taking expectations and using the dominated convergence theorem, we obtain + +$$ +\frac{1}{T}E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] \to h_{KL}(f, g) \quad \text{as } T \to \infty +$$ + +Here we invite readers to pause and compare this result with {eq}`eq:kl_likelihood_link`. + +Let's confirm this in the simulation below. + +### Simulations + +Let's implement simulations to illustrate these concepts with a three-state Markov chain. + +We start by writing functions to compute the stationary distribution and the KL divergence rate for Markov chain models. + +```{code-cell} ipython3 +:tags: [hide-input] + +def compute_stationary_dist(P): + """ + Compute stationary distribution of transition matrix P + """ + eigenvalues, eigenvectors = np.linalg.eig(P.T) + idx = np.argmax(np.abs(eigenvalues)) + stationary = np.real(eigenvectors[:, idx]) + return stationary / stationary.sum() + +def markov_kl_divergence(P_f, P_g, pi_f): + """ + Compute KL divergence rate between two Markov chains + """ + if np.any((P_f > 0) & (P_g == 0)): + return np.inf + + valid_mask = (P_f > 0) & (P_g > 0) + log_ratios = np.zeros_like(P_f) + log_ratios[valid_mask] = np.log(P_f[valid_mask] / P_g[valid_mask]) + + # Weight by stationary probabilities and sum + kl_rate = np.sum(pi_f[:, np.newaxis] * P_f * log_ratios) + return kl_rate + +def simulate_markov_chain(P, pi_0, T, N_paths=1000): + """ + Simulate N_paths sample paths from a Markov chain + """ + mc = qe.MarkovChain(P, state_values=None) + initial_states = np.random.choice(len(P), size=N_paths, p=pi_0) + paths = np.zeros((N_paths, T+1), dtype=int) + + for i in range(N_paths): + path = mc.simulate(T+1, init=initial_states[i]) + paths[i, :] = path + + return paths + +def compute_likelihood_ratio_markov(paths, P_f, P_g, π_0_f, π_0_g): + """ + Compute likelihood ratio process for Markov chain paths + """ + N_paths, T_plus_1 = paths.shape + T = T_plus_1 - 1 + L_ratios = np.ones((N_paths, T+1)) + + # Initial likelihood ratio + L_ratios[:, 0] = π_0_f[paths[:, 0]] / π_0_g[paths[:, 0]] + + # Compute sequential likelihood ratios + for t in range(1, T+1): + prev_states = paths[:, t-1] + curr_states = paths[:, t] + + transition_ratios = (P_f[prev_states, curr_states] / + P_g[prev_states, curr_states]) + L_ratios[:, t] = L_ratios[:, t-1] * transition_ratios + + return L_ratios + +def analyze_markov_chains(P_f, P_g, + T=500, N_paths=1000, plot_paths=True, n_show=50): + """ + Complete analysis of two Markov chains + """ + # Compute stationary distributions + π_f = compute_stationary_dist(P_f) + π_g = compute_stationary_dist(P_g) + + print(f"Stationary distribution (f): {π_f}") + print(f"Stationary distribution (g): {π_g}") + + # Compute KL divergence rates + kl_rate_fg = markov_kl_divergence(P_f, P_g, π_f) + kl_rate_gf = markov_kl_divergence(P_g, P_f, π_g) + + print(f"\nKL divergence rate h(f, g): {kl_rate_fg:.4f}") + print(f"KL divergence rate h(g, f): {kl_rate_gf:.4f}") + + if plot_paths: + # Simulate and plot paths + paths_from_f = simulate_markov_chain(P_f, π_f, T, N_paths) + L_ratios_f = compute_likelihood_ratio_markov( + paths_from_f, P_f, P_g, π_f, π_g) + + plt.figure(figsize=(10, 6)) + + # Plot individual paths + for i in range(min(n_show, N_paths)): + plt.plot(np.log(L_ratios_f[i, :]), alpha=0.3, color='blue', lw=0.8) + + # Plot theoretical expectation + theory_line = kl_rate_fg * np.arange(T+1) + plt.plot(theory_line, 'k--', linewidth=2.5, + label=r'$T \times h_{KL}(f,g)$') + + # Plot empirical mean + avg_log_L = np.mean(np.log(L_ratios_f), axis=0) + plt.plot(avg_log_L, 'r-', linewidth=2.5, + label='empirical average', alpha=0.7) + + plt.axhline(y=0, color='gray', linestyle='--', alpha=0.5) + plt.xlabel(r'$T$') + plt.ylabel(r'$\log L_T$') + plt.title('Markov chain likelihood ratios (nature = f)') + plt.legend() + plt.show() + + return { + 'stationary_f': π_f, + 'stationary_g': π_g, + 'kl_rate_fg': kl_rate_fg, + 'kl_rate_gf': kl_rate_gf + } + +def compute_markov_selection_error(T_values, P_f, P_g, π_0_f, π_0_g, N_sim=1000): + """ + Compute model selection error probability for Markov chains + """ + errors = [] + + for T in T_values: + # Simulate from both models + paths_f = simulate_markov_chain(P_f, π_0_f, T, N_sim//2) + paths_g = simulate_markov_chain(P_g, π_0_g, T, N_sim//2) + + # Compute likelihood ratios + L_f = compute_likelihood_ratio_markov(paths_f, P_f, P_g, π_0_f, π_0_g) + L_g = compute_likelihood_ratio_markov(paths_g, P_f, P_g, π_0_f, π_0_g) + + # Decision rule: choose f if L_T >= 1 + error_f = np.mean(L_f[:, -1] < 1) # Type I error + error_g = np.mean(L_g[:, -1] >= 1) # Type II error + + total_error = 0.5 * (error_f + error_g) + errors.append(total_error) + + return np.array(errors) +``` + +Now let's create an example with two different 3-state Markov chains. + +We are now ready to simulate paths and visualize how likelihood ratios evolve. + +We verify $\frac{1}{T}E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = h_{KL}(f, g)$ starting from the stationary distribution by plotting both the empirical average and the line predicted by the theory + +```{code-cell} ipython3 +# Define example Markov chain transition matrices +P_f = np.array([[0.7, 0.2, 0.1], + [0.3, 0.5, 0.2], + [0.1, 0.3, 0.6]]) + +P_g = np.array([[0.5, 0.3, 0.2], + [0.2, 0.6, 0.2], + [0.2, 0.2, 0.6]]) + +markov_results = analyze_markov_chains(P_f, P_g) +``` ## Related Lectures Likelihood processes play an important role in Bayesian learning, as described in {doc}`likelihood_bayes` and as applied in {doc}`odu`. -Likelihood ratio processes appear again in {doc}`advanced:additive_functionals`, which contains another illustration -of the **peculiar property** of likelihood ratio processes described above. +Likelihood ratio processes are central to Lawrence Blume and David Easley's answer to their question ''If you're so smart, why aren't you rich?'' {cite}`blume2006if`, the subject of the lecture{doc}`likelihood_ratio_process_2`. + +Likelihood ratio processes also appear in {doc}`advanced:additive_functionals`, which contains another illustration of the **peculiar property** of likelihood ratio processes described above. ## Exercises @@ -2496,7 +1775,7 @@ $$ Now, from the definition of Kullback-Leibler divergence $$ -K_f = KL(h, f) = \int h(w) \log \frac{h(w)}{f(w)} dw = E_h[\log h(w)] - E_h[\log f(w)] +K_f = \int h(w) \log \frac{h(w)}{f(w)} dw = E_h[\log h(w)] - E_h[\log f(w)] $$ This gives us @@ -2532,7 +1811,7 @@ Building on {ref}`lr_ex1`, use the result to explain what happens to $L_t$ as $t 1. When $K_g > K_f$ (i.e., $f$ is "closer" to $h$ than $g$ is) 2. When $K_g < K_f$ (i.e., $g$ is "closer" to $h$ than $f$ is) -Relate your answer to the simulation results shown in the {ref}`Kullback-Leibler Divergence ` section. +Relate your answer to the simulation results shown in {ref}`this section `. ``` ```{solution-start} lr_ex2 @@ -2569,94 +1848,3 @@ Therefore by similar reasoning $L_t \to 0$ almost surely. ```{solution-end} ``` - -```{exercise} -:label: lr_ex3 - -Starting from {eq}`eq:priceequation1`, show that the competitive equilibrium prices can be expressed as - -$$ -p_t(s^t) = \frac{\delta^t}{\lambda(1-\lambda)} \pi_t^2(s^t) \bigl[1 - \lambda + \lambda l_t(s^t)\bigr] -$$ - -``` - -```{solution-start} lr_ex3 -:class: dropdown -``` - -Starting from - -$$ -p_t(s^t) = \frac{\delta^t \pi_t^i(s^t)}{\mu_i c_t^i(s^t)}, \qquad i=1,2. -$$ - -Since both expressions equal the same price, we can equate them - -$$ -\frac{\pi_t^1(s^t)}{\mu_1 c_t^1(s^t)} = \frac{\pi_t^2(s^t)}{\mu_2 c_t^2(s^t)} -$$ - -Rearranging gives - -$$ -\frac{c_t^1(s^t)}{c_t^2(s^t)} = \frac{\mu_2}{\mu_1} l_t(s^t) -$$ - -where $l_t(s^t) \equiv \pi_t^1(s^t)/\pi_t^2(s^t)$ is the likelihood ratio process. - -Using $c_t^2(s^t) = 1 - c_t^1(s^t)$: - -$$ -\frac{c_t^1(s^t)}{1 - c_t^1(s^t)} = \frac{\mu_2}{\mu_1} l_t(s^t) -$$ - -Solving for $c_t^1(s^t)$ - -$$ -c_t^1(s^t) = \frac{\mu_2 l_t(s^t)}{\mu_1 + \mu_2 l_t(s^t)} -$$ - - -The planner's solution gives - -$$ -c_t^1(s^t) = \frac{\lambda l_t(s^t)}{1 - \lambda + \lambda l_t(s^t)} -$$ - -To match them, we need the following equality to hold - -$$ -\frac{\mu_2}{\mu_1} = \frac{\lambda}{1 - \lambda} -$$ - -Hence we have - -$$ -\mu_1 = 1 - \lambda, \qquad \mu_2 = \lambda -$$ - - -With $\mu_1 = 1-\lambda$ and $c_t^1(s^t) = \frac{\lambda l_t(s^t)}{1-\lambda+\lambda l_t(s^t)}$, -we have - -$$ -\begin{aligned} -p_t(s^t) &= \frac{\delta^t \pi_t^1(s^t)}{(1-\lambda) c_t^1(s^t)} \\ -&= \frac{\delta^t \pi_t^1(s^t)}{(1-\lambda)} \cdot \frac{1 - \lambda + \lambda l_t(s^t)}{\lambda l_t(s^t)} \\ -&= \frac{\delta^t \pi_t^1(s^t)}{(1-\lambda)\lambda l_t(s^t)} \bigl[1 - \lambda + \lambda l_t(s^t)\bigr]. -\end{aligned} -$$ - -Since $\pi_t^1(s^t) = l_t(s^t) \pi_t^2(s^t)$, we have - -$$ -\begin{aligned} -p_t(s^t) &= \frac{\delta^t l_t(s^t) \pi_t^2(s^t)}{(1-\lambda)\lambda l_t(s^t)} \bigl[1 - \lambda + \lambda l_t(s^t)\bigr] \\ -&= \frac{\delta^t \pi_t^2(s^t)}{(1-\lambda)\lambda} \bigl[1 - \lambda + \lambda l_t(s^t)\bigr] \\ -&= \frac{\delta^t}{\lambda(1-\lambda)} \pi_t^2(s^t) \bigl[1 - \lambda + \lambda l_t(s^t)\bigr]. -\end{aligned} -$$ - -```{solution-end} -``` diff --git a/lectures/likelihood_ratio_process_2.md b/lectures/likelihood_ratio_process_2.md new file mode 100644 index 000000000..0fad9ebd9 --- /dev/null +++ b/lectures/likelihood_ratio_process_2.md @@ -0,0 +1,930 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.17.1 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +(likelihood_ratio_process_2)= +```{raw} jupyter + +``` + +# Heterogeneous Beliefs and Financial Markets + +```{contents} Contents +:depth: 2 +``` + +(overview)= +## Overview + +A likelihood ratio process lies behind Lawrence Blume and David Easley's answer to their question +''If you're so smart, why aren't you rich?'' {cite}`blume2006if`. + +Blume and Easley constructed formal models to study how differences of opinions about probabilities governing risky income processes would influence outcomes and be reflected in prices of stocks, bonds, and insurance policies that individuals use to share and hedge risks. + +```{note} +{cite}`alchian1950uncertainty` and {cite}`friedman1953essays` conjectured that, by rewarding traders with more realistic probability models, competitive markets in financial securities put wealth in the hands of better informed traders and help +make prices of risky assets reflect realistic probability assessments. +``` + + +Here we'll provide an example that illustrates basic components of Blume and Easley's analysis. + +We'll focus only on their analysis of an environment with complete markets in which trades in all conceivable risky securities are possible. + +We'll study two alternative arrangements: + +* perfect socialism in which individuals surrender their endowments of consumption goods each period to a central planner who then dictatorially allocates those goods +* a decentralized system of competitive markets in which selfish price-taking individuals voluntarily trade with each other in competitive markets + +The fundamental theorems of welfare economics will apply and assure us that these two arrangements end up producing exactly the same allocation of consumption goods to individuals **provided** that the social planner assigns an appropriate set of **Pareto weights**. + +```{note} +You can learn about how the two welfare theorems are applied in modern macroeconomic models in {doc}`this lecture on a planning problem ` and {doc}`this lecture on a related competitive equilibrium `. +``` + + + + +Let's start by importing some Python tools. + +```{code-cell} ipython3 +import matplotlib.pyplot as plt +import numpy as np +from numba import vectorize, jit +from math import gamma +from scipy.integrate import quad +from scipy.optimize import brentq, minimize_scalar +import pandas as pd +from IPython.display import display, Math +import quantecon as qe +``` + +## Review: Likelihood Ratio Processes + +We'll begin by reminding ourselves definitions and properties of likelihood ratio processes. + +A nonnegative random variable $W$ has one of two probability density functions, either +$f$ or $g$. + +Before the beginning of time, nature once and for all decides whether she will draw a sequence of IID draws from either +$f$ or $g$. + +We will sometimes let $q$ be the density that nature chose once and for all, so +that $q$ is either $f$ or $g$, permanently. + +Nature knows which density it permanently draws from, but we the observers do not. + +We know both $f$ and $g$ but we don't know which density nature +chose. + +But we want to know. + +To do that, we use observations. + +We observe a sequence $\{w_t\}_{t=1}^T$ of $T$ IID draws that we know came from either $f$ or $g$. + +We want to use these observations to infer whether nature chose $f$ or $g$. + +A **likelihood ratio process** is a useful tool for this task. + +To begin, we define a key component of a likelihood ratio process, namely, the time $t$ likelihood ratio as the random variable + +$$ +\ell (w_t)=\frac{f\left(w_t\right)}{g\left(w_t\right)},\quad t\geq1. +$$ + +We assume that $f$ and $g$ both put positive probabilities on the +same intervals of possible realizations of the random variable $W$. + +That means that under the $g$ density, $\ell (w_t)= +\frac{f\left(w_{t}\right)}{g\left(w_{t}\right)}$ +is a nonnegative random variable with mean $1$. + +A **likelihood ratio process** for sequence +$\left\{ w_{t}\right\} _{t=1}^{\infty}$ is defined as + +$$ +L\left(w^{t}\right)=\prod_{i=1}^{t} \ell (w_i), +$$ + +where $w^t=\{ w_1,\dots,w_t\}$ is a history of +observations up to and including time $t$. + +Sometimes for shorthand we'll write $L_t = L(w^t)$. + +Notice that the likelihood process satisfies the *recursion* + +$$ +L(w^t) = \ell (w_t) L (w^{t-1}) . +$$ + +The likelihood ratio and its logarithm are key tools for making +inferences using a classic frequentist approach due to Neyman and +Pearson {cite}`Neyman_Pearson`. + +To help us appreciate how things work, the following Python code evaluates $f$ and $g$ as two different +Beta distributions, then computes and simulates an associated likelihood +ratio process by generating a sequence $w^t$ from one of the two +probability distributions, for example, a sequence of IID draws from $g$. + +```{code-cell} ipython3 +# Parameters in the two Beta distributions. +F_a, F_b = 1, 1 +G_a, G_b = 3, 1.2 + +@vectorize +def p(x, a, b): + r = gamma(a + b) / (gamma(a) * gamma(b)) + return r * x** (a-1) * (1 - x) ** (b-1) + +# The two density functions. +f = jit(lambda x: p(x, F_a, F_b)) +g = jit(lambda x: p(x, G_a, G_b)) +``` + +```{code-cell} ipython3 +@jit +def simulate(a, b, T=50, N=500): + ''' + Generate N sets of T observations of the likelihood ratio, + return as N x T matrix. + ''' + + l_arr = np.empty((N, T)) + + for i in range(N): + for j in range(T): + w = np.random.beta(a, b) + l_arr[i, j] = f(w) / g(w) + + return l_arr +``` + +## Blume and Easley's Setting + +Let the random variable $s_t \in (0,1)$ at time $t =0, 1, 2, \ldots$ be distributed according to the same Beta distribution with parameters +$\theta = \{\theta_1, \theta_2\}$. + +We'll denote this probability density as + +$$ +\pi(s_t|\theta) +$$ + +Below, we'll often just write $\pi(s_t)$ instead of $\pi(s_t|\theta)$ to save space. + +Let $s_t \equiv y_t^1$ be the endowment of a nonstorable consumption good that a person we'll call "agent 1" receives at time $t$. + +Let a history $s^t = [s_t, s_{t-1}, \ldots, s_0]$ be a sequence of i.i.d. random variables with joint distribution + +$$ +\pi_t(s^t) = \pi(s_t) \pi(s_{t-1}) \cdots \pi(s_0) +$$ + +So in our example, the history $s^t$ is a comprehensive record of agent $1$'s endowments of the consumption good from time $0$ up to time $t$. + +If agent $1$ were to live on an island by himself, agent $1$'s consumption $c^1(s_t)$ at time $t$ is + +$$c^1(s_t) = y_t^1 = s_t. $$ + +But in our model, agent 1 is not alone. + +## Nature and Agents' Beliefs + +Nature draws i.i.d. sequences $\{s_t\}_{t=0}^\infty$ from $\pi_t(s^t)$. + +* so $\pi$ without a superscript is nature's model +* but in addition to nature, there are other entities inside our model -- artificial people that we call "agents" +* each agent has a sequence of probability distributions over $s^t$ for $t=0, \ldots$ +* agent $i$ thinks that nature draws i.i.d. sequences $\{s_t\}_{t=0}^\infty$ from $\{\pi_t^i(s^t)\}_{t=0}^\infty$ + * agent $i$ is mistaken unless $\pi_t^i(s^t) = \pi_t(s^t)$ + +```{note} +A **rational expectations** model would set $\pi_t^i(s^t) = \pi_t(s^t)$ for all agents $i$. +``` + +There are two agents named $i=1$ and $i=2$. + +At time $t$, agent $1$ receives an endowment + +$$ +y_t^1 = s_t +$$ + +of a nonstorable consumption good, while agent $2$ receives an endowment of + +$$ +y_t^2 = 1 - s_t +$$ + +The aggregate endowment of the consumption good is + +$$ +y_t^1 + y_t^2 = 1 +$$ + +at each date $t \geq 0$. + +At date $t$ agent $i$ consumes $c_t^i(s^t)$ of the good. + +A (non wasteful) feasible allocation of the aggregate endowment of $1$ each period satisfies + +$$ +c_t^1 + c_t^2 = 1 . +$$ + +## A Socialist Risk-Sharing Arrangement + +In order to share risks, a benevolent social planner dictates a history-dependent consumption allocation that takes the form of a sequence of functions + +$$ +c_t^i = c_t^i(s^t) +$$ + +that satisfy + +$$ +c_t^1(s^t) + c_t^2(s^t) = 1 +$$ (eq:feasibility) + +for all $s^t$ for all $t \geq 0$. + +To design a socially optimal allocation, the social planner wants to know what agent $1$ believes about the endowment sequence and how they feel about bearing risks. + +As for the endowment sequences, agent $i$ believes that nature draws i.i.d. sequences from joint densities + +$$ +\pi_t^i(s^t) = \pi(s_t)^i \pi^i(s_{t-1}) \cdots \pi^i(s_0) +$$ + +As for attitudes toward bearing risks, agent $i$ has a one-period utility function + +$$ +u(c_t^i) = \ln (c_t^i) +$$ + +with marginal utility of consumption in period $i$ + +$$ +u'(c_t^i) = \frac{1}{c_t^i} +$$ + +Putting its beliefs about its random endowment sequence and its attitudes toward bearing risks together, agent $i$ has intertemporal utility function + +$$ +V^i = \sum_{t=0}^{\infty} \sum_{s^t} \delta^t u(c_t^i(s^t)) \pi_t^i(s^t) , +$$ (eq:objectiveagenti) + +where $\delta \in (0,1)$ is an intertemporal discount factor, and $u(\cdot)$ is a strictly increasing, concave one-period utility function. + + +## Social Planner's Allocation Problem + +The benevolent dictator has all the information it requires to choose a consumption allocation that maximizes the social welfare criterion + +$$ +W = \lambda V^1 + (1-\lambda) V^2 +$$ (eq:welfareW) + +where $\lambda \in [0,1]$ is a Pareto weight that tells how much the planner likes agent $1$ and $1 - \lambda$ is a Pareto weight that tells how much the social planner likes agent $2$. + +Setting $\lambda = .5$ expresses ''egalitarian'' social preferences. + +Notice how social welfare criterion {eq}`eq:welfareW` takes into account both agents' preferences as represented by formula {eq}`eq:objectiveagenti`. + +This means that the social planner knows and respects + +* each agent's one period utility function $u(\cdot) = \ln(\cdot)$ +* each agent $i$'s probability model $\{\pi_t^i(s^t)\}_{t=0}^\infty$ + +Consequently, we anticipate that these objects will appear in the social planner's rule for allocating the aggregate endowment each period. + + +First-order necessary conditions for maximizing welfare criterion {eq}`eq:welfareW` subject to the feasibility constraint {eq}`eq:feasibility` are + +$$\frac{\pi_t^2(s^t)}{\pi_t^1(s^t)} \frac{(1/c_t^2(s^t))}{(1/c_t^1(s^t))} = \frac{\lambda}{1-\lambda}$$ + +which can be rearranged to become + + + + +$$ +\frac{c_t^1(s^t)}{c_t^2(s^t)} = \frac{\lambda}{1-\lambda} l_t(s^t) +$$ (eq:allocationrule0) + + +where + +$$ l_t(s^t) = \frac{\pi_t^1(s^t)}{\pi_t^2(s^t)} $$ + +is the likelihood ratio of agent 1's joint density to agent 2's joint density. + +Using + +$$c_t^1(s^t) + c_t^2(s^t) = 1$$ + +we can rewrite allocation rule {eq}`eq:allocationrule0` as + + + +$$\frac{c_t^1(s^t)}{1 - c_t^1(s^t)} = \frac{\lambda}{1-\lambda} l_t(s^t)$$ + +or + +$$c_t^1(s^t) = \frac{\lambda}{1-\lambda} l_t(s^t)(1 - c_t^1(s^t))$$ + +which implies that the social planner's allocation rule is + +$$ +c_t^1(s^t) = \frac{\lambda l_t(s^t)}{1-\lambda + \lambda l_t(s^t)} +$$ (eq:allocationrule1) + +If we define a temporary or **continuation Pareto weight** process as + +$$ +\lambda_t(s^t) = \frac{\lambda l_t(s^t)}{1-\lambda + \lambda l_t(s^t)}, +$$ + +then we can represent the social planner's allocation rule as + +$$ +c_t^1(s^t) = \lambda_t(s^t) . +$$ + + + + +## If You're So Smart, $\ldots$ + + +Let's compute some values of limiting allocations {eq}`eq:allocationrule1` for some interesting possible limiting +values of the likelihood ratio process $l_t(s^t)$: + + $$l_\infty (s^\infty)= 1; \quad c_\infty^1 = \lambda$$ + + * In the above case, both agents are equally smart (or equally not smart) and the consumption allocation stays put at a $\lambda, 1 - \lambda$ split between the two agents. + +$$l_\infty (s^\infty) = 0; \quad c_\infty^1 = 0$$ + +* In the above case, agent 2 is ''smarter'' than agent 1, and agent 1's share of the aggregate endowment converges to zero. + + + + +$$l_\infty (s^\infty)= \infty; \quad c_\infty^1 = 1$$ + +* In the above case, agent 1 is smarter than agent 2, and agent 1's share of the aggregate endowment converges to 1. + +```{note} +These three cases are somehow telling us about how relative **wealths** of the agents evolve as time passes. +* when the two agents are equally smart and $\lambda \in (0,1)$, agent 1's wealth share stays at $\lambda$ perpetually. +* when agent 1 is smarter and $\lambda \in (0,1)$, agent 1 eventually "owns" the entire continuation endowment and agent 2 eventually "owns" nothing. +* when agent 2 is smarter and $\lambda \in (0,1)$, agent 2 eventually "owns" the entire continuation endowment and agent 1 eventually "owns" nothing. +Continuation wealths can be defined precisely after we introduce a competitive equilibrium **price** system below. +``` + + +Soon we'll do some simulations that will shed further light on possible outcomes. + +But before we do that, let's take a detour and study some "shadow prices" for the social planning problem that can readily be +converted to "equilibrium prices" for a competitive equilibrium. + +Doing this will allow us to connect our analysis with an argument of {cite}`alchian1950uncertainty` and {cite}`friedman1953essays` that competitive market processes can make prices of risky assets better reflect realistic probability assessments. + + + +## Competitive Equilibrium Prices + +Two fundamental welfare theorems for general equilibrium models lead us to anticipate that there is a connection between the allocation that solves the social planning problem we have been studying and the allocation in a **competitive equilibrium** with complete markets in history-contingent commodities. + +```{note} +For the two welfare theorems and their history, see . +Again, for applications to a classic macroeconomic growth model, see {doc}`this lecture on a planning problem ` and {doc}`this lecture on a related competitive equilibrium ` +``` + +Such a connection prevails for our model. + +We'll sketch it now. + +In a competitive equilibrium, there is no social planner that dictatorially collects everybody's endowments and then reallocates them. + +Instead, there is a comprehensive centralized market that meets at one point in time. + +There are **prices** at which price-taking agents can buy or sell whatever goods that they want. + +Trade is multilateral in the sense that that there is a "Walrasian auctioneer" who lives outside the model and whose job is to verify that +each agent's budget constraint is satisfied. + +That budget constraint involves the total value of the agent's endowment stream and the total value of its consumption stream. + +These values are computed at price vectors that the agents take as given -- they are "price-takers" who assume that they can buy or sell +whatever quantities that they want at those prices. + +Suppose that at time $-1$, before time $0$ starts, agent $i$ can purchase one unit $c_t(s^t)$ of consumption at time $t$ after history +$s^t$ at price $p_t(s^t)$. + +Notice that there is (very long) **vector** of prices. + + * there is one price $p_t(s^t)$ for each history $s^t$ at every date $t = 0, 1, \ldots, $. + * so there are as many prices as there are histories and dates. + +These prices determined at time $-1$ before the economy starts. + +The market meets once at time $-1$. + +At times $t =0, 1, 2, \ldots$ trades made at time $-1$ are executed. + + + +* in the background, there is an "enforcement" procedure that forces agents to carry out the exchanges or "deliveries" that they agreed to at time $-1$. + + + +We want to study how agents' beliefs influence equilibrium prices. + +Agent $i$ faces a **single** intertemporal budget constraint + +$$ +\sum_{t=0}^\infty\sum_{s^t} p_t(s^t) c_t^i (s^t) \leq \sum_{t=0}^\infty\sum_{s^t} p_t(s^t) y_t^i (s^t) +$$ (eq:budgetI) + +According to budget constraint {eq}`eq:budgetI`, trade is **multilateral** in the following sense + +* we can imagine that agent $i$ first sells his random endowment stream $\{y_t^i (s^t)\}$ and then uses the proceeds (i.e., his "wealth") to purchase a random consumption stream $\{c_t^i (s^t)\}$. + +Agent $i$ puts a Lagrange multiplier $\mu_i$ on {eq}`eq:budgetI` and once-and-for-all chooses a consumption plan $\{c^i_t(s^t)\}_{t=0}^\infty$ +to maximize criterion {eq}`eq:objectiveagenti` subject to budget constraint {eq}`eq:budgetI`. + +This means that the agent $i$ chooses many objects, namely, $c_t^i(s^t)$ for all $s^t$ for $t = 0, 1, 2, \ldots$. + + +For convenience, let's remind ourselves of criterion $V^i$ defined in {eq}`eq:objectiveagenti`: + +$$ +V^i = \sum_{t=0}^{\infty} \sum_{s^t} \delta^t u(c_t^i(s^t)) \pi_t^i(s^t) +$$ + +First-order necessary conditions for maximizing objective $V^i$ defined in {eq}`eq:objectiveagenti` with respect to $c_t^i(s^t)$ are + +$$ +\delta^t u'(c^i_t(s^t)) \pi_t^i(s^t) = \mu_i p_t(s^t) , +$$ + +which we can rearrange to obtain + +$$ +p_t(s^t) = \frac{ \delta^t \pi_t^i(s^t)}{\mu_i c^i_t(s^t)} +$$ (eq:priceequation1) + +for $i=1,2$. + +If we divide equation {eq}`eq:priceequation1` for agent $1$ by the appropriate version of equation {eq}`eq:priceequation1` for agent 2, use +$c^2_t(s^t) = 1 - c^1_t(s^t)$, and do some algebra, we'll obtain + +$$ +c_t^1(s^t) = \frac{\mu_1 l_t(s^t)}{\mu_2 + \mu_1 l_t(s^t)} . +$$ (eq:allocationce) + +We now engage in an extended "guess-and-verify" exercise that involves matching objects in our competitive equilibrium with objects in +our social planning problem. + +* we'll match consumption allocations in the planning problem with equilibrium consumption allocations in the competitive equilibrium +* we'll match "shadow" prices in the planning problem with competitive equilibrium prices. + +Notice that if we set $\mu_1 = 1-\lambda$ and $\mu_2 = \lambda$, then formula {eq}`eq:allocationce` agrees with formula +{eq}`eq:allocationrule1`. + + * doing this amounts to choosing a **numeraire** or normalization for the price system $\{p_t(s^t)\}_{t=0}^\infty$ + +```{note} +For information about how a numeraire must be chosen to pin down the absolute price level in a model like ours that determines only +relative prices, see . +``` + +If we substitute formula {eq}`eq:allocationce` for $c_t^1(s^t)$ into formula {eq}`eq:priceequation1` and rearrange, we obtain + +$$ +p_t(s^t) = \frac{\delta^t}{\lambda(1-\lambda)} \pi_t^2(s^t) \bigl[1 - \lambda + \lambda l_t(s^t)\bigr] +$$ + +or + +$$ +p_t(s^t) = \frac{\delta^t}{\lambda(1-\lambda)} \bigl[(1 - \lambda) \pi_t^2(s^t) + \lambda \pi_t^1(s^t)\bigr] +$$ (eq:pformulafinal) + +According to formula {eq}`eq:pformulafinal`, we have the following possible limiting cases: + +* when $l_\infty = 0$, $c_\infty^1 = 0 $ and tails of competitive equilibrium prices reflect agent $2$'s probability model $\pi_t^2(s^t)$ according to $p_t(s^t) \propto \delta^t \pi_t^2(s^t) $ +* when $l_\infty = \infty$, $c_\infty^1 = 1 $ and tails of competitive equilibrium prices reflect agent $1$'s probability model $\pi_t^1(s^t)$ according to $p_t(s^t) \propto \delta^t \pi_t^1(s^t) $ +* for small $t$'s, competitive equilibrium prices reflect both agents' probability models. + +## Simulations + +Now let's implement some simulations when agent $1$ believes marginal density + +$$\pi^1(s_t) = f(s_t) $$ + +and agent $2$ believes marginal density + +$$ \pi^2(s_t) = g(s_t) $$ + +where $f$ and $g$ are Beta distributions like ones that we used in earlier sections of this lecture. + +Meanwhile, we'll assume that nature believes a marginal density + +$$ +\pi(s_t) = h(s_t) +$$ + +where $h(s_t)$ is perhaps a mixture of $f$ and $g$. + + +First, we write a function to compute the likelihood ratio process + +```{code-cell} ipython3 +def compute_likelihood_ratios(sequences, f, g): + """Compute likelihood ratios and cumulative products.""" + l_ratios = f(sequences) / g(sequences) + L_cumulative = np.cumprod(l_ratios, axis=1) + return l_ratios, L_cumulative +``` + +Let's compute the Kullback–Leibler discrepancies by quadrature +integration. + +```{code-cell} ipython3 +def compute_KL(f, g): + """ + Compute KL divergence KL(f, g) + """ + integrand = lambda w: f(w) * np.log(f(w) / g(w)) + val, _ = quad(integrand, 1e-5, 1-1e-5) + return val +``` + +We also create a helper function to compute KL divergence with respect to a reference distribution $h$ + +```{code-cell} ipython3 +def compute_KL_h(h, f, g): + """ + Compute KL divergence with reference distribution h + """ + + Kf = compute_KL(h, f) + Kg = compute_KL(h, g) + + return Kf, Kg +``` + +Let's write a Python function that computes agent 1's consumption share + +```{code-cell} ipython3 +def simulate_blume_easley(sequences, f_belief=f, g_belief=g, λ=0.5): + """Simulate Blume-Easley model consumption shares.""" + l_ratios, l_cumulative = compute_likelihood_ratios(sequences, f_belief, g_belief) + c1_share = λ * l_cumulative / (1 - λ + λ * l_cumulative) + return l_cumulative, c1_share +``` + +Now let's use this function to generate sequences in which + +* nature draws from $f$ each period, or +* nature draws from $g$ each period, or +* nature flips a fair coin each period to decide whether to draw from $f$ or $g$ + +```{code-cell} ipython3 +λ = 0.5 +T = 100 +N = 10000 + +# Nature follows f, g, or mixture +s_seq_f = np.random.beta(F_a, F_b, (N, T)) +s_seq_g = np.random.beta(G_a, G_b, (N, T)) + +h = jit(lambda x: 0.5 * f(x) + 0.5 * g(x)) +model_choices = np.random.rand(N, T) < 0.5 +s_seq_h = np.empty((N, T)) +s_seq_h[model_choices] = np.random.beta(F_a, F_b, size=model_choices.sum()) +s_seq_h[~model_choices] = np.random.beta(G_a, G_b, size=(~model_choices).sum()) + +l_cum_f, c1_f = simulate_blume_easley(s_seq_f) +l_cum_g, c1_g = simulate_blume_easley(s_seq_g) +l_cum_h, c1_h = simulate_blume_easley(s_seq_h) +``` + +Before looking at the figure below, have some fun by guessing whether agent 1 or agent 2 will have a larger and larger consumption share as time passes in our three cases. + +To make better guesses, let's visualize instances of the likelihood ratio processes in the three cases. + +```{code-cell} ipython3 +fig, axes = plt.subplots(2, 3, figsize=(15, 10)) + +titles = ["Nature = f", "Nature = g", "Nature = mixture"] +data_pairs = [(l_cum_f, c1_f), (l_cum_g, c1_g), (l_cum_h, c1_h)] + +for i, ((l_cum, c1), title) in enumerate(zip(data_pairs, titles)): + # Likelihood ratios + ax = axes[0, i] + for j in range(min(50, l_cum.shape[0])): + ax.plot(l_cum[j, :], alpha=0.3, color='blue') + ax.set_yscale('log') + ax.set_xlabel('time') + ax.set_ylabel('Likelihood ratio $l_t$') + ax.set_title(title) + ax.axhline(y=1, color='red', linestyle='--', alpha=0.5) + + # Consumption shares + ax = axes[1, i] + for j in range(min(50, c1.shape[0])): + ax.plot(c1[j, :], alpha=0.3, color='green') + ax.set_xlabel('time') + ax.set_ylabel("Agent 1's consumption share") + ax.set_ylim([0, 1]) + ax.axhline(y=λ, color='red', linestyle='--', alpha=0.5) + +plt.tight_layout() +plt.show() +``` + +In the left panel, nature chooses $f$. Agent 1's consumption reaches $1$ very quickly. + +In the middle panel, nature chooses $g$. Agent 1's consumption ratio tends to move towards $0$ but not as fast as in the first case. + +In the right panel, nature flips coins each period. We see a very similar pattern to the processes in the left panel. + +The figures in the top panel remind us of the discussion in [this section](KL_link). + +We invite readers to revisit [that section](llr_h) and try to infer the relationships among $D_{KL}(f\|g)$, $D_{KL}(g\|f)$, $D_{KL}(h\|f)$, and $D_{KL}(h\|g)$. + + +Let's compute values of KL divergence + +```{code-cell} ipython3 +shares = [np.mean(c1_f[:, -1]), np.mean(c1_g[:, -1]), np.mean(c1_h[:, -1])] +Kf_g, Kg_f = compute_KL(f, g), compute_KL(g, f) +Kf_h, Kg_h = compute_KL_h(h, f, g) + +print(f"Final shares: f={shares[0]:.3f}, g={shares[1]:.3f}, mix={shares[2]:.3f}") +print(f"KL divergences: \nKL(f,g)={Kf_g:.3f}, KL(g,f)={Kg_f:.3f}") +print(f"KL(h,f)={Kf_h:.3f}, KL(h,g)={Kg_h:.3f}") +``` + +We find that $KL(f,g) > KL(g,f)$ and $KL(h,g) > KL(h,f)$. + +The first inequality tells us that the average "surprise" from having belief $g$ when nature chooses $f$ is greater than the "surprise" from having belief $f$ when nature chooses $g$. + +This explains the difference between the first two panels we noted above. + +The second inequality tells us that agent 1's belief distribution $f$ is closer to nature's pick than agent 2's belief $g$. + ++++ + +To make this idea more concrete, let's compare two cases: + +- agent 1's belief distribution $f$ is close to agent 2's belief distribution $g$; +- agent 1's belief distribution $f$ is far from agent 2's belief distribution $g$. + + +We use the two distributions visualized below + +```{code-cell} ipython3 +def plot_distribution_overlap(ax, x_range, f_vals, g_vals, + f_label='f', g_label='g', + f_color='blue', g_color='red'): + """Plot two distributions with their overlap region.""" + ax.plot(x_range, f_vals, color=f_color, linewidth=2, label=f_label) + ax.plot(x_range, g_vals, color=g_color, linewidth=2, label=g_label) + + overlap = np.minimum(f_vals, g_vals) + ax.fill_between(x_range, 0, overlap, alpha=0.3, color='purple', label='Overlap') + ax.set_xlabel('x') + ax.set_ylabel('Density') + ax.legend() + +# Define close and far belief distributions +f_close = jit(lambda x: p(x, 1, 1)) +g_close = jit(lambda x: p(x, 1.1, 1.05)) + +f_far = jit(lambda x: p(x, 1, 1)) +g_far = jit(lambda x: p(x, 3, 1.2)) + +# Visualize the belief distributions +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) + +x_range = np.linspace(0.001, 0.999, 200) + +# Close beliefs +f_close_vals = [f_close(x) for x in x_range] +g_close_vals = [g_close(x) for x in x_range] +plot_distribution_overlap(ax1, x_range, f_close_vals, g_close_vals, + f_label='f (Beta(1, 1))', g_label='g (Beta(1.1, 1.05))') +ax1.set_title(f'Close Beliefs') + +# Far beliefs +f_far_vals = [f_far(x) for x in x_range] +g_far_vals = [g_far(x) for x in x_range] +plot_distribution_overlap(ax2, x_range, f_far_vals, g_far_vals, + f_label='f (Beta(1, 1))', g_label='g (Beta(3, 1.2))') +ax2.set_title(f'Far Beliefs') + +plt.tight_layout() +plt.show() +``` + +Let's draw the same consumption ratio plots as above for agent 1. + +We replace the simulation paths with median and percentiles to make the figure cleaner. + +Staring at the figure below, can we infer the relation between $KL(f,g)$ and $KL(g,f)$? + +From the right panel, can we infer the relation between $KL(h,g)$ and $KL(h,f)$? + +```{code-cell} ipython3 +fig, axes = plt.subplots(2, 3, figsize=(15, 10)) +nature_params = {'close': [(1, 1), (1.1, 1.05), (2, 1.5)], + 'far': [(1, 1), (3, 1.2), (2, 1.5)]} +nature_labels = ["Nature = f", "Nature = g", "Nature = h"] +colors = {'close': 'blue', 'far': 'red'} + +threshold = 1e-5 # "close to zero" cutoff + +for row, (f_belief, g_belief, label) in enumerate([ + (f_close, g_close, 'close'), + (f_far, g_far, 'far')]): + + for col, nature_label in enumerate(nature_labels): + params = nature_params[label][col] + s_seq = np.random.beta(params[0], params[1], (1000, 200)) + _, c1 = simulate_blume_easley(s_seq, f_belief, g_belief, λ) + + median_c1 = np.median(c1, axis=0) + p10, p90 = np.percentile(c1, [10, 90], axis=0) + + ax = axes[row, col] + color = colors[label] + ax.plot(median_c1, color=color, linewidth=2, label='Median') + ax.fill_between(range(len(median_c1)), p10, p90, alpha=0.3, color=color, label='10–90%') + ax.set_xlabel('time') + ax.set_ylabel("Agent 1's share") + ax.set_ylim([0, 1]) + ax.set_title(nature_label) + ax.axhline(y=λ, color='gray', linestyle='--', alpha=0.5) + below = np.where(median_c1 < threshold)[0] + above = np.where(median_c1 > 1-threshold)[0] + if below.size > 0: first_zero = (below[0], True) + elif above.size > 0: first_zero = (above[0], False) + else: first_zero = None + if first_zero is not None: + ax.axvline(x=first_zero[0], color='black', linestyle='--', + alpha=0.7, + label=fr'Median $\leq$ {threshold}' if first_zero[1] + else fr'Median $\geq$ 1-{threshold}') + ax.legend() + +plt.tight_layout() +plt.show() +``` + +Holding to our guesses, let's calculate the four values + +```{code-cell} ipython3 +# Close case +Kf_g, Kg_f = compute_KL(f_close, g_close), compute_KL(g_close, f_close) +Kf_h, Kg_h = compute_KL_h(h, f_close, g_close) + +print(f"KL divergences (close): \nKL(f,g)={Kf_g:.3f}, KL(g,f)={Kg_f:.3f}") +print(f"KL(h,f)={Kf_h:.3f}, KL(h,g)={Kg_h:.3f}") + +# Far case +Kf_g, Kg_f = compute_KL(f_far, g_far), compute_KL(g_far, f_far) +Kf_h, Kg_h = compute_KL_h(h, f_far, g_far) + +print(f"KL divergences (far): \nKL(f,g)={Kf_g:.3f}, KL(g,f)={Kg_f:.3f}") +print(f"KL(h,f)={Kf_h:.3f}, KL(h,g)={Kg_h:.3f}") +``` + +We find that in the first case, $KL(f,g) \approx KL(g,f)$ and both are relatively small, so although either agent 1 or agent 2 will eventually consume everything, convergence displayed in the first two panels on the top is pretty slow. + +In the first two panels at the bottom, we see convergence occurring faster (as indicated by the black dashed line) because the divergence gaps $KL(f, g)$ and $KL(g, f)$ are larger. + +Since $KL(f,g) > KL(g,f)$, we see faster convergence in the first panel at the bottom when nature chooses $f$ than in the second panel where nature chooses $g$. + +This ties in nicely with {eq}`eq:kl_likelihood_link`. + + + +## Related Lectures + +Likelihood processes play an important role in Bayesian learning, as described in {doc}`likelihood_bayes` +and as applied in {doc}`odu`. + +Likelihood ratio processes appear again in {doc}`advanced:additive_functionals`. + + +## Exercise + +```{exercise} +:label: lr_ex3 + +Starting from {eq}`eq:priceequation1`, show that the competitive equilibrium prices can be expressed as + +$$ +p_t(s^t) = \frac{\delta^t}{\lambda(1-\lambda)} \pi_t^2(s^t) \bigl[1 - \lambda + \lambda l_t(s^t)\bigr] +$$ + +``` + +```{solution-start} lr_ex3 +:class: dropdown +``` + +Starting from + +$$ +p_t(s^t) = \frac{\delta^t \pi_t^i(s^t)}{\mu_i c_t^i(s^t)}, \qquad i=1,2. +$$ + +Since both expressions equal the same price, we can equate them + +$$ +\frac{\pi_t^1(s^t)}{\mu_1 c_t^1(s^t)} = \frac{\pi_t^2(s^t)}{\mu_2 c_t^2(s^t)} +$$ + +Rearranging gives + +$$ +\frac{c_t^1(s^t)}{c_t^2(s^t)} = \frac{\mu_2}{\mu_1} l_t(s^t) +$$ + +where $l_t(s^t) \equiv \pi_t^1(s^t)/\pi_t^2(s^t)$ is the likelihood ratio process. + +Using $c_t^2(s^t) = 1 - c_t^1(s^t)$: + +$$ +\frac{c_t^1(s^t)}{1 - c_t^1(s^t)} = \frac{\mu_2}{\mu_1} l_t(s^t) +$$ + +Solving for $c_t^1(s^t)$ + +$$ +c_t^1(s^t) = \frac{\mu_2 l_t(s^t)}{\mu_1 + \mu_2 l_t(s^t)} +$$ + + +The planner's solution gives + +$$ +c_t^1(s^t) = \frac{\lambda l_t(s^t)}{1 - \lambda + \lambda l_t(s^t)} +$$ + +To match them, we need the following equality to hold + +$$ +\frac{\mu_2}{\mu_1} = \frac{\lambda}{1 - \lambda} +$$ + +Hence we have + +$$ +\mu_1 = 1 - \lambda, \qquad \mu_2 = \lambda +$$ + + +With $\mu_1 = 1-\lambda$ and $c_t^1(s^t) = \frac{\lambda l_t(s^t)}{1-\lambda+\lambda l_t(s^t)}$, +we have + +$$ +\begin{aligned} +p_t(s^t) &= \frac{\delta^t \pi_t^1(s^t)}{(1-\lambda) c_t^1(s^t)} \\ +&= \frac{\delta^t \pi_t^1(s^t)}{(1-\lambda)} \cdot \frac{1 - \lambda + \lambda l_t(s^t)}{\lambda l_t(s^t)} \\ +&= \frac{\delta^t \pi_t^1(s^t)}{(1-\lambda)\lambda l_t(s^t)} \bigl[1 - \lambda + \lambda l_t(s^t)\bigr]. +\end{aligned} +$$ + +Since $\pi_t^1(s^t) = l_t(s^t) \pi_t^2(s^t)$, we have + +$$ +\begin{aligned} +p_t(s^t) &= \frac{\delta^t l_t(s^t) \pi_t^2(s^t)}{(1-\lambda)\lambda l_t(s^t)} \bigl[1 - \lambda + \lambda l_t(s^t)\bigr] \\ +&= \frac{\delta^t \pi_t^2(s^t)}{(1-\lambda)\lambda} \bigl[1 - \lambda + \lambda l_t(s^t)\bigr] \\ +&= \frac{\delta^t}{\lambda(1-\lambda)} \pi_t^2(s^t) \bigl[1 - \lambda + \lambda l_t(s^t)\bigr]. +\end{aligned} +$$ + +```{solution-end} +``` diff --git a/lectures/likelihood_var.md b/lectures/likelihood_var.md new file mode 100644 index 000000000..2bacf4812 --- /dev/null +++ b/lectures/likelihood_var.md @@ -0,0 +1,766 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.17.1 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +(var_likelihood)= +```{raw} jupyter + +``` + +# Likelihood Processes for VAR Models + +```{contents} Contents +:depth: 2 +``` + +## Overview + +This lecture extends our analysis of likelihood ratio processes to Vector Autoregressions (VARs). + +We'll + +* Construct likelihood functions for VAR models +* Form likelihood ratio processes for comparing two VAR models +* Visualize the evolution of likelihood ratios over time +* Connect VAR likelihood ratios to the Samuelson multiplier-accelerator model + +Our analysis builds on concepts from: +- {doc}`likelihood_ratio_process` +- {doc}`linear_models` +- {doc}`samuelson` + +Let's start by importing helpful libraries: + +```{code-cell} ipython3 +import numpy as np +import matplotlib.pyplot as plt +from scipy import linalg +from scipy.stats import multivariate_normal as mvn +from quantecon import LinearStateSpace +import quantecon as qe +from numba import jit +from typing import NamedTuple, Optional, Tuple +from collections import namedtuple +``` + +## VAR model setup + +Consider a VAR model of the form: + +$$ +\begin{aligned} +x_{t+1} & = A x_t + C w_{t+1} \\ +x_0 & \sim \mathcal{N}(\mu_0, \Sigma_0) +\end{aligned} +$$ + +where: +- $x_t$ is an $n \times 1$ state vector +- $w_{t+1} \sim \mathcal{N}(0, I)$ is an $m \times 1$ vector of shocks +- $A$ is an $n \times n$ transition matrix +- $C$ is an $n \times m$ volatility matrix + +Let's define the necessary data structures for the VAR model + +```{code-cell} ipython3 +VARModel = namedtuple('VARModel', ['A', 'C', 'μ_0', 'Σ_0', + 'CC', 'CC_inv', 'log_det_CC', + 'Σ_0_inv', 'log_det_Σ_0']) +def compute_stationary_var(A, C): + """ + Compute stationary mean and covariance for VAR model + """ + n = A.shape[0] + + # Check stability + eigenvalues = np.linalg.eigvals(A) + if np.max(np.abs(eigenvalues)) >= 1: + raise ValueError("VAR is not stationary") + + μ_0 = np.zeros(n) + + # Stationary covariance: solve discrete Lyapunov equation + # Σ_0 = A @ Σ_0 @ A.T + C @ C.T + CC = C @ C.T + Σ_0 = linalg.solve_discrete_lyapunov(A, CC) + + return μ_0, Σ_0 + +def create_var_model(A, C, μ_0=None, Σ_0=None, stationary=True): + """ + Create a VAR model with parameters and precomputed matrices + """ + A = np.asarray(A) + C = np.asarray(C) + n = A.shape[0] + CC = C @ C.T + + if stationary: + μ_0_comp, Σ_0_comp = compute_stationary_var(A, C) + else: + μ_0_comp = μ_0 if μ_0 is not None else np.zeros(n) + Σ_0_comp = Σ_0 if Σ_0 is not None else np.eye(n) + + # Check if CC is singular + det_CC = np.linalg.det(CC) + if np.abs(det_CC) < 1e-10: + # Use pseudo-inverse for singular case + CC_inv = np.linalg.pinv(CC) + CC_reg = CC + 1e-10 * np.eye(CC.shape[0]) + log_det_CC = np.log(np.linalg.det(CC_reg)) + else: + CC_inv = np.linalg.inv(CC) + log_det_CC = np.log(det_CC) + + # Same check for Σ_0 + det_Σ_0 = np.linalg.det(Σ_0_comp) + if np.abs(det_Σ_0) < 1e-10: + Σ_0_inv = np.linalg.pinv(Σ_0_comp) + Σ_0_reg = Σ_0_comp + 1e-10 * np.eye(Σ_0_comp.shape[0]) + log_det_Σ_0 = np.log(np.linalg.det(Σ_0_reg)) + else: + Σ_0_inv = np.linalg.inv(Σ_0_comp) + log_det_Σ_0 = np.log(det_Σ_0) + + return VARModel(A=A, C=C, μ_0=μ_0_comp, Σ_0=Σ_0_comp, + CC=CC, CC_inv=CC_inv, log_det_CC=log_det_CC, + Σ_0_inv=Σ_0_inv, log_det_Σ_0=log_det_Σ_0) +``` + +### Joint distribution + +The joint probability distribution $f(x_T, x_{T-1}, \ldots, x_0)$ can be factored as: + +$$ +f(x_T, \ldots, x_0) = f(x_T | x_{T-1}) f(x_{T-1} | x_{T-2}) \cdots f(x_1 | x_0) f(x_0) +$$ + +Since the VAR is Markovian, $f(x_{t+1} | x_t, \ldots, x_0) = f(x_{t+1} | x_t)$. + +### Conditional densities + +Given the Gaussian structure, the conditional distribution $f(x_{t+1} | x_t)$ is Gaussian with: +- Mean: $A x_t$ +- Covariance: $CC'$ + +The log conditional density is: + +$$ +\log f(x_{t+1} | x_t) = -\frac{n}{2} \log(2\pi) - \frac{1}{2} \log \det(CC') - \frac{1}{2} (x_{t+1} - A x_t)' (CC')^{-1} (x_{t+1} - A x_t) +$$ + +```{code-cell} ipython3 +def log_likelihood_transition(x_next, x_curr, model): + """ + Compute log likelihood of transition from x_curr to x_next + """ + x_next = np.atleast_1d(x_next) + x_curr = np.atleast_1d(x_curr) + n = len(x_next) + diff = x_next - model.A @ x_curr + return -0.5 * (n * np.log(2 * np.pi) + model.log_det_CC + + diff @ model.CC_inv @ diff) +``` + +The log density of the initial state is: + +$$ +\log f(x_0) = -\frac{n}{2} \log(2\pi) - \frac{1}{2} \log \det(\Sigma_0) - \frac{1}{2} (x_0 - \mu_0)' \Sigma_0^{-1} (x_0 - \mu_0) +$$ + +```{code-cell} ipython3 +def log_likelihood_initial(x_0, model): + """ + Compute log likelihood of initial state + """ + x_0 = np.atleast_1d(x_0) + n = len(x_0) + diff = x_0 - model.μ_0 + return -0.5 * (n * np.log(2 * np.pi) + model.log_det_Σ_0 + + diff @ model.Σ_0_inv @ diff) +``` + +Now let's group the likelihood computations into a single function that computes the log likelihood of an entire path + +```{code-cell} ipython3 +def log_likelihood_path(X, model): + """ + Compute log likelihood of entire path + """ + + T = X.shape[0] - 1 + log_L = log_likelihood_initial(X[0], model) + + for t in range(T): + log_L += log_likelihood_transition(X[t+1], X[t], model) + + return log_L + +def simulate_var(model, T, N_paths=1): + """ + Simulate paths from the VAR model + """ + n = model.A.shape[0] + m = model.C.shape[1] + paths = np.zeros((N_paths, T+1, n)) + + for i in range(N_paths): + # Draw initial state + x = mvn.rvs(mean=model.μ_0, cov=model.Σ_0) + x = np.atleast_1d(x) + paths[i, 0] = x + + # Simulate forward + for t in range(T): + w = np.random.randn(m) + x = model.A @ x + model.C @ w + paths[i, t+1] = x + + return paths if N_paths > 1 else paths[0] +``` + +## Likelihood ratio process + +Now let's compute likelihood ratio processes for comparing two VAR models + +```{code-cell} ipython3 +def compute_likelihood_ratio_var(paths, model_f, model_g): + """ + Compute likelihood ratio process for VAR models + """ + if paths.ndim == 2: + paths = paths[np.newaxis, :] + + N_paths, T_plus_1, n = paths.shape + T = T_plus_1 - 1 + log_L_ratios = np.zeros((N_paths, T+1)) + + for i in range(N_paths): + X = paths[i] + + # Initial log likelihood ratio + log_L_f_0 = log_likelihood_initial(X[0], model_f) + log_L_g_0 = log_likelihood_initial(X[0], model_g) + log_L_ratios[i, 0] = log_L_f_0 - log_L_g_0 + + # Recursive computation + for t in range(1, T+1): + log_L_f_t = log_likelihood_transition(X[t], X[t-1], model_f) + log_L_g_t = log_likelihood_transition(X[t], X[t-1], model_g) + + # Update log likelihood ratio + log_diff = log_L_f_t - log_L_g_t + + log_L_prev = log_L_ratios[i, t-1] + log_L_new = log_L_prev + log_diff + log_L_ratios[i, t] = log_L_new + + return log_L_ratios if N_paths > 1 else log_L_ratios[0] +``` + +## Example 1: Two AR(1) processes + +Let's start with a simple example comparing two univariate AR(1) processes with $A_f = 0.8$, $A_g = 0.5$, and $C_f = 0.3$, $C_g = 0.4$ + +```{code-cell} ipython3 +# Model f: AR(1) with persistence ρ = 0.8 +A_f = np.array([[0.8]]) +C_f = np.array([[0.3]]) + +# Model g: AR(1) with persistence ρ = 0.5 +A_g = np.array([[0.5]]) +C_g = np.array([[0.4]]) + +# Create VAR models +model_f = create_var_model(A_f, C_f) +model_g = create_var_model(A_g, C_g) +``` + +Let's generate 100 paths of length 200 from model $f$ and compute the likelihood ratio processes + +```{code-cell} ipython3 +# Simulate from model f +T = 200 +N_paths = 100 +paths_from_f = simulate_var(model_f, T, N_paths) + +L_ratios_f = compute_likelihood_ratio_var(paths_from_f, model_f, model_g) + +fig, ax = plt.subplots() + +for i in range(min(20, N_paths)): + ax.plot(L_ratios_f[i], alpha=0.3, color='C0', lw=2) + +ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5) +ax.set_ylabel(r'$\log L_t$') +ax.set_title('log likelihood ratio processes (nature = f)') + +plt.tight_layout() +plt.show() +``` + +As we expected, the likelihood ratio processes goes to $+\infty$ as $T$ increases, indicating that model $f$ is chosen correctly by our algorithm. + +## Example 2: Bivariate VAR models + +Now let's consider an example with bivariate VAR models with + +$$ +A_f & = \begin{bmatrix} 0.7 & 0.2 \\ 0.1 & 0.6 \end{bmatrix}, \quad C_f = \begin{bmatrix} 0.3 & 0.1 \\ 0.1 & 0.3 \end{bmatrix} +$$ + +and + +$$ +A_g & = \begin{bmatrix} 0.5 & 0.3 \\ 0.2 & 0.5 \end{bmatrix}, \quad C_g = \begin{bmatrix} 0.4 & 0.0 \\ 0.0 & 0.4 \end{bmatrix} +$$ + +```{code-cell} ipython3 +A_f = np.array([[0.7, 0.2], + [0.1, 0.6]]) + +C_f = np.array([[0.3, 0.1], + [0.1, 0.3]]) + +A_g = np.array([[0.5, 0.3], + [0.2, 0.5]]) + +C_g = np.array([[0.4, 0.0], + [0.0, 0.4]]) + +# Create VAR models +model2_f = create_var_model(A_f, C_f) +model2_g = create_var_model(A_g, C_g) + +# Check stationarity +print("model f eigenvalues:", np.linalg.eigvals(A_f)) +print("model g eigenvalues:", np.linalg.eigvals(A_g)) +``` + +Let's generate 50 paths of length 50 from both models and compute the likelihood ratio processes + +```{code-cell} ipython3 +# Simulate from both models +T = 50 +N_paths = 50 + +paths_from_f = simulate_var(model2_f, T, N_paths) +paths_from_g = simulate_var(model2_g, T, N_paths) + +# Compute likelihood ratios +L_ratios_ff = compute_likelihood_ratio_var(paths_from_f, model2_f, model2_g) +L_ratios_gf = compute_likelihood_ratio_var(paths_from_g, model2_f, model2_g) +``` + +We can see that for paths generated from model $f$, the likelihood ratio processes tend to go to $+\infty$, while for paths from model $g$, they tend to go to $-\infty$. + +```{code-cell} ipython3 +# Visualize the results +fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + +ax = axes[0] +for i in range(min(20, N_paths)): + ax.plot(L_ratios_ff[i], alpha=0.5, color='C0', lw=2) +ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5, lw=2) +ax.set_title(r'$\log L_t$ (nature = f)') +ax.set_ylabel(r'$\log L_t$') + +ax = axes[1] +for i in range(min(20, N_paths)): + ax.plot(L_ratios_gf[i], alpha=0.5, color='C1', lw=2) +ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5, lw=2) +ax.set_title(r'$\log L_t$ (nature = g)') +plt.tight_layout() +plt.show() +``` + +Let's apply a Neyman-Pearson frequentist decision rule described in {doc}`likelihood_ratio_process` that selects model $f$ when $\log L_T \geq 0$ and model $g$ when $\log L_T < 0$ + +```{code-cell} ipython3 +fig, ax = plt.subplots() +T_values = np.arange(0, T+1) +accuracy_f = np.zeros(len(T_values)) +accuracy_g = np.zeros(len(T_values)) + +for i, t in enumerate(T_values): + # Correct selection when data from f + accuracy_f[i] = np.mean(L_ratios_ff[:, t] > 0) + # Correct selection when data from g + accuracy_g[i] = np.mean(L_ratios_gf[:, t] < 0) + +ax.plot(T_values, accuracy_f, 'C0', linewidth=2, label='accuracy (nature = f)') +ax.plot(T_values, accuracy_g, 'C1', linewidth=2, label='accuracy (nature = g)') +ax.axhline(y=0.5, color='gray', linestyle='--', alpha=0.5) +ax.set_xlabel('T') +ax.set_ylabel('accuracy') +ax.legend() + +plt.tight_layout() +plt.show() +``` + +Evidently, the accuracy approaches $1$ as $T$ increases, and it does so very quickly. + +Let's also check the type I and type II errors as functions of $T$ + +```{code-cell} ipython3 +def model_selection_analysis(T_values, model_f, model_g, N_sim=500): + """ + Analyze model selection performance for different sample sizes + """ + errors_f = [] # Type I errors + errors_g = [] # Type II errors + + for T in T_values: + # Simulate from model f + paths_f = simulate_var(model_f, T, N_sim//2) + L_ratios_f = compute_likelihood_ratio_var(paths_f, model_f, model_g) + + # Simulate from model g + paths_g = simulate_var(model_g, T, N_sim//2) + L_ratios_g = compute_likelihood_ratio_var(paths_g, model_f, model_g) + + # Decision rule: choose f if log L_T >= 0 + errors_f.append(np.mean(L_ratios_f[:, -1] < 0)) + errors_g.append(np.mean(L_ratios_g[:, -1] >= 0)) + + return np.array(errors_f), np.array(errors_g) + +T_values = np.arange(1, 50, 1) +errors_f, errors_g = model_selection_analysis(T_values, model2_f, model2_g, N_sim=400) + +fig, ax = plt.subplots() + +ax.plot(T_values, errors_f, 'C0', linewidth=2, label='type I error') +ax.plot(T_values, errors_g, 'C1', linewidth=2, label='type II error') +ax.plot(T_values, 0.5 * (errors_f + errors_g), 'g--', +linewidth=2, label='average error') +ax.set_xlabel('$T$') +ax.set_ylabel('error probability') +ax.set_title('model selection errors') +plt.tight_layout() +plt.show() +``` + +## Application: Samuelson multiplier-accelerator + +Now let's connect to the Samuelson multiplier-accelerator model. + +The model consists of: + +- Consumption: $C_t = \gamma + a Y_{t-1}$ where $a \in (0,1)$ is the marginal propensity to consume +- Investment: $I_t = b(Y_{t-1} - Y_{t-2})$ where $b > 0$ is the accelerator coefficient +- Government spending: $G_t = G$ (constant) + +We have the national income identity + +$$ +Y_t = C_t + I_t + G_t +$$ + +Equations yields the second-order difference equation: + +$$ +Y_t = (\gamma + G) + (a + b)Y_{t-1} - b Y_{t-2} + \sigma \epsilon_t +$$ + +With $\rho_1 = a + b$ and $\rho_2 = -b$, we have: + +$$ +Y_t = (\gamma + G) + \rho_1 Y_{t-1} + \rho_2 Y_{t-2} + \sigma \epsilon_t +$$ + +To fit into our discussion, we write it into state-space representation. + +To handle the constant term properly, we use an augmented state vector $\mathbf{x}_t = [1, Y_t, Y_{t-1}]'$: + +$$ +\mathbf{x}_{t+1} = \begin{bmatrix} +1 \\ +Y_{t+1} \\ +Y_t +\end{bmatrix} = \begin{bmatrix} +1 & 0 & 0 \\ +\gamma + G & \rho_1 & \rho_2 \\ +0 & 1 & 0 +\end{bmatrix} \begin{bmatrix} +1 \\ +Y_t \\ +Y_{t-1} +\end{bmatrix} + \begin{bmatrix} +0 \\ +\sigma \\ +0 +\end{bmatrix} \epsilon_{t+1} +$$ + +The observation equation extracts the economic variables: + +$$ +\mathbf{y}_t = \begin{bmatrix} +Y_t \\ +C_t \\ +I_t +\end{bmatrix} = \begin{bmatrix} +\gamma + G & \rho_1 & \rho_2 \\ +\gamma & a & 0 \\ +0 & b & -b +\end{bmatrix} \begin{bmatrix} +1 \\ +Y_t \\ +Y_{t-1} +\end{bmatrix} +$$ + +This gives us: + +- $Y_t = (\gamma + G) \cdot 1 + \rho_1 Y_{t-1} + \rho_2 Y_{t-2}$ (total output) +- $C_t = \gamma \cdot 1 + a Y_{t-1}$ (consumption) +- $I_t = b(Y_{t-1} - Y_{t-2})$ (investment) + +```{code-cell} ipython3 +def samuelson_to_var(a, b, γ, G, σ): + """ + Convert Samuelson model parameters to VAR form with augmented state + + Samuelson model: + - Y_t = C_t + I_t + G + - C_t = γ + a*Y_{t-1} + - I_t = b*(Y_{t-1} - Y_{t-2}) + + Reduced form: Y_t = (γ+G) + (a+b)*Y_{t-1} - b*Y_{t-2} + σ*ε_t + + State vector is [1, Y_t, Y_{t-1}]' + """ + ρ_1 = a + b + ρ_2 = -b + + # State transition matrix for augmented state + A = np.array([[1, 0, 0], + [γ + G, ρ_1, ρ_2], + [0, 1, 0]]) + + # Shock loading matrix + C = np.array([[0], + [σ], + [0]]) + + # Observation matrix (extracts Y_t, C_t, I_t) + G_obs = np.array([[γ + G, ρ_1, ρ_2], # Y_t + [γ, a, 0], # C_t + [0, b, -b]]) # I_t + + return A, C, G_obs +``` + +We define functions in the code cell below to get the initial conditions and check stability + +```{code-cell} ipython3 +:tags: [hide-input] + +def get_samuelson_initial_conditions(a, b, γ, G, y_0=None, y_m1=None, + stationary_init=False): + """ + Get initial conditions for Samuelson model + """ + # Calculate steady state + y_ss = (γ + G) / (1 - a - b) + + if y_0 is None: + y_0 = y_ss + if y_m1 is None: + y_m1 = y_ss if stationary_init else y_0 * 0.95 + + # Initial mean + μ_0 = np.array([1.0, y_0, y_m1]) + + if stationary_init: + Σ_0 = np.array([[0, 0, 0], + [0, 1, 0.5], + [0, 0.5, 1]]) + else: + Σ_0 = np.array([[0, 0, 0], + [0, 25, 15], + [0, 15, 25]]) + + return μ_0, Σ_0 + +def check_samuelson_stability(a, b): + """ + Check stability of Samuelson model and return characteristic roots + """ + ρ_1 = a + b + ρ_2 = -b + + roots = np.roots([1, -ρ_1, -ρ_2]) + max_abs_root = np.max(np.abs(roots)) + is_stable = max_abs_root < 1 + + # Determine type of dynamics + if np.iscomplex(roots[0]): + if max_abs_root < 1: + dynamics = "Damped oscillations" + else: + dynamics = "Explosive oscillations" + else: + if max_abs_root < 1: + dynamics = "Smooth convergence" + else: + if np.max(roots) > 1: + dynamics = "Explosive growth" + else: + dynamics = "Explosive oscillations (real roots)" + + return is_stable, roots, max_abs_root, dynamics +``` + +Let's implement it and inspect the likelihood ratio processes induced by two Samuelson models with different parameters. + +```{code-cell} ipython3 +def create_samuelson_var_model(a, b, γ, G, σ, stationary_init=False, + y_0=None, y_m1=None): + """ + Create a VAR model from Samuelson parameters + """ + A, C, G_obs = samuelson_to_var(a, b, γ, G, σ) + + μ_0, Σ_0 = get_samuelson_initial_conditions( + a, b, γ, G, y_0, y_m1, stationary_init + ) + + # Create VAR model + model = create_var_model(A, C, μ_0, Σ_0, stationary=False) + is_stable, roots, max_root, dynamics = check_samuelson_stability(a, b) + info = { + 'a': a, 'b': b, 'γ': γ, 'G': G, 'σ': σ, + 'ρ_1': a + b, 'ρ_2': -b, + 'steady_state': (γ + G) / (1 - a - b), + 'is_stable': is_stable, + 'roots': roots, + 'max_abs_root': max_root, + 'dynamics': dynamics + } + + return model, G_obs, info + +def simulate_samuelson(model, G_obs, T, N_paths=1): + """ + Simulate Samuelson model + """ + # Simulate state paths + states = simulate_var(model, T, N_paths) + + # Extract observables using G matrix + if N_paths == 1: + # Single path: states is (T+1, 3) + observables = (G_obs @ states.T).T + else: + # Multiple paths: states is (N_paths, T+1, 3) + observables = np.zeros((N_paths, T+1, 3)) + for i in range(N_paths): + observables[i] = (G_obs @ states[i].T).T + + return states, observables +``` + +Now let's simulate two Samuelson models with different accelerator coefficients and plot their sample paths + +```{code-cell} ipython3 +# Model f: Higher accelerator coefficient +a_f, b_f = 0.98, 0.9 +γ_f, G_f, σ_f = 10, 10, 0.5 + +# Model g: Lower accelerator coefficient +a_g, b_g = 0.98, 0.85 +γ_g, G_g, σ_g = 10, 10, 0.5 + + +model_sam_f, G_obs_f, info_f = create_samuelson_var_model( + a_f, b_f, γ_f, G_f, σ_f, + stationary_init=False, + y_0=100, y_m1=95 +) + +model_sam_g, G_obs_g, info_g = create_samuelson_var_model( + a_g, b_g, γ_g, G_g, σ_g, + stationary_init=False, + y_0=100, y_m1=95 +) + +T = 50 +N_paths = 50 + +# Get both states and observables +states_f, obs_f = simulate_samuelson(model_sam_f, G_obs_f, T, N_paths) +states_g, obs_g = simulate_samuelson(model_sam_g, G_obs_g, T, N_paths) + +output_paths_f = obs_f[:, :, 0] +output_paths_g = obs_g[:, :, 0] + +print("model f:") +print(f" ρ_1 = a + b = {info_f['ρ_1']:.2f}") +print(f" ρ_2 = -b = {info_f['ρ_2']:.2f}") +print(f" roots: {info_f['roots']}") +print(f" dynamics: {info_f['dynamics']}") + +print("\nmodel g:") +print(f" ρ_1 = a + b = {info_g['ρ_1']:.2f}") +print(f" ρ_2 = -b = {info_g['ρ_2']:.2f}") +print(f" roots: {info_g['roots']}") +print(f" dynamics: {info_g['dynamics']}") + + +fig, ax = plt.subplots(1, 1) + +for i in range(min(20, N_paths)): + ax.plot(output_paths_f[i], alpha=0.6, color='C0', linewidth=0.8) + ax.plot(output_paths_g[i], alpha=0.6, color='C1', linewidth=0.8) +ax.set_xlabel('$t$') +ax.set_ylabel('$Y_t$') +ax.legend(['model f', 'model g'], loc='upper left') +plt.tight_layout() +plt.show() +``` + +```{code-cell} ipython3 +# Compute likelihood ratios +L_ratios_ff = compute_likelihood_ratio_var(states_f, model_sam_f, model_sam_g) +L_ratios_gf = compute_likelihood_ratio_var(states_g, model_sam_f, model_sam_g) + +fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + +ax = axes[0] +for i in range(min(20, N_paths)): + ax.plot(L_ratios_ff[i], alpha=0.5, color='C0', lw=0.8) +ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5) +ax.set_title(r'$\log L_t$ (nature = f)') +ax.set_ylabel(r'$\log L_t$') + +ax = axes[1] +for i in range(min(20, N_paths)): + ax.plot(L_ratios_gf[i], alpha=0.5, color='C1', lw=0.8) +ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5) +ax.set_title(r'$\log L_t$ (nature = g)') +plt.show() +``` + +In the figure on the left, data are generated by $f$ and the likelihood ratio diverges to plus infinity. + +In the figure on the right, data are generated by $g$ and the likelihood ratio diverges to negative infinity. + +In both cases, we applied a lower and upper threshold for the log likelihood ratio process for numerical stability since they grow unbounded very quickly. + +In both cases, the likelihood ratio processes eventually lead us to select the correct model. diff --git a/lectures/wald_friedman.md b/lectures/wald_friedman.md index 802f7908d..75b341e32 100644 --- a/lectures/wald_friedman.md +++ b/lectures/wald_friedman.md @@ -66,13 +66,14 @@ We'll begin with some imports: ```{code-cell} ipython3 import numpy as np import matplotlib.pyplot as plt -from numba import njit, prange +from numba import njit, prange, vectorize, jit from numba.experimental import jitclass from math import gamma from scipy.integrate import quad from scipy.stats import beta from collections import namedtuple import pandas as pd +import scipy as sp ``` This lecture uses ideas studied in {doc}`the lecture on likelihood ratio processes` and {doc}`the lecture on Bayesian learning`. @@ -123,7 +124,7 @@ Realizing that, they told Abraham Wald about the problem. That set Wald on a path that led him to create *Sequential Analysis* {cite}`Wald47`. -## Neyman-Pearson Formulation +## Neyman-Pearson formulation It is useful to begin by describing the theory underlying the test that the U.S. Navy told Captain G. S. Schuyler to use. @@ -225,7 +226,7 @@ interpret $\alpha$ and $\beta$: > Thus, we can say that in the long run [ here Wald applies law of > large numbers by driving $M \rightarrow \infty$ (our comment, > not Wald's) ] the proportion of wrong statements will be -> $\alpha$ if $H_0$is true and $\beta$ if +> $\alpha$ if $H_0$ is true and $\beta$ if > $H_1$ is true. The quantity $\alpha$ is called the *size* of the critical region, @@ -275,7 +276,7 @@ Here is how Wald introduces the notion of a sequential test > a random variable, since the value of $n$ depends on the outcome of the > observations. -## Wald's Sequential Formulation +## Wald's sequential formulation By way of contrast to Neyman and Pearson's formulation of the problem, in Wald's formulation @@ -290,38 +291,48 @@ A decision-maker can observe a sequence of draws of a random variable $z$. He (or she) wants to know which of two probability distributions $f_0$ or $f_1$ governs $z$. +We use beta distributions as examples. -To illustrate, let's inspect some beta distributions. - -The density of a Beta probability distribution with parameters $a$ and $b$ is - -$$ -f(z; a, b) = \frac{\Gamma(a+b) z^{a-1} (1-z)^{b-1}}{\Gamma(a) \Gamma(b)} -\quad \text{where} \quad -\Gamma(p) := \int_{0}^{\infty} x^{p-1} e^{-x} dx -$$ - -The next figure shows two beta distributions. +We will also work with Jensen-Shannon divergence introduced in {doc}`divergence_measures`. ```{code-cell} ipython3 -@njit +@vectorize def p(x, a, b): + """Beta distribution density function.""" r = gamma(a + b) / (gamma(a) * gamma(b)) - return r * x**(a-1) * (1 - x)**(b-1) + return r * x** (a-1) * (1 - x) ** (b-1) -f0 = lambda x: p(x, 1, 1) -f1 = lambda x: p(x, 9, 9) -grid = np.linspace(0, 1, 50) +def create_beta_density(a, b): + """Create a beta density function with specified parameters.""" + return jit(lambda x: p(x, a, b)) -fig, ax = plt.subplots(figsize=(10, 8)) +def compute_KL(f, g): + """Compute KL divergence KL(f, g)""" + integrand = lambda w: f(w) * np.log(f(w) / g(w)) + val, _ = quad(integrand, 1e-5, 1-1e-5) + return val -ax.set_title("Original Distributions") +def compute_JS(f, g): + """Compute Jensen-Shannon divergence""" + def m(w): + return 0.5 * (f(w) + g(w)) + + js_div = 0.5 * compute_KL(f, m) + 0.5 * compute_KL(g, m) + return js_div +``` + +The next figure shows two beta distributions + +```{code-cell} ipython3 +f0 = create_beta_density(1, 1) +f1 = create_beta_density(9, 9) +grid = np.linspace(0, 1, 50) + +fig, ax = plt.subplots() ax.plot(grid, f0(grid), lw=2, label="$f_0$") ax.plot(grid, f1(grid), lw=2, label="$f_1$") - ax.legend() ax.set(xlabel="$z$ values", ylabel="probability of $z_k$") - plt.tight_layout() plt.show() ``` @@ -341,7 +352,7 @@ Consequently, the observer has something to learn, namely, whether the observati The decision maker wants to decide which of the two distributions is generating outcomes. -### Type I and Type II Errors +### Type I and type II errors If we regard $f=f_0$ as a null hypothesis and $f=f_1$ as an alternative hypothesis, then @@ -392,7 +403,7 @@ The following figure illustrates aspects of Wald's procedure. ``` -## Links Between $A,B$ and $\alpha, \beta$ +## Links between $A,B$ and $\alpha, \beta$ In chapter 3 of **Sequential Analysis** {cite}`Wald47` Wald establishes the inequalities @@ -448,6 +459,8 @@ We will focus on the case where $f_0$ and $f_1$ are beta distributions since it First, we define a namedtuple to store all the parameters we need for our simulation studies. +We also compute Wald's recommended thresholds $A$ and $B$ based on the target type I and type II errors $\alpha$ and $\beta$ + ```{code-cell} ipython3 SPRTParams = namedtuple('SPRTParams', ['α', 'β', # Target type I and type II errors @@ -455,6 +468,13 @@ SPRTParams = namedtuple('SPRTParams', 'a1', 'b1', # Shape parameters for f_1 'N', # Number of simulations 'seed']) + +@njit +def compute_wald_thresholds(α, β): + """Compute Wald's recommended thresholds.""" + A = (1 - β) / α + B = β / (1 - α) + return A, B, np.log(A), np.log(B) ``` Now we can run the simulation following Wald's recommendation. @@ -490,23 +510,14 @@ def sprt_single_run(a0, b0, a1, b1, logA, logB, true_f0, seed): """Run a single SPRT until a decision is reached.""" log_L = 0.0 n = 0 - - # Set seed for this run np.random.seed(seed) while True: - # Draw a random variable from the appropriate distribution - if true_f0: - z = np.random.beta(a0, b0) - else: - z = np.random.beta(a1, b1) - + z = np.random.beta(a0, b0) if true_f0 else np.random.beta(a1, b1) n += 1 - # Update the log-likelihood ratio - log_f1_z = np.log(p(z, a1, b1)) - log_f0_z = np.log(p(z, a0, b0)) - log_L += log_f1_z - log_f0_z + # Update log-likelihood ratio + log_L += np.log(p(z, a1, b1)) - np.log(p(z, a0, b0)) # Check stopping conditions if log_L >= logA: @@ -516,31 +527,21 @@ def sprt_single_run(a0, b0, a1, b1, logA, logB, true_f0, seed): @njit(parallel=True) def run_sprt_simulation(a0, b0, a1, b1, α, β, N, seed): - """SPRT simulation described by the algorithm.""" + """SPRT simulation.""" + A, B, logA, logB = compute_wald_thresholds(α, β) - # Calculate thresholds - A = (1 - β) / α - B = β / (1 - α) - logA = np.log(A) - logB = np.log(B) - - # Pre-allocate arrays stopping_times = np.zeros(N, dtype=np.int64) - - # Store decision and ground truth as boolean arrays decisions_h0 = np.zeros(N, dtype=np.bool_) truth_h0 = np.zeros(N, dtype=np.bool_) - # Run simulations in parallel for i in prange(N): true_f0 = (i % 2 == 0) truth_h0[i] = true_f0 n, accept_f0 = sprt_single_run( - a0, b0, a1, b1, - logA, logB, - true_f0, seed + i) - + a0, b0, a1, b1, + logA, logB, + true_f0, seed + i) stopping_times[i] = n decisions_h0[i] = accept_f0 @@ -548,7 +549,6 @@ def run_sprt_simulation(a0, b0, a1, b1, α, β, N, seed): def run_sprt(params): """Run SPRT simulations with given parameters.""" - stopping_times, decisions_h0, truth_h0 = run_sprt_simulation( params.a0, params.b0, params.a1, params.b1, params.α, params.β, params.N, params.seed @@ -557,36 +557,27 @@ def run_sprt(params): # Calculate error rates truth_h0_bool = truth_h0.astype(bool) decisions_h0_bool = decisions_h0.astype(bool) - - # For type I error: P(reject H0 | H0 is true) - type_I = np.sum(truth_h0_bool - & ~decisions_h0_bool) / np.sum(truth_h0_bool) - # For type II error: P(accept H0 | H0 is false) - type_II = np.sum(~truth_h0_bool - & decisions_h0_bool) / np.sum(~truth_h0_bool) - - # Create scipy distributions for compatibility - f0 = beta(params.a0, params.b0) - f1 = beta(params.a1, params.b1) + type_I = np.sum(truth_h0_bool & ~decisions_h0_bool) \ + / np.sum(truth_h0_bool) + type_II = np.sum(~truth_h0_bool & decisions_h0_bool) \ + / np.sum(~truth_h0_bool) return { 'stopping_times': stopping_times, 'decisions_h0': decisions_h0_bool, 'truth_h0': truth_h0_bool, 'type_I': type_I, - 'type_II': type_II, - 'f0': f0, - 'f1': f1 + 'type_II': type_II } - + # Run simulation params = SPRTParams(α=0.05, β=0.10, a0=2, b0=5, a1=5, b1=2, N=20000, seed=1) results = run_sprt(params) print(f"Average stopping time: {results['stopping_times'].mean():.2f}") -print(f"Empirical type I error: {results['type_I']:.3f} (target = {params.α})") -print(f"Empirical type II error: {results['type_II']:.3f} (target = {params.β})") +print(f"Empirical type I error: {results['type_I']:.3f} (target = {params.α})") +print(f"Empirical type II error: {results['type_II']:.3f} (target = {params.β})") ``` As anticipated in the passage above in which Wald discussed the quality of @@ -598,133 +589,59 @@ we find that the algorithm actually gives For recent work on the quality of approximation {eq}`eq:Waldrule`, see, e.g., {cite}`fischer2024improving`. ``` -The following code constructs a graph that lets us visualize two distributions and the distribution of times to reach a decision. - -```{code-cell} ipython3 -fig, axes = plt.subplots(1, 2, figsize=(14, 5)) - -z_grid = np.linspace(0, 1, 200) -axes[0].plot(z_grid, results['f0'].pdf(z_grid), 'b-', - lw=2, label=f'$f_0 = \\text{{Beta}}({params.a0},{params.b0})$') -axes[0].plot(z_grid, results['f1'].pdf(z_grid), 'r-', - lw=2, label=f'$f_1 = \\text{{Beta}}({params.a1},{params.b1})$') -axes[0].fill_between(z_grid, 0, - np.minimum(results['f0'].pdf(z_grid), - results['f1'].pdf(z_grid)), - alpha=0.3, color='purple', label='overlap region') -axes[0].set_xlabel('z') -axes[0].set_ylabel('density') -axes[0].legend() - -axes[1].hist(results['stopping_times'], - bins=np.arange(1, results['stopping_times'].max() + 1.5) - 0.5, - color="steelblue", alpha=0.8, edgecolor="black") -axes[1].set_title("distribution of stopping times $n$") -axes[1].set_xlabel("$n$") -axes[1].set_ylabel("frequency") - -plt.show() -``` - -In this example, the stopping time stays below 10. - -We can construct a $2 \times 2$ "confusion matrix" whose diagonal elements -count the number of times that Wald's decision rule correctly accepts and -rejects the null hypothesis. - -```{code-cell} ipython3 -# Accept H0 when H0 is true (correct) -f0_correct = np.sum(results['truth_h0'] & results['decisions_h0']) - -# Reject H0 when H0 is true (incorrect) -f0_incorrect = np.sum(results['truth_h0'] & (~results['decisions_h0'])) - -# Reject H0 when H1 is true (correct) -f1_correct = np.sum((~results['truth_h0']) & (~results['decisions_h0'])) - -# Accept H0 when H1 is true (incorrect) -f1_incorrect = np.sum((~results['truth_h0']) & results['decisions_h0']) - -# First row is when f0 is the true distribution -# Second row is when f1 is true -confusion_data = np.array([[f0_correct, f0_incorrect], - [f1_incorrect, f1_correct]]) - -row_totals = confusion_data.sum(axis=1, keepdims=True) - -print("Confusion Matrix:") -print(confusion_data) - -fig, ax = plt.subplots() -ax.imshow(confusion_data, cmap='Blues', aspect='equal') -ax.set_xticks([0, 1]) -ax.set_xticklabels(['accept $H_0$', 'reject $H_0$']) -ax.set_yticks([0, 1]) -ax.set_yticklabels(['true $f_0$', 'true $f_1$']) - -for i in range(2): - for j in range(2): - percent = confusion_data[i, j] / row_totals[i, 0] if row_totals[i, 0] > 0 else 0 - color = 'white' if confusion_data[i, j] > confusion_data.max() * 0.5 else 'black' - ax.text(j, i, f'{confusion_data[i, j]}\n({percent:.1%})', - ha="center", va="center", - color=color, fontweight='bold') -plt.tight_layout() -plt.show() -``` - -Next we use our code to study three different $f_0, f_1$ pairs having different discrepancies between distributions. - -We plot the same three graphs we used above for each pair of distributions - -```{code-cell} ipython3 -params_1 = SPRTParams(α=0.05, β=0.10, a0=2, b0=8, a1=8, b1=2, N=5000, seed=42) -results_1 = run_sprt(params_1) - -params_2 = SPRTParams(α=0.05, β=0.10, a0=4, b0=5, a1=5, b1=4, N=5000, seed=42) -results_2 = run_sprt(params_2) - -params_3 = SPRTParams(α=0.05, β=0.10, a0=0.5, b0=0.4, a1=0.4, - b1=0.5, N=5000, seed=42) -results_3 = run_sprt(params_3) -``` +The following code creates a few graphs that illustrate the results of our simulation. ```{code-cell} ipython3 :tags: [hide-input] +@njit +def compute_wald_thresholds(α, β): + """Compute Wald's recommended thresholds.""" + A = (1 - β) / α + B = β / (1 - α) + return A, B, np.log(A), np.log(B) + def plot_sprt_results(results, params, title=""): - """Plot SPRT simulation results.""" - fig, axes = plt.subplots(1, 3, figsize=(22, 8)) + """Plot SPRT results.""" + fig, axes = plt.subplots(1, 3, figsize=(20, 6)) # Distribution plots z_grid = np.linspace(0, 1, 200) - axes[0].plot(z_grid, results['f0'].pdf(z_grid), 'b-', lw=2, - label=f'$f_0 = \\text{{Beta}}({params.a0},{params.b0})$') - axes[0].plot(z_grid, results['f1'].pdf(z_grid), 'r-', lw=2, - label=f'$f_1 = \\text{{Beta}}({params.a1},{params.b1})$') + f0 = create_beta_density(params.a0, params.b0) + f1 = create_beta_density(params.a1, params.b1) + + axes[0].plot(z_grid, f0(z_grid), 'b-', lw=2, + label=f'$f_0 = \\text{{Beta}}({params.a0},{params.b0})$') + axes[0].plot(z_grid, f1(z_grid), 'r-', lw=2, + label=f'$f_1 = \\text{{Beta}}({params.a1},{params.b1})$') axes[0].fill_between(z_grid, 0, - np.minimum(results['f0'].pdf(z_grid), results['f1'].pdf(z_grid)), - alpha=0.3, color='purple', label='overlap') + np.minimum(f0(z_grid), f1(z_grid)), + alpha=0.3, color='purple', label='overlap') if title: - axes[0].set_title(title, fontsize=25) - axes[0].set_xlabel('z', fontsize=25) - axes[0].set_ylabel('density', fontsize=25) - axes[0].legend(fontsize=18) - axes[0].tick_params(axis='both', which='major', labelsize=18) + axes[0].set_title(title, fontsize=20) + axes[0].set_xlabel('z', fontsize=16) + axes[0].set_ylabel('density', fontsize=16) + axes[0].legend(fontsize=14) # Stopping times - max_n = max(results['stopping_times'].max(), 101) - bins = np.arange(1, min(max_n, 101)) - 0.5 + max_n = min(results['stopping_times'].max(), 101) + bins = np.arange(1, max_n) - 0.5 axes[1].hist(results['stopping_times'], bins=bins, - color="steelblue", alpha=0.8, edgecolor="black") - axes[1].set_title(f'stopping times (mean={results["stopping_times"].mean():.1f})', - fontsize=25) - axes[1].set_xlabel('n', fontsize=25) - axes[1].set_ylabel('frequency', fontsize=25) + color="steelblue", alpha=0.8, edgecolor="black") + axes[1].set_title(f'stopping times (μ={results["stopping_times"].mean():.1f})', + fontsize=16) + axes[1].set_xlabel('n', fontsize=16) + axes[1].set_ylabel('frequency', fontsize=16) axes[1].set_xlim(0, 100) - axes[1].tick_params(axis='both', which='major', labelsize=18) # Confusion matrix + plot_confusion_matrix(results, axes[2]) + + plt.tight_layout() + plt.show() + +def plot_confusion_matrix(results, ax): + """Plot confusion matrix for SPRT results.""" f0_correct = np.sum(results['truth_h0'] & results['decisions_h0']) f0_incorrect = np.sum(results['truth_h0'] & (~results['decisions_h0'])) f1_correct = np.sum((~results['truth_h0']) & (~results['decisions_h0'])) @@ -734,27 +651,57 @@ def plot_sprt_results(results, params, title=""): [f1_incorrect, f1_correct]]) row_totals = confusion_data.sum(axis=1, keepdims=True) - im = axes[2].imshow(confusion_data, cmap='Blues', aspect='equal') - axes[2].set_title(f'errors: I={results["type_I"]:.3f} '+ - f'II={results["type_II"]:.3f}', fontsize=25) - axes[2].set_xticks([0, 1]) - axes[2].set_xticklabels(['accept $H_0$', 'reject $H_0$'], fontsize=22) - axes[2].set_yticks([0, 1]) - axes[2].set_yticklabels(['true $f_0$', 'true $f_1$'], fontsize=22) - axes[2].tick_params(axis='both', which='major', labelsize=18) - + im = ax.imshow(confusion_data, cmap='Blues', aspect='equal') + ax.set_title(f'errors: I={results["type_I"]:.3f} II={results["type_II"]:.3f}', + fontsize=16) + ax.set_xticks([0, 1]) + ax.set_xticklabels(['accept $H_0$', 'reject $H_0$'], fontsize=14) + ax.set_yticks([0, 1]) + ax.set_yticklabels(['true $f_0$', 'true $f_1$'], fontsize=14) for i in range(2): for j in range(2): - percent = confusion_data[i, j] / row_totals[i, 0] if row_totals[i, 0] > 0 else 0 - color = 'white' if confusion_data[i, j] > confusion_data.max() * 0.5 else 'black' - axes[2].text(j, i, f'{confusion_data[i, j]}\n({percent:.1%})', - ha="center", va="center", - color=color, fontweight='bold', - fontsize=18) + percent = confusion_data[i, j] / row_totals[i, 0] \ + if row_totals[i, 0] > 0 else 0 + color = 'white' if confusion_data[i, j] > confusion_data.max() * 0.5 \ + else 'black' + ax.text(j, i, f'{confusion_data[i, j]}\n({percent:.1%})', + ha="center", va="center", color=color, fontweight='bold', + fontsize=14) +``` - plt.tight_layout() - plt.show() +Let's plot the results of our simulation + +```{code-cell} ipython3 +plot_sprt_results(results, params) +``` + +In this example, the stopping time stays below 10. + +We can construct a $2 \times 2$ "confusion matrix" whose diagonal elements +count the number of times that Wald's decision rule correctly accepts and +rejects the null hypothesis. + +```{code-cell} ipython3 +print("Confusion Matrix data:") +print(f"Type I error: {results['type_I']:.3f}") +print(f"Type II error: {results['type_II']:.3f}") +``` + +Next we use our code to study three different $f_0, f_1$ pairs having different discrepancies between distributions. + +We plot the same three graphs we used above for each pair of distributions + +```{code-cell} ipython3 +params_1 = SPRTParams(α=0.05, β=0.10, a0=2, b0=8, a1=8, b1=2, N=5000, seed=42) +results_1 = run_sprt(params_1) + +params_2 = SPRTParams(α=0.05, β=0.10, a0=4, b0=5, a1=5, b1=4, N=5000, seed=42) +results_2 = run_sprt(params_2) + +params_3 = SPRTParams(α=0.05, β=0.10, a0=0.5, b0=0.4, a1=0.4, + b1=0.5, N=5000, seed=42) +results_3 = run_sprt(params_3) ``` ```{code-cell} ipython3 @@ -789,20 +736,14 @@ That is what we shall do now. We shall compute Jensen-Shannon distance and plot it against the average stopping times. ```{code-cell} ipython3 -def kl_div(h, f): - """KL divergence""" - integrand = lambda w: h(w) * np.log(h(w) / f(w)) - val, _ = quad(integrand, 0, 1) - return val - def js_dist(a0, b0, a1, b1): """Jensen–Shannon distance""" - f0 = lambda w: p(w, a0, b0) - f1 = lambda w: p(w, a1, b1) + f0 = create_beta_density(a0, b0) + f1 = create_beta_density(a1, b1) # Mixture m = lambda w: 0.5*(f0(w) + f1(w)) - return np.sqrt(0.5*kl_div(m, f0) + 0.5*kl_div(m, f1)) + return np.sqrt(0.5*compute_KL(m, f0) + 0.5*compute_KL(m, f1)) def generate_β_pairs(N=100, T=10.0, d_min=0.5, d_max=9.5): ds = np.linspace(d_min, d_max, N) @@ -831,11 +772,10 @@ for a0, b0, a1, b1 in param_comb: param_list.append((a0, b0, a1, b1)) # Create the plot -fig, ax = plt.subplots(figsize=(6, 6)) +fig, ax = plt.subplots() scatter = ax.scatter(js_dists, mean_stopping_times, - s=80, alpha=0.7, c=range(len(js_dists)), - linewidth=0.5) + s=80, alpha=0.7, linewidth=0.5) ax.set_xlabel('Jensen–Shannon distance', fontsize=14) ax.set_ylabel('mean stopping time', fontsize=14) @@ -851,47 +791,41 @@ As Jensen-Shannon divergence increases (distributions become more separated), t Below are sampled examples from the experiments we have above ```{code-cell} ipython3 -selected_indices = [0, - len(param_comb)//6, - len(param_comb)//3, - len(param_comb)//2, - 2*len(param_comb)//3, - -1] - -fig, axes = plt.subplots(2, 3, figsize=(15, 8)) - -for i, idx in enumerate(selected_indices): - row = i // 3 - col = i % 3 +def plot_beta_distributions_grid(param_list, js_dists, mean_stopping_times, + selected_indices=None): + """Plot grid of beta distributions with JS distance and stopping times.""" + if selected_indices is None: + selected_indices = [0, len(param_list)//6, len(param_list)//3, + len(param_list)//2, 2*len(param_list)//3, -1] - a0, b0, a1, b1 = param_list[idx] - js_dist = js_dists[idx] - mean_time = mean_stopping_times[idx] - - # Plot the distributions + fig, axes = plt.subplots(2, 3, figsize=(15, 8)) z_grid = np.linspace(0, 1, 200) - f0_dist = beta(a0, b0) - f1_dist = beta(a1, b1) - - axes[row, col].plot(z_grid, f0_dist.pdf(z_grid), 'b-', - lw=2, label='$f_0$') - axes[row, col].plot(z_grid, f1_dist.pdf(z_grid), 'r-', - lw=2, label='$f_1$') - axes[row, col].fill_between(z_grid, 0, - np.minimum(f0_dist.pdf(z_grid), - f1_dist.pdf(z_grid)), - alpha=0.3, color='purple') - axes[row, col].set_title(f'JS dist: {js_dist:.3f}' - +f'\nMean time: {mean_time:.1f}', fontsize=12) - axes[row, col].set_xlabel('z', fontsize=10) - if i == 0: - axes[row, col].set_ylabel('density', fontsize=10) - axes[row, col].legend(fontsize=10) + for i, idx in enumerate(selected_indices): + row, col = i // 3, i % 3 + a0, b0, a1, b1 = param_list[idx] + + f0 = create_beta_density(a0, b0) + f1 = create_beta_density(a1, b1) + + axes[row, col].plot(z_grid, f0(z_grid), 'b-', lw=2, label='$f_0$') + axes[row, col].plot(z_grid, f1(z_grid), 'r-', lw=2, label='$f_1$') + axes[row, col].fill_between(z_grid, 0, + np.minimum(f0(z_grid), f1(z_grid)), + alpha=0.3, color='purple') + + axes[row, col].set_title(f'JS dist: {js_dists[idx]:.3f}' + f'\nMean time: {mean_stopping_times[idx]:.1f}', + fontsize=12) + axes[row, col].set_xlabel('z', fontsize=10) + if i == 0: + axes[row, col].set_ylabel('density', fontsize=10) + axes[row, col].legend(fontsize=10) + plt.tight_layout() + plt.show() -plt.tight_layout() -plt.show() +plot_beta_distributions_grid(param_list, js_dists, mean_stopping_times) ``` Again, we find that the stopping time is shorter when the distributions are more separated, as @@ -901,18 +835,14 @@ Let's visualize individual likelihood ratio processes to see how they evolve tow ```{code-cell} ipython3 def plot_likelihood_paths(params, n_highlight=10, n_background=200): - """Plot likelihood ratio paths""" - - A = (1 - params.β) / params.α - B = params.β / (1 - params.α) - logA, logB = np.log(A), np.log(B) - - f0 = beta(params.a0, params.b0) - f1 = beta(params.a1, params.b1) + """visualize likelihood ratio paths.""" + A, B, logA, logB = compute_wald_thresholds(params.α, params.β) + f0, f1 = map(lambda ab: create_beta_density(*ab), + [(params.a0, params.b0), + (params.a1, params.b1)]) fig, axes = plt.subplots(1, 2, figsize=(14, 7)) - # Generate and plot paths for each distribution for dist_idx, (true_f0, ax, title) in enumerate([ (True, axes[0], 'true distribution: $f_0$'), (False, axes[1], 'true distribution: $f_1$') @@ -920,47 +850,42 @@ def plot_likelihood_paths(params, n_highlight=10, n_background=200): rng = np.random.default_rng(seed=42 + dist_idx) paths_data = [] + # Generate paths for path in range(n_background + n_highlight): - log_L_path = [0.0] # Start at 0 - log_L = 0.0 - n = 0 + log_L_path, log_L, n = [0.0], 0.0, 0 while True: - z = f0.rvs(random_state=rng) if true_f0 else f1.rvs(random_state=rng) + z = rng.beta(params.a0, params.b0) if true_f0 \ + else rng.beta(params.a1, params.b1) n += 1 - log_L += np.log(f1.pdf(z)) - np.log(f0.pdf(z)) + log_L += np.log(f1(z)) - np.log(f0(z)) log_L_path.append(log_L) - # Check stopping conditions if log_L >= logA or log_L <= logB: - # True = reject H0, False = accept H0 - decision = log_L >= logA + paths_data.append((log_L_path, n, log_L >= logA)) break - - paths_data.append((log_L_path, n, decision)) - for i, (path, n, decision) in enumerate(paths_data[:n_background]): - color = 'C1' if decision else 'C0' - ax.plot(range(len(path)), path, - color=color, alpha=0.2, linewidth=0.5) + # Plot background paths + for path, _, decision in paths_data[:n_background]: + ax.plot(range(len(path)), path, color='C1' if decision else 'C0', + alpha=0.2, linewidth=0.5) - for i, (path, n, decision) in enumerate(paths_data[n_background:]): - # Color code by decision - color = 'C1' if decision else 'C0' - ax.plot(range(len(path)), path, color=color, - alpha=0.8, linewidth=1.5, - label='reject $H_0$' if decision and i == 0 else ( - 'accept $H_0$' if not decision and i == 0 else '')) + # Plot highlighted paths with labels + for i, (path, _, decision) in enumerate(paths_data[n_background:]): + ax.plot(range(len(path)), path, color='C1' if decision else 'C0', + alpha=0.8, linewidth=1.5, + label='reject $H_0$' if decision and i == 0 else ( + 'accept $H_0$' if not decision and i == 0 else '')) + # Add threshold lines and formatting ax.axhline(y=logA, color='C1', linestyle='--', linewidth=2, label=f'$\\log A = {logA:.2f}$') ax.axhline(y=logB, color='C0', linestyle='--', linewidth=2, label=f'$\\log B = {logB:.2f}$') - ax.axhline(y=0, color='black', linestyle='-', - alpha=0.5, linewidth=1) + ax.axhline(y=0, color='black', linestyle='-', alpha=0.5, linewidth=1) - ax.set_xlabel(r'$n$') - ax.set_ylabel(r'$log(L_n)$') + ax.set_xlabel(r'$n$') + ax.set_ylabel(r'$\log(L_n)$') ax.set_title(title, fontsize=20) ax.legend(fontsize=18, loc='center right') @@ -1072,7 +997,7 @@ This increases the probability of Type II errors. The table confirms this intuition: as $A$ decreases and $B$ increases from their optimal Wald values, both Type I and Type II error rates increase, while the mean stopping time decreases. -## Related Lectures +## Related lectures We'll dig deeper into some of the ideas used here in the following earlier and later lectures: @@ -1081,3 +1006,368 @@ We'll dig deeper into some of the ideas used here in the following earlier and l * For a deeper understanding of likelihood ratio processes and their role in frequentist and Bayesian statistical theories, see {doc}`likelihood_ratio_process`. * Building on that foundation, {doc}`likelihood_bayes` examines the role of likelihood ratio processes in **Bayesian learning**. * Finally, {doc}`this later lecture ` revisits the subject discussed here and examines whether the frequentist decision rule that the Navy ordered the captain to use would perform better or worse than Abraham Wald's sequential decision rule. + +## Exercises + +In the two exercises below, please try to rewrite the entire SPRT suite in this lecture. + +```{exercise} +:label: wald_friedman_ex1 + +In the first exercise, we apply the sequential probability ratio test to distinguish two models generated by 3-state Markov chains + +(For a review on likelihood ratio processes for Markov chains, see [this section](lrp_markov).) + +Consider distinguishing between two 3-state Markov chain models using Wald's sequential probability ratio test. + +You have competing hypotheses about the transition probabilities: + +- $H_0$: The chain follows transition matrix $P^{(0)}$ +- $H_1$: The chain follows transition matrix $P^{(1)}$ + +Given transition matrices: + +$$ +P^{(0)} = \begin{bmatrix} +0.7 & 0.2 & 0.1 \\ +0.3 & 0.5 & 0.2 \\ +0.1 & 0.3 & 0.6 +\end{bmatrix}, \quad +P^{(1)} = \begin{bmatrix} +0.5 & 0.3 & 0.2 \\ +0.2 & 0.6 & 0.2 \\ +0.2 & 0.2 & 0.6 +\end{bmatrix} +$$ + +For a sequence of observations $(x_0, x_1, \ldots, x_t)$, the likelihood ratio is: + +$$ +\Lambda_t = \frac{\pi_{x_0}^{(1)}}{\pi_{x_0}^{(0)}} \prod_{s=1}^t \frac{P_{x_{s-1},x_s}^{(1)}}{P_{x_{s-1},x_s}^{(0)}} +$$ + +where $\pi^{(i)}$ is the stationary distribution under hypothesis $i$. + +Tasks: +1. Implement the likelihood ratio computation for Markov chains +2. Implement Wald's sequential test with Type I error $\alpha = 0.05$ and Type II error $\beta = 0.10$ +3. Run 1000 simulations under each hypothesis and compute empirical error rates +4. Analyze the distribution of stopping times + +The test stops when: +- $\Lambda_t \geq A = \frac{1-\beta}{\alpha} = 18$: Reject $H_0$ +- $\Lambda_t \leq B = \frac{\beta}{1-\alpha} = 0.105$: Accept $H_0$ +``` + + +```{solution-start} wald_friedman_ex1 +:class: dropdown +``` + +Below is one solution to the exercise. + +In the lecture, we write the code more verbosely to illustrate the concepts clearly. + +In the code below, we simplified some of the code structure for a shorter presentation. + +First we define the parameters for the Markov chain SPRT + +```{code-cell} ipython3 +MarkovSPRTParams = namedtuple('MarkovSPRTParams', + ['α', 'β', 'P_0', 'P_1', 'N', 'seed']) + +def compute_stationary_distribution(P): + """Compute stationary distribution of transition matrix P.""" + eigenvalues, eigenvectors = np.linalg.eig(P.T) + idx = np.argmin(np.abs(eigenvalues - 1)) + pi = np.real(eigenvectors[:, idx]) + return pi / pi.sum() + +@njit +def simulate_markov_chain(P, pi_0, T, seed): + """Simulate a Markov chain path.""" + np.random.seed(seed) + path = np.zeros(T, dtype=np.int32) + + cumsum_pi = np.cumsum(pi_0) + path[0] = np.searchsorted(cumsum_pi, np.random.uniform()) + + for t in range(1, T): + cumsum_row = np.cumsum(P[path[t-1]]) + path[t] = np.searchsorted(cumsum_row, np.random.uniform()) + + return path +``` + +Here we define the function that runs SPRT for Markov chains + +```{code-cell} ipython3 +@njit +def markov_sprt_single_run(P_0, P_1, π_0, π_1, + logA, logB, true_P, true_π, seed): + """Run single SPRT for Markov chains.""" + max_n = 10000 + path = simulate_markov_chain(true_P, true_π, max_n, seed) + + log_L = np.log(π_1[path[0]] / π_0[path[0]]) + if log_L >= logA: return 1, False + if log_L <= logB: return 1, True + + for t in range(1, max_n): + prev_state, curr_state = path[t-1], path[t] + p_1, p_0 = P_1[prev_state, curr_state], P_0[prev_state, curr_state] + + if p_0 > 0: + log_L += np.log(p_1 / p_0) + elif p_1 > 0: + log_L = np.inf + + if log_L >= logA: return t+1, False + if log_L <= logB: return t+1, True + + return max_n, log_L < 0 + +def run_markov_sprt(params): + """Run SPRT for Markov chains.""" + π_0 = compute_stationary_distribution(params.P_0) + π_1 = compute_stationary_distribution(params.P_1) + A, B, logA, logB = compute_wald_thresholds(params.α, params.β) + + stopping_times = np.zeros(params.N, dtype=np.int64) + decisions_h0 = np.zeros(params.N, dtype=bool) + truth_h0 = np.zeros(params.N, dtype=bool) + + for i in range(params.N): + true_P, true_π = (params.P_0, π_0) if i % 2 == 0 else (params.P_1, π_1) + truth_h0[i] = i % 2 == 0 + + n, accept_h0 = markov_sprt_single_run( + params.P_0, params.P_1, π_0, π_1, logA, logB, + true_P, true_π, params.seed + i) + + stopping_times[i] = n + decisions_h0[i] = accept_h0 + + type_I = np.sum(truth_h0 & ~decisions_h0) / np.sum(truth_h0) + type_II = np.sum(~truth_h0 & decisions_h0) / np.sum(~truth_h0) + + return { + 'stopping_times': stopping_times, 'decisions_h0': decisions_h0, + 'truth_h0': truth_h0, 'type_I': type_I, 'type_II': type_II + } +``` + +Now we can run the SPRT for the Markov chain models and visualize the results + +```{code-cell} ipython3 +# Run Markov chain SPRT +P_0 = np.array([[0.7, 0.2, 0.1], + [0.3, 0.5, 0.2], + [0.1, 0.3, 0.6]]) + +P_1 = np.array([[0.5, 0.3, 0.2], + [0.2, 0.6, 0.2], + [0.2, 0.2, 0.6]]) + +params_markov = MarkovSPRTParams(α=0.05, β=0.10, + P_0=P_0, P_1=P_1, N=1000, seed=42) +results_markov = run_markov_sprt(params_markov) + + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5)) + +ax1.hist(results_markov['stopping_times'], + bins=50, color="steelblue", alpha=0.8) +ax1.set_title("stopping times") +ax1.set_xlabel("n") +ax1.set_ylabel("frequency") + +plot_confusion_matrix(results_markov, ax2) + +plt.tight_layout() +plt.show() +``` + +```{solution-end} +``` + + +```{exercise} +:label: wald_friedman_ex2 + +In this exercise, apply Wald's sequential test to distinguish between two VAR(1) models with different dynamics and noise structures. + +For a review of the likelihood ratio process with VAR models, see {doc}`likelihood_var`. + +Given VAR models under each hypothesis: +- $H_0$: $x_{t+1} = A^{(0)} x_t + C^{(0)} w_{t+1}$ +- $H_1$: $x_{t+1} = A^{(1)} x_t + C^{(1)} w_{t+1}$ + +where $w_t \sim \mathcal{N}(0, I)$ and: + +$$ +A^{(0)} = \begin{bmatrix} 0.8 & 0.1 \\ 0.2 & 0.7 \end{bmatrix}, \quad +C^{(0)} = \begin{bmatrix} 0.3 & 0.1 \\ 0.1 & 0.3 \end{bmatrix} +$$ + +$$ +A^{(1)} = \begin{bmatrix} 0.6 & 0.2 \\ 0.3 & 0.5 \end{bmatrix}, \quad +C^{(1)} = \begin{bmatrix} 0.4 & 0 \\ 0 & 0.4 \end{bmatrix} +$$ + +Tasks: +1. Implement the VAR likelihood ratio using the functions from the VAR lecture +2. Implement Wald's sequential test with $\alpha = 0.05$ and $\beta = 0.10$ +3. Analyze performance under both hypotheses and with model misspecification +4. Compare with the Markov chain case in terms of stopping times and accuracy + +``` + +```{solution-start} wald_friedman_ex2 +:class: dropdown +``` + +Below is one solution to the exercise. + +First we define the parameters for the VAR models and simulator + +```{code-cell} ipython3 +VARSPRTParams = namedtuple('VARSPRTParams', + ['α', 'β', 'A_0', 'C_0', 'A_1', 'C_1', 'N', 'seed']) + +def create_var_model(A, C): + """Create VAR model.""" + μ_0 = np.zeros(A.shape[0]) + CC = C @ C.T + Σ_0 = sp.linalg.solve_discrete_lyapunov(A, CC) + + CC_inv = np.linalg.inv(CC + 1e-10 * np.eye(CC.shape[0])) + Σ_0_inv = np.linalg.inv(Σ_0 + 1e-10 * np.eye(Σ_0.shape[0])) + + return { + 'A': A, 'C': C, 'μ_0': μ_0, 'Σ_0': Σ_0, + 'CC_inv': CC_inv, 'Σ_0_inv': Σ_0_inv, + 'log_det_CC': np.log( + np.linalg.det(CC + 1e-10 * np.eye(CC.shape[0]))), + 'log_det_Σ_0': np.log( + np.linalg.det(Σ_0 + 1e-10 * np.eye(Σ_0.shape[0]))) + } +``` + +Now we define the likelihood ratio for the VAR models and the SPRT function similar to the +Markov chain case + +```{code-cell} ipython3 +def var_log_likelihood(x_curr, x_prev, model, initial=False): + """Compute VAR log-likelihood.""" + n = len(x_curr) + if initial: + diff = x_curr - model['μ_0'] + return -0.5 * (n * np.log(2 * np.pi) + model['log_det_Σ_0'] + + diff @ model['Σ_0_inv'] @ diff) + else: + diff = x_curr - model['A'] @ x_prev + return -0.5 * (n * np.log(2 * np.pi) + model['log_det_CC'] + + diff @ model['CC_inv'] @ diff) + +def var_sprt_single_run(model_0, model_1, model_true, + logA, logB, seed): + """Single VAR SPRT run.""" + np.random.seed(seed) + max_T = 500 + + # Generate VAR path + Σ_chol = np.linalg.cholesky(model_true['Σ_0']) + x = model_true['μ_0'] + Σ_chol @ np.random.randn( + len(model_true['μ_0'])) + + # Initial likelihood ratio + log_L = (var_log_likelihood(x, None, model_1, True) - + var_log_likelihood(x, None, model_0, True)) + + if log_L >= logA: return 1, False + if log_L <= logB: return 1, True + + # Sequential updates + for t in range(1, max_T): + x_prev = x.copy() + w = np.random.randn(model_true['C'].shape[1]) + x = model_true['A'] @ x + model_true['C'] @ w + + log_L += (var_log_likelihood(x, x_prev, model_1) - + var_log_likelihood(x, x_prev, model_0)) + + if log_L >= logA: return t+1, False + if log_L <= logB: return t+1, True + + return max_T, log_L < 0 + +def run_var_sprt(params): + """Run VAR SPRT.""" + + model_0 = create_var_model(params.A_0, params.C_0) + model_1 = create_var_model(params.A_1, params.C_1) + A, B, logA, logB = compute_wald_thresholds(params.α, params.β) + + stopping_times = np.zeros(params.N) + decisions_h0 = np.zeros(params.N, dtype=bool) + truth_h0 = np.zeros(params.N, dtype=bool) + + for i in range(params.N): + model_true = model_0 if i % 2 == 0 else model_1 + truth_h0[i] = i % 2 == 0 + + n, accept_h0 = var_sprt_single_run(model_0, model_1, model_true, + logA, logB, params.seed + i) + stopping_times[i] = n + decisions_h0[i] = accept_h0 + + type_I = np.sum(truth_h0 & ~decisions_h0) / np.sum(truth_h0) + type_II = np.sum(~truth_h0 & decisions_h0) / np.sum(~truth_h0) + + return {'stopping_times': stopping_times, + 'decisions_h0': decisions_h0, + 'truth_h0': truth_h0, + 'type_I': type_I, 'type_II': type_II} +``` + +Let's run SPRT and visualize the results + +```{code-cell} ipython3 +# Run VAR SPRT +A_0 = np.array([[0.8, 0.1], + [0.2, 0.7]]) +C_0 = np.array([[0.3, 0.1], + [0.1, 0.3]]) +A_1 = np.array([[0.6, 0.2], + [0.3, 0.5]]) +C_1 = np.array([[0.4, 0.0], + [0.0, 0.4]]) + +params_var = VARSPRTParams(α=0.05, β=0.10, + A_0=A_0, C_0=C_0, A_1=A_1, C_1=C_1, + N=1000, seed=42) +results_var = run_var_sprt(params_var) + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) + +ax1.boxplot([results_markov['stopping_times'], + results_var['stopping_times']], + tick_labels=['Markov Chain', 'VAR(1)']) +ax1.set_ylabel('stopping time') + +x = np.arange(2) +ax2.bar(x - 0.2, [results_markov['type_I'], results_var['type_I']], + 0.4, label='Type I', alpha=0.7) +ax2.bar(x + 0.2, [results_markov['type_II'], results_var['type_II']], + 0.4, label='Type II', alpha=0.7) +ax2.axhline(y=0.05, linestyle='--', alpha=0.5, color='C0') +ax2.axhline(y=0.10, linestyle='--', alpha=0.5, color='C1') +ax2.set_xticks(x), ax2.set_xticklabels(['Markov', 'VAR']) +ax2.legend() +plt.tight_layout() +plt.show() +``` + +```{solution-end} +``` \ No newline at end of file diff --git a/lectures/wald_friedman_2.md b/lectures/wald_friedman_2.md index 6d5142dc6..ec7937cac 100644 --- a/lectures/wald_friedman_2.md +++ b/lectures/wald_friedman_2.md @@ -266,7 +266,7 @@ parameter such as $c$ or $L_0$ on $A$ or $B$. Let $J(\pi)$ be the total loss for a decision-maker with current belief $\pi$ who chooses optimally. -With some thought, you will agree that $J$ should satisfy the Bellman equation +Principles of **dynamic programming** teach us that an optimal loss function $J$ satisfies the following the Bellman functional equation ```{math} :label: new1 @@ -340,7 +340,7 @@ $$ Our aim is to compute the cost function $J$ as well as the associated cutoffs $A$ and $B$. -To make our computations manageable, using {eq}`optdec`, we can write the continuation cost $h(\pi)$ as +To help make our computations more manageable, we can use {eq}`optdec` to write the continuation cost $h(\pi)$ as ```{math} :label: optdec2