Econometrics I¶

TA Christian Alemán

Session 1: Friday 21, January 2022

Activity 1: Simulations

Simulate a Normal Distribution:

$$x\sim\mathcal{N}(\mu,\sigma^{2})$$

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

plt.rcParams['figure.figsize'] = (7, 8)

rng = np.random.default_rng(1234)   # Set seed for reproducibility

N = 500                              # Number of observations
sigma = np.array([0.1, 1, 2])        # Standard deviation sigma
mu    = np.array([0, 0, 0])          # Mean mu

# This loop simulates the random variable
x_vec = np.empty((N, 3))
for i in range(3):
    x_vec[:, i] = mu[i] + sigma[i] * rng.standard_normal(N)

# Generate a density histogram
num_bins_grid = [100, 50, 10]        # Number of bins to try (columns)

fig, axes = plt.subplots(3, 3, figsize=(11, 8))
for i in range(3):                       # rows: sigma / mu combination
    for j, num_bins in enumerate(num_bins_grid):   # columns: bin count
        ax = axes[i, j]
        ax.hist(x_vec[:, i], bins=num_bins, density=True, color='steelblue', edgecolor='white')
        ax.set_ylabel('pdf')
        ax.set_title(rf'Bins $=\,${num_bins} $\mu=\,${mu[i]} $\sigma^{{2}}=\,${round(sigma[i]**2, 6)}')
plt.tight_layout()
plt.show()

Activity 2: Kernel Density (Estimation)

Histogram:

$$f^{hist}(x)= \frac{1}{N}\sum^{N}_{i=1}\frac{1}{b} \textbf{1}_{x_{i}\in[x-0.5b,x+0.5b]}$$

Where $b$ is the bin width.

Kernel:

$$f^{Kernel}(x)= \frac{1}{N\hat{b}}\sum^{N}_{i=1}K\left(\frac{x_{i}-x}{b}\right)$$

Where $\hat{b}$ is the bandwidth and $K(\cdot)$ is the kernel function. The most common types are:

Normal
Box
Triangle
Epanechnikov

Practical Example

N2 = 300
x2 = rng.normal(loc=mu[1], scale=sigma[1], size=N2)   # Simulate again

support_grid = np.linspace(-4, 4, 100)

# Default-bandwidth kernel density
kde_default = stats.gaussian_kde(x2)
density_default = kde_default(support_grid)
default_bw = kde_default.factor * x2.std(ddof=1)      # effective bandwidth used

# Theoretical (true) density
true_density = stats.norm.pdf(support_grid, loc=mu[1], scale=sigma[1])

fig, ax = plt.subplots(figsize=(7, 4.5))
ax.hist(x2, bins=50, density=True, color='lightgray', edgecolor='white', label='_nolegend_')
ax.plot(support_grid, density_default, 'r-o', markersize=3, linewidth=1.2, label='Kernel Density')
ax.plot(support_grid, true_density, 'b-x', markersize=3, linewidth=1.2, label='Underlying Density')
ax.set_ylabel('pdf')
ax.legend()
plt.tight_layout()
plt.show()

# Now compare a few explicit bandwidths
bws = [default_bw, 1.8, 3]
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.hist(x2, bins=50, density=True, color='lightgray', edgecolor='white', label='_nolegend_')
styles = ['r-', 'r--', 'r-.']
for bw, style in zip(bws, styles):
    kde = stats.gaussian_kde(x2, bw_method=bw / x2.std(ddof=1))
    ax.plot(support_grid, kde(support_grid), style, linewidth=1.2, label=f'Kernel BW = {bw:.4g}')
ax.set_ylabel('pdf')
ax.legend()
plt.tight_layout()
plt.show()

Simulate Data from any distribution

The inverse CDF method

$$x = F^{-1}(u)$$

Where $F$ is the CDF of the desired distribution and $u$ are draws from $U(0,1)$.

Repeat the above example for the Generalized Pareto distribution, using our own inverse-CDF simulator, and check it against the theoretical density.

