16. Simultaneous Equations Models

# Import required libraries
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import wooldridge as wool
from IPython.display import display
from linearmodels.iv import IV2SLS
from linearmodels.system import IV3SLS

plt.style.use("theme.mplstyle")

rng = np.random.default_rng(2026)
n = 800

# Exogenous demand shifter (income) and supply shifter (weather).
income = rng.normal(0, 1, size=n)
weather = rng.normal(0, 1, size=n)
u_d = rng.normal(0, 1, size=n)
u_s = rng.normal(0, 1, size=n)

# True structural parameters
alpha = {"intercept": 5.0, "price": -1.5, "income": 1.0}  # demand: Q = α0 + α1 P + α2 I + u_d
beta = {"intercept": 1.0, "price": 0.8, "weather": 1.2}   # supply: Q = β0 + β1 P + β2 W + u_s

# Reduced form derived from setting D = S and solving for price.
denom = beta["price"] - alpha["price"]
price = (
    (alpha["intercept"] - beta["intercept"])
    + alpha["income"] * income
    - beta["weather"] * weather
    + (u_d - u_s)
) / denom
quantity = alpha["intercept"] + alpha["price"] * price + alpha["income"] * income + u_d

sim = pd.DataFrame({"price": price, "quantity": quantity, "income": income, "weather": weather})

# Structural OLS (biased: price is endogenous in both equations).
demand_ols = smf.ols("quantity ~ price + income", data=sim).fit()
supply_ols = smf.ols("quantity ~ price + weather", data=sim).fit()

# Reduced form: regress each endogenous variable on ALL exogenous variables.
rf_price = smf.ols("price ~ income + weather", data=sim).fit()
rf_quantity = smf.ols("quantity ~ income + weather", data=sim).fit()

summary = pd.DataFrame(
    {
        "True α₁ (price in demand)": [alpha["price"]],
        "OLS demand β̂_price": [demand_ols.params["price"]],
        "True β₁ (price in supply)": [beta["price"]],
        "OLS supply β̂_price": [supply_ols.params["price"]],
    },
).T.rename(columns={0: "value"})
display(summary.round(3))

print("Reduced-form regressions (consistent: all RHS exogenous):")
display(
    pd.DataFrame(
        {
            "price ~": [rf_price.params["income"], rf_price.params["weather"]],
            "quantity ~": [rf_quantity.params["income"], rf_quantity.params["weather"]],
        },
        index=["income", "weather"],
    ).round(3),
)

# Simulate supply and demand with simultaneity
np.random.seed(42)
n = 200

# Exogenous variables
income = np.random.normal(50, 10, n)  # Demand shifter
weather = np.random.normal(0, 5, n)  # Supply shifter

# Structural errors
u_demand = np.random.normal(0, 3, n)
u_supply = np.random.normal(0, 3, n)

# True structural parameters
alpha_0, alpha_1, alpha_2 = 100, 2, 0.5  # Demand: Q = 100 - 2P + 0.5I + u_D
beta_0, beta_1, beta_2 = 50, 1.5, 0.3  # Supply: Q = 50 + 1.5P + 0.3W + u_S

# Solve for equilibrium price and quantity
# Q_D = Q_S -> alpha_0 - alpha_1*P + alpha_2*I + u_D = beta_0 + beta_1*P + beta_2*W + u_S
# P = (alpha_0 - beta_0 + alpha_2*I - beta_2*W + u_D - u_S) / (alpha_1 + beta_1)

price_equilibrium = (
    alpha_0 - beta_0 + alpha_2 * income - beta_2 * weather + u_demand - u_supply
) / (alpha_1 + beta_1)
quantity_equilibrium = (
    alpha_0 - alpha_1 * price_equilibrium + alpha_2 * income + u_demand
)

# Create DataFrame
sim_data = pd.DataFrame(
    {
        "quantity": quantity_equilibrium,
        "price": price_equilibrium,
        "income": income,
        "weather": weather,
    },
)

print("SIMULATED SUPPLY AND DEMAND DATA")
print("=" * 70)
print("\nTrue Structural Parameters:")
print(f"Demand:  Q = {alpha_0} - {alpha_1}*P + {alpha_2}*I + u_D")
print(f"Supply:  Q = {beta_0} + {beta_1}*P + {beta_2}*W + u_S")
print("\nObserved equilibrium data (Q*, P*) where supply = demand")

display(sim_data.describe().round(2))

mroz = wool.data("mroz")
labor = mroz.dropna(subset=["lwage"]).copy()
labor["expersq"] = labor["exper"] ** 2

