import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

rng = np.random.default_rng(777)
plt.rcParams['figure.figsize'] = (8, 4.5)

url = 'https://vincentarelbundock.github.io/Rdatasets/csv/sampleSelection/Mroz87.csv'
try:
    mroz = pd.read_csv(url)                       # Rdatasets Mroz87 (needs a network)
except Exception:
    # Offline fallback: the identical Mroz (1987) data ships with the `wooldridge` package.
    import wooldridge
    mroz = (wooldridge.data('mroz')
            .rename(columns={'inlf': 'lfp', 'kidslt6': 'kids5', 'kidsge6': 'kids618'})
            .copy())
print(f'Number of observations: {len(mroz)}')
print(f'Working women (lfp=1):  {(mroz["lfp"]==1).sum()}')
mroz.head(3)

Number of observations: 753
Working women (lfp=1):  428

mroz_w = mroz.loc[mroz['lfp'] == 1].copy()
mroz_w['lwage'] = np.log(mroz_w['wage'])
print(f'I = {len(mroz_w)}')
print(mroz_w[['wage', 'lwage', 'educ', 'exper']].describe().round(2))

I = 428
         wage   lwage    educ   exper
count  428.00  428.00  428.00  428.00
mean     4.18    1.19   12.66   13.04
std      3.31    0.72    2.29    8.06
min      0.13   -2.05    5.00    0.00
25%      2.26    0.82   12.00    7.00
50%      3.48    1.25   12.00   12.00
75%      4.97    1.60   14.00   18.00
max     25.00    3.22   17.00   38.00

def clt_simulation(sampler, mu, sigma, I=200, R=20_000):
    X_r_i = sampler(size=(R, I))     # R Monte Carlo samples (r in [R]) of I individuals (i in [I])
    mean_r = X_r_i.mean(axis=1)
    Z_r = np.sqrt(I) * (mean_r - mu) / sigma
    return Z_r

specs_clt = {
    'Uniform(0,1)':       (lambda size: rng.random(size), 0.5,  np.sqrt(1/12)),
    'Exponential(1)':     (lambda size: rng.exponential(1.0, size), 1.0, 1.0),
    'Bernoulli(0.1)':     (lambda size: (rng.random(size) < 0.1).astype(float), 0.1, np.sqrt(0.1*0.9)),
    'Lognormal(0,1)':     (lambda size: rng.lognormal(0.0, 1.0, size), np.exp(0.5), np.sqrt((np.exp(1)-1)*np.exp(1))),
}

I = 200
R = 20_000
S = len(specs_clt)
Z_s_r = np.empty((S, R))
for s, (_, (sampler, mu, sigma)) in enumerate(specs_clt.items()):
    Z_s_r[s] = clt_simulation(sampler, mu, sigma, I=I, R=R)

z_g = np.linspace(-4, 4, 400)

fig, axes = plt.subplots(2, 2, figsize=(11, 6.5))
for s, (ax, name) in enumerate(zip(axes.flat, specs_clt)):
    ax.hist(Z_s_r[s], bins=60, density=True, alpha=0.55, edgecolor='white')
    ax.plot(z_g, stats.norm.pdf(z_g), 'k-', lw=1.5, label=r'$\mathcal{N}(0,1)$')
    ax.set_xlim(-4, 4)
    ax.set_title(f'{name},  I = {I}')
    ax.legend(loc='upper right', fontsize=9)
plt.tight_layout()
plt.show()

# Berry-Esseen bounds the Kolmogorov distance  sup_x |P(Z_I <= x) - Phi(x)|  by
# C * rho / (sigma**3 * sqrt(I)),  with rho = E|X - mu|**3 the third absolute moment.
# We estimate that sup by the Kolmogorov-Smirnov (KS) distance below, then check that
# it decays at the predicted I^(-1/2) rate (a straight line of slope -1/2 on log-log axes).