def sim_gp(n, k, sigma, theta, rng):
    '''
    This function simulates a vector data for the generalized Pareto using the inverse
    CDF method:
    '''
    u = rng.random(n)
    return sigma / k * ((1 - u) ** (-k) - 1) + theta

N3 = 5000
sigma_gp, k_gp, theta_gp = 1.0, 5.0, 0.0

x3 = sim_gp(N3, k_gp, sigma_gp, theta_gp, rng)
x3 = x3[x3 <= 10]           # trim the extreme tail, as in the original
N3 = x3.size

support_gp = np.linspace(0, 10, 100)
true_density_gp = stats.genpareto.pdf(support_gp, c=k_gp, loc=theta_gp, scale=sigma_gp)

kde_default = stats.gaussian_kde(x3)
default_bw3 = kde_default.factor * x3.std(ddof=1)

fig, ax = plt.subplots(figsize=(7, 4.5))
ax.hist(x3, bins=50, density=True, color='lightgray', edgecolor='white', label='_nolegend_')
ax.plot(support_gp, kde_default(support_gp), 'r-o', markersize=3, linewidth=1.2, label='Kernel Density')
ax.plot(support_gp, true_density_gp, 'b-x', markersize=3, linewidth=1.2, label='Underlying Density')
ax.set_ylabel('pdf')
ax.legend()
plt.tight_layout()
plt.show()

bws = [default_bw3, 1.8, 3]
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.hist(x3, bins=50, density=True, color='lightgray', edgecolor='white', label='_nolegend_')
for bw, style in zip(bws, ['r-', 'r--', 'r-.']):
    kde = stats.gaussian_kde(x3, bw_method=bw / x3.std(ddof=1))
    ax.plot(support_gp, kde(support_gp), style, linewidth=1.2, label=f'Kernel BW = {bw:.4g}')
ax.set_ylabel('pdf')
ax.legend()
plt.tight_layout()
plt.show()

Activity 3: Spurious relations

$$x_{1} = \beta_{0}+\beta_{1}x_{2}+\epsilon_{1} \qquad x_{3} = \alpha_{0}+\alpha_{1}x_{2}+\epsilon_{3}$$

rng2 = np.random.default_rng(4567)  # fresh seed for this activity

n = 100
x2_mu, x2_sigma = 3.0, 1.0
eps_sigma1, eps_sigma2 = 0.1, 0.3
beta0, beta1 = 4.0, 0.3
alpha0, alpha1 = -1.0, 2.0

x2_vec = x2_mu + x2_sigma * rng2.standard_normal(n)
eps_1 = eps_sigma1 * rng2.standard_normal(n)
eps_2 = eps_sigma2 * rng2.standard_normal(n)
x1_vec = beta0 + beta1 * x2_vec + eps_1
x3_vec = alpha0 + alpha1 * x2_vec + eps_2

corr_table = pd.DataFrame({
    'rho(x1,x2)': [np.corrcoef(x1_vec, x2_vec)[0, 1]],
    'rho(x3,x2)': [np.corrcoef(x3_vec, x2_vec)[0, 1]],
    'rho(x1,x3)': [np.corrcoef(x1_vec, x3_vec)[0, 1]],
})
corr_table

fig, axes = plt.subplots(1, 3, figsize=(13, 4))
axes[0].plot(x1_vec, x2_vec, 'ko', markerfacecolor='k')
axes[0].set_xlabel('$x_1$', fontsize=14)
axes[0].set_ylabel('$x_2$', fontsize=14)

axes[1].plot(x3_vec, x2_vec, 'ko', markerfacecolor='k')
axes[1].set_xlabel('$x_3$', fontsize=14)
axes[1].set_ylabel('$x_2$', fontsize=14)

axes[2].plot(x1_vec, x3_vec, 'ko', markerfacecolor='k')
axes[2].set_xlabel('$x_1$', fontsize=14)
axes[2].set_ylabel('$x_3$', fontsize=14)

plt.tight_layout()
plt.show()