# OLS treats log(wage) as exogenous (biased under simultaneity).
ols_labor = smf.ols("hours ~ lwage + educ + age + kidslt6 + nwifeinc", data=labor).fit()

# 2SLS: instrument log(wage) with exper and exper^2.
iv_labor = IV2SLS.from_formula(
    "hours ~ 1 + educ + age + kidslt6 + nwifeinc + [lwage ~ exper + expersq]",
    data=labor,
).fit(cov_type="robust")

# First stage to confirm relevance.
first_stage = smf.ols(
    "lwage ~ exper + expersq + educ + age + kidslt6 + nwifeinc", data=labor,
).fit()
joint_F = first_stage.f_test("exper = 0, expersq = 0").fvalue

comparison = pd.DataFrame(
    {
        "OLS β̂": [ols_labor.params[k] for k in ["lwage", "educ", "kidslt6"]],
        "OLS SE": [ols_labor.bse[k] for k in ["lwage", "educ", "kidslt6"]],
        "2SLS β̂": [iv_labor.params[k] for k in ["lwage", "educ", "kidslt6"]],
        "2SLS SE": [iv_labor.std_errors[k] for k in ["lwage", "educ", "kidslt6"]],
    },
    index=["lwage (endogenous)", "educ", "kidslt6"],
)
display(comparison.round(2))
print(f"\nFirst-stage joint F on excluded instruments (exper, exper²): {float(joint_F):.2f}")
print(f"Sargan overid p-value (H₀: instruments valid): {iv_labor.sargan.pval:.3f}")

# Load married women's labor force participation data
mroz = wool.data("mroz")

# Restrict to working women (non-missing wages)
mroz_working = mroz.dropna(subset=["lwage"])

print("SIMULTANEOUS EQUATIONS: Labor Supply and Wages")
print("=" * 70)
print("\nSample: Married women with observed wages")
print(f"Observations: {len(mroz_working)}")

print("\nSTRUCTURAL EQUATIONS:")
print("\n1. Labor Supply Equation:")
print(
    "   hours = beta_1_0 + beta_1_2*log(wage) + beta_1_1*educ + beta_1_2*age + beta_1_3*kidslt6 + beta_1_4*nwifeinc + u_1"
)
print("   Endogenous: log(wage)")
print("   Instruments: exper, exper^2")

print("\n2. Wage Equation:")
print("   log(wage) = beta_2_0 + beta_2_1*hours + beta_2_2*educ + beta_2_3*exper + beta_2_4*exper^2 + u_2")
print("   Endogenous: hours")
print("   Instruments: age, kidslt6, nwifeinc")

print("\nKEY QUESTION: Are hours and wages simultaneously determined?")
print("- Labor supply theory: Higher wages -> More hours worked (substitution effect)")
print("- Wage determination: More hours -> Lower hourly wage? (fatigue, bargaining)")

# Summary statistics
key_vars = ["hours", "lwage", "educ", "age", "kidslt6", "nwifeinc", "exper"]
display(mroz_working[key_vars].describe().round(2))

# EQUATION-BY-EQUATION 2SLS ESTIMATION
print("\n2SLS ESTIMATION (Equation-by-Equation)")
print("=" * 70)

# Equation 1: Labor Supply (hours depends on wage)
print("\nEQUATION 1: Labor Supply")
print("-" * 70)
print("Dependent variable: hours")
print("Endogenous regressor: log(wage)")
print("Instruments: exper, exper^2")

labor_supply_2sls = IV2SLS.from_formula(
    formula="hours ~ 1 + educ + age + kidslt6 + nwifeinc + "
    "[np.log(wage) ~ exper + I(exper**2)]",
    data=mroz_working,
).fit(cov_type="unadjusted", debiased=True)

table_supply = pd.DataFrame(
    {
        "Coefficient": labor_supply_2sls.params,
        "Std. Error": labor_supply_2sls.std_errors,
        "t-statistic": labor_supply_2sls.tstats,
        "p-value": labor_supply_2sls.pvalues,
    },
)
display(table_supply.round(4))

print(f"\nR-squared: {labor_supply_2sls.rsquared:.4f}")
print("First-stage F-stat: Check that instruments are relevant")

# Equation 2: Wage Equation (wage depends on hours)
print("\nEQUATION 2: Wage Equation")
print("-" * 70)
print("Dependent variable: log(wage)")
print("Endogenous regressor: hours")
print("Instruments: age, kidslt6, nwifeinc")