def ks_distance(Z_r):
    """Two-sided Kolmogorov distance  sup_x |F_emp(x) - Phi(x)|.

    The empirical CDF is a step function, so at each order statistic z_(r) the gap has to be
    measured on BOTH sides of the jump: r/R - Phi(z_(r)) just after it, and Phi(z_(r)) -
    (r-1)/R just before it. Taking only the first understates the supremum.
    """
    Z_sorted_r = np.sort(Z_r)
    R = len(Z_sorted_r)
    F_th_r = stats.norm.cdf(Z_sorted_r)
    d_plus  = np.max(np.arange(1, R + 1) / R - F_th_r)      # above the step
    d_minus = np.max(F_th_r - np.arange(0, R) / R)          # below the step
    return max(d_plus, d_minus)

# Cross-check the implementation against scipy's two-sided one-sample KS statistic.
_z_check = rng.standard_normal(4_000)
assert abs(ks_distance(_z_check)
           - stats.kstest(_z_check, 'norm').statistic) < 1e-12, 'KS distance mismatch'
print('ks_distance agrees with scipy.stats.kstest to machine precision')

I_s = np.array([10, 30, 100, 300, 1000, 3000, 10_000])   # grid of sample sizes, s in [S]
R = 5_000
sampler = lambda size: rng.exponential(1.0, size)
mu, sigma = 1.0, 1.0

ks_s = np.empty(len(I_s))
for s, I in enumerate(I_s):
    Z_r = clt_simulation(sampler, mu, sigma, I=I, R=R)
    ks_s[s] = ks_distance(Z_r)

ks_distance agrees with scipy.stats.kstest to machine precision

fig, ax = plt.subplots(figsize=(7, 4.2))
ax.loglog(I_s, ks_s, 'o-', label=r'KS distance to $\mathcal{N}(0,1)$')
ax.loglog(I_s, ks_s[0] * np.sqrt(I_s[0] / I_s), 'k--', label=r'reference $I^{-1/2}$')
ax.set_xlabel('I')
ax.set_ylabel('KS distance')
ax.set_title('Berry–Esseen rate: parent = Exponential(1)')
ax.legend()
plt.tight_layout()
plt.show()

I_s = np.unique(np.round(np.logspace(1, 4, 30)).astype(int))   # grid of sample sizes, s in [S]
R = 5
sigma2_true = (np.exp(1) - 1) * np.exp(1)
s2_r_s = np.empty((R, len(I_s)))
for r in range(R):
    x_i = rng.lognormal(0.0, 1.0, size=I_s.max())
    s2_r_s[r] = np.array([x_i[:I].var(ddof=0) for I in I_s])

fig, ax = plt.subplots(figsize=(8, 4.2))
for r in range(R):
    ax.plot(I_s, s2_r_s[r], lw=0.9, alpha=0.85)
ax.axhline(sigma2_true, color='k', ls='--', label=fr'$\sigma^2 = {sigma2_true:.4f}$')
ax.set_xscale('log')
ax.set_xlabel('I (log scale)')
ax.set_ylabel(r'$S_I^2$')
ax.set_title('Sample variance of Lognormal(0,1) draws')
ax.legend()
plt.tight_layout()
plt.show()

def t_simulation(sampler, mu, I=20, R=20_000):
    X_r_i = sampler(size=(R, I))
    mean_r = X_r_i.mean(axis=1)
    sd_r = X_r_i.std(axis=1, ddof=1)
    T_r = np.sqrt(I) * (mean_r - mu) / sd_r
    return T_r

I = 20
R = 20_000
sampler = lambda size: rng.exponential(1.0, size)
T_r = t_simulation(sampler, mu=1.0, I=I, R=R)

# Empirical CDF vs N(0,1) and t_{I-1}
T_sorted_r = np.sort(T_r)
F_emp_r = np.arange(1, R+1) / R
z_g = np.linspace(-5, 5, 400)