wage_eq_2sls = IV2SLS.from_formula(
    formula="np.log(wage) ~ 1 + educ + exper + I(exper**2) + "
    "[hours ~ age + kidslt6 + nwifeinc]",
    data=mroz_working,
).fit(cov_type="unadjusted", debiased=True)

table_wage = pd.DataFrame(
    {
        "Coefficient": wage_eq_2sls.params,
        "Std. Error": wage_eq_2sls.std_errors,
        "t-statistic": wage_eq_2sls.tstats,
        "p-value": wage_eq_2sls.pvalues,
    },
)
display(table_wage.round(4))

print(f"\nR-squared: {wage_eq_2sls.rsquared:.4f}")

# Check for cross-equation error correlation
print("\nCHECKING CROSS-EQUATION ERROR CORRELATION")
print("=" * 70)
print("If residuals are correlated, 3SLS will be more efficient than 2SLS")

# Get residuals from both equations
u1_hat = labor_supply_2sls.resids.values
u2_hat = wage_eq_2sls.resids.values

# Calculate correlation
cor_u1u2 = np.corrcoef(u1_hat, u2_hat)[0, 1]

print(f"\nCorrelation between residuals: {cor_u1u2:.4f}")

if abs(cor_u1u2) > 0.1:
    print(f"YES Substantial correlation (|{cor_u1u2:.4f}| > 0.1)")
    print("  -> 3SLS will be more efficient than 2SLS")
    print("  -> Errors across equations are related")
else:
    print(f"NO Low correlation (|{cor_u1u2:.4f}| <= 0.1)")
    print("  -> Little efficiency gain from 3SLS")
    print("  -> 2SLS and 3SLS will give similar results")

# Interpretation of 2SLS results
print("\nINTERPRETATION OF 2SLS RESULTS")
print("=" * 70)

# Labor supply equation
wage_effect = labor_supply_2sls.params["np.log(wage)"]
wage_se = labor_supply_2sls.std_errors["np.log(wage)"]

print("\n1. LABOR SUPPLY EQUATION:")
print(f"   Effect of wage on hours: {wage_effect:.2f}")
if wage_effect > 0:
    print(f"   A 1% increase in wage -> {wage_effect:.2f} more hours worked per year")
    print("   YES Positive substitution effect (higher wage -> more labor supply)")
else:
    print(f"   A 1% increase in wage -> {abs(wage_effect):.2f} fewer hours worked")
    print("   YES Negative income effect dominates (higher wage -> less need to work)")

print(f"\n   Standard error: {wage_se:.4f}")
if abs(wage_effect / wage_se) > 1.96:
    print(f"   YES Statistically significant (t = {wage_effect / wage_se:.2f})")
else:
    print(f"   NO Not statistically significant (t = {wage_effect / wage_se:.2f})")

# Wage equation
hours_effect = wage_eq_2sls.params["hours"]
hours_se = wage_eq_2sls.std_errors["hours"]

print("\n2. WAGE EQUATION:")
print(f"   Effect of hours on log(wage): {hours_effect:.6f}")
if hours_effect > 0:
    print(f"   Each additional hour worked -> {100 * hours_effect:.4f}% higher wage")
    print("   Interpretation: Positive correlation (more productive workers work more)")
else:
    print(f"   Each additional hour worked -> {100 * abs(hours_effect):.4f}% lower wage")
    print("   Interpretation: Fatigue effect or part-time premium")

print(f"\n   Standard error: {hours_se:.6f}")
if abs(hours_effect / hours_se) > 1.96:
    print(f"   YES Statistically significant (t = {hours_effect / hours_se:.2f})")
else:
    print(f"   NO Not statistically significant (t = {hours_effect / hours_se:.2f})")

# THREE-STAGE LEAST SQUARES (3SLS)
print("\n\n3SLS ESTIMATION (System Estimation)")
print("=" * 70)
print("Estimates entire system jointly, accounting for cross-equation correlations")

# Define system of equations
system_formula = {
    "labor_supply": "hours ~ 1 + educ + age + kidslt6 + nwifeinc + "
    "[np.log(wage) ~ exper + I(exper**2)]",
    "wage_eq": "np.log(wage) ~ 1 + educ + exper + I(exper**2) + "
    "[hours ~ age + kidslt6 + nwifeinc]",
}

# Estimate system
system_3sls = IV3SLS.from_formula(system_formula, data=mroz_working)
results_3sls = system_3sls.fit(cov_type="unadjusted", debiased=True)