# Q-Q plot vs N(0,1)
u_q = np.linspace(0.005, 0.995, 200)      # quantile levels, q in [Q]
qemp_q = np.quantile(T_r, u_q)
qth_q = stats.norm.ppf(u_q)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.2))

axes[0].plot(T_sorted_r, F_emp_r, lw=1.4, label=f'empirical CDF, I={I}')
axes[0].plot(z_g, stats.norm.cdf(z_g), 'k--', lw=1.0, label=r'$\mathcal{N}(0,1)$')
axes[0].plot(z_g, stats.t.cdf(z_g, df=I-1), 'r:', lw=1.2, label=fr'$t_{{{I-1}}}$')
axes[0].set_xlim(-5, 5)
axes[0].set_xlabel(r'$T_I$')
axes[0].set_ylabel('CDF')
axes[0].legend(loc='lower right')
axes[0].set_title('CDF of the studentized mean')

axes[1].plot(qth_q, qemp_q, 'o', ms=2.5)
axes[1].plot([-3, 3], [-3, 3], 'k--', lw=1.0)
axes[1].set_xlabel(r'theoretical $\mathcal{N}(0,1)$ quantiles')
axes[1].set_ylabel(r'empirical $T_I$ quantiles')
axes[1].set_title(fr'Q–Q plot of $T_I$ vs $\mathcal{{N}}(0,1)$, I = {I}')
plt.tight_layout()
plt.show()

x_i = mroz_w['lwage'].values
I = len(x_i)
xbar = x_i.mean()
s = x_i.std(ddof=1)
se = s / np.sqrt(I)

alpha = 0.05
z = stats.norm.ppf(1 - alpha/2)
t_crit = stats.t.ppf(1 - alpha/2, df=I-1)

ci_asy = (xbar - z*se, xbar + z*se)
ci_exact = (xbar - t_crit*se, xbar + t_crit*se)

print(f'  I          = {I}')
print(f'  mean lwage = {xbar:.4f}')
print(f'  s.e.       = {se:.4f}')
print(f'  asymptotic 95% CI (CLT+Slutsky) = [{ci_asy[0]:.4f}, {ci_asy[1]:.4f}]')
print(f'  exact      95% CI (t_{I-1})       = [{ci_exact[0]:.4f}, {ci_exact[1]:.4f}]')

  I          = 428
  mean lwage = 1.1902
  s.e.       = 0.0350
  asymptotic 95% CI (CLT+Slutsky) = [1.1217, 1.2587]
  exact      95% CI (t_427)       = [1.1215, 1.2589]

	lfp	hours	kids5	kids618	age	educ	wage	repwage	hushrs	husage	...	faminc	mtr	motheduc	fatheduc	unem	city	exper	nwifeinc	lwage	expersq
0	1	1610	1	0	32	12	3.3540	2.65	2708	34	...	16310.0	0.7215	12	7	5.0	0	14	10.910060	1.210154	196
1	1	1656	0	2	30	12	1.3889	2.65	2310	30	...	21800.0	0.6615	7	7	11.0	1	5	19.499981	0.328512	25
2	1	1980	1	3	35	12	4.5455	4.04	3072	40	...	21040.0	0.6915	12	7	5.0	0	15	12.039910	1.514138	225

Asymptotic foundations¶

Alfred Galichon (NYU)¶

'math+econ+code' masterclass series: econometrics¶

CLT, CMT, Slutsky, and a first confidence interval¶

Learning objectives¶

References¶

Motivation: completing the asymptotic toolkit¶

Loading our libraries¶

Loading our data¶

1. The Central Limit Theorem¶

1.1 Statement¶

1.2 Simulation across parent distributions¶

1.3 Rate of convergence: the Berry–Esseen theorem¶

2. The Continuous Mapping Theorem¶

2.1 Application: consistency of the sample variance¶

3. Slutsky's theorem¶

3.1 The studentized sample mean¶

3.2 Simulation¶

4. Synthesis: an asymptotic confidence interval for the mean log-wage¶

Looking ahead¶