print("\n3SLS SYSTEM RESULTS")
print(results_3sls)

# Compare 2SLS and 3SLS estimates
print("\n\nCOMPARISON: 2SLS vs 3SLS")
print("=" * 70)

print("\nLabor Supply Equation (hours):")
print("-" * 70)

comparison_supply = pd.DataFrame(
    {
        "2SLS": labor_supply_2sls.params,
        "2SLS SE": labor_supply_2sls.std_errors,
        "3SLS": results_3sls.params.filter(like="labor_supply"),
        "3SLS SE": results_3sls.std_errors.filter(like="labor_supply"),
    },
)

# Clean up index names for 3SLS
comparison_supply.index = [
    idx.replace("labor_supply_", "") for idx in comparison_supply.index
]
display(comparison_supply.round(4))

print("\nWage Equation (log(wage)):")
print("-" * 70)

comparison_wage = pd.DataFrame(
    {
        "2SLS": wage_eq_2sls.params,
        "2SLS SE": wage_eq_2sls.std_errors,
        "3SLS": results_3sls.params.filter(like="wage_eq"),
        "3SLS SE": results_3sls.std_errors.filter(like="wage_eq"),
    },
)

# Clean up index names for 3SLS
comparison_wage.index = [idx.replace("wage_eq_", "") for idx in comparison_wage.index]
display(comparison_wage.round(4))

print("\nKEY INSIGHTS:")
print("1. Coefficients: 2SLS and 3SLS estimates are similar (as expected)")
print("2. Standard Errors: 3SLS should have smaller SEs if errors are correlated")
print("3. Both methods account for simultaneity (unlike OLS which is biased)")
print("\nWHICH TO USE?")
print("- 3SLS: More efficient if system is correctly specified")
print("- 2SLS: More robust to misspecification of individual equations")
print("- In practice: Report both and compare for robustness")

# Conceptual example of macro SEM
print("TIME SERIES SEM: Macroeconomic Model")
print("=" * 70)

print("\nIS-LM MODEL (Conceptual)")
print("\nIS Equation (Goods Market Equilibrium):")
print("  Y = alpha_0 + alpha_1*r + alpha_2*G + u")
print("  Output (Y) depends negatively on interest rate (r)")
print("  Government spending (G) is exogenous policy variable")

print("\nLM Equation (Money Market Equilibrium):")
print("  r = beta_0 + beta_1*Y + beta_2*M + v")
print("  Interest rate (r) depends positively on output (Y)")
print("  Money supply (M) is exogenous policy variable")

print("\nIDENTIFICATION:")
print("YES IS equation: Use M (excluded from IS) as instrument for r")
print("YES LM equation: Use G (excluded from LM) as instrument for Y")

print("\nESTIMATION:")
print("- Apply 2SLS to each equation using appropriate instruments")
print("- Or use 3SLS for system estimation")
print("- With time series: Check stationarity, use robust SEs")

print("\nEXTENSIONS:")
print("- Add dynamics: Y_{t-1}, r_{t-1} as regressors")
print("- Use lags as instruments (if predetermined)")
print("- Model expectations: Forward-looking behavior")

# Conceptual panel SEM example
print("PANEL DATA SEM: Firm Productivity and Training")
print("=" * 70)

print("\nSTRUCTURAL EQUATIONS:")

print("\n1. Productivity Equation:")
print("   log(output)_i_t = beta_1*training_i_t + controls + alpha_i + u_1_i_t")
print("   Productivity depends on training")

print("\n2. Training Decision:")
print("   training_i_t = beta_2*log(output)_i_t_-_1 + controls + c_i + u_2_i_t")
print("   Training depends on past productivity")

print("\nSIMULTANEITY:")
print("- Firms choose training based on productivity")
print("- Training affects future productivity")
print("- Need to account for both directions of causality")

print("\nIDENTIFICATION STRATEGY:")
print("YES Use lagged productivity as instrument for current training")
print("YES Use policy variables (grants) as instrument for training")
print("YES Fixed effects control for time-invariant firm heterogeneity")

print("\nESTIMATION:")
print("- First-difference to remove fixed effects")
print("- Apply 2SLS using lagged variables as instruments")
print("- Or use system GMM for efficiency")

print("\nADVANTAGES OF PANEL SEM:")
print("YES More instruments available (lags, differences)")
print("YES Control for unobserved heterogeneity (fixed effects)")
print("YES Richer identification from time variation")
print("YES Can test Granger causality")