13. Pooling Cross Sections across Time: Simple Panel Data Methods

# 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 import FirstDifferenceOLS, PanelOLS

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

# Load the CPS data for 1978 and 1985
cps78_85 = wool.data("cps78_85")

# Display basic information
print("Dataset shape:", cps78_85.shape)
print("\nVariable names:", cps78_85.columns.tolist())
print("\nSummary statistics:")
display(cps78_85.describe())

# Check the distribution of observations across years
year_counts = cps78_85["year"].value_counts().sort_index()
print("Observations by year:")
print(year_counts)

# Create year dummy (1985 = 1, 1978 = 0)
cps78_85["y85"] = (cps78_85["year"] == 85).astype(int)

# Estimate model with interaction terms
# The formula y85*(educ+female) expands to:
# y85 + educ + female + y85:educ + y85:female
model_interact = smf.ols(
    formula="lwage ~ y85*(educ+female) + exper + I((exper**2)/100) + union",
    data=cps78_85,
)
results_interact = model_interact.fit()

# Create formatted results table
results_table = pd.DataFrame(
    {
        "Coefficient": results_interact.params,
        "Std. Error": results_interact.bse,
        "t-statistic": results_interact.tvalues,
        "p-value": results_interact.pvalues,
    },
)
print("Pooled Regression with Interaction Terms")
print("=" * 70)
display(results_table.round(4))
print(f"\nR-squared: {results_interact.rsquared:.4f}")
print(f"Observations: {int(results_interact.nobs)}")

# Interpretation of key coefficients
print("INTERPRETATION OF KEY RESULTS:")
print("=" * 70)

# Base return to education in 1978
base_educ = results_interact.params["educ"]
print(f"\n1. Return to Education in 1978: {base_educ:.4f}")
print(f"   A one-year increase in education raises log(wage) by {base_educ:.4f}")
print(f"   Approximately {100 * base_educ:.2f}% increase in wages")

# Change in return to education from 1978 to 1985
change_educ = results_interact.params["y85:educ"]
print(f"\n2. Change in Return to Education (1985 vs 1978): {change_educ:.4f}")
if change_educ > 0:
    print(
        f"   The return to education INCREASED by {100 * change_educ:.2f} percentage points"
    )
else:
    print(
        f"   The return to education DECREASED by {-100 * change_educ:.2f} percentage points"
    )

# Total return to education in 1985
total_educ_1985 = base_educ + change_educ
print(f"\n3. Return to Education in 1985: {total_educ_1985:.4f}")
print(f"   This is {100 * total_educ_1985:.2f}% wage increase per year of education")

# Gender wage gap
female_coef = results_interact.params["female"]
print(f"\n4. Gender Wage Gap in 1978: {female_coef:.4f}")
print(f"   Women earned approximately {100 * female_coef:.2f}% less than men")

# Change in gender wage gap
if "y85:female" in results_interact.params:
    change_female = results_interact.params["y85:female"]
    print(f"\n5. Change in Gender Gap (1985 vs 1978): {change_female:.4f}")
    total_female_1985 = female_coef + change_female
    print(f"   Gender gap in 1985: {100 * total_female_1985:.2f}%")

# Load housing price data from before and after incinerator announcement
kielmc = wool.data("kielmc")

print("Dataset Information:")
print(f"Total observations: {len(kielmc)}")
print(f"\nYears in dataset: {sorted(kielmc['year'].unique())}")
print("\nVariable nearinc: 1 if near incinerator, 0 otherwise")

# Examine the distribution of houses by year and proximity
cross_tab = pd.crosstab(
    kielmc["year"],
    kielmc["nearinc"],
    margins=True,
    margins_name="Total",
)
print("Houses by Year and Proximity to Incinerator:")
display(cross_tab)

# Summary statistics by group
print("\nMean Housing Prices by Group:")
summary_stats = (
    kielmc.groupby(["year", "nearinc"])["rprice"].agg(["count", "mean", "std"]).round(2)
)
display(summary_stats)

# Regression for 1978 (before incinerator announcement)
year_1978 = kielmc["year"] == 1978
model_1978 = smf.ols(formula="rprice ~ nearinc", data=kielmc, subset=year_1978)
results_1978 = model_1978.fit()

# Regression for 1981 (after incinerator announcement)
year_1981 = kielmc["year"] == 1981
model_1981 = smf.ols(formula="rprice ~ nearinc", data=kielmc, subset=year_1981)
results_1981 = model_1981.fit()

# Display results for 1978
print("REGRESSION RESULTS FOR 1978 (Before Incinerator)")
print("=" * 70)
table_1978 = pd.DataFrame(
    {
        "Coefficient": results_1978.params,
        "Std. Error": results_1978.bse,
        "t-statistic": results_1978.tvalues,
        "p-value": results_1978.pvalues,
    },
)
display(table_1978.round(4))
print(f"R-squared: {results_1978.rsquared:.4f}\n")

# Display results for 1981
print("REGRESSION RESULTS FOR 1981 (After Incinerator)")
print("=" * 70)
table_1981 = pd.DataFrame(
    {
        "Coefficient": results_1981.params,
        "Std. Error": results_1981.bse,
        "t-statistic": results_1981.tvalues,
        "p-value": results_1981.pvalues,
    },
)
display(table_1981.round(4))
print(f"R-squared: {results_1981.rsquared:.4f}")

# Calculate simple difference-in-differences manually
print("\nMANUAL DIFFERENCE-IN-DIFFERENCES CALCULATION:")
print("=" * 70)

# Get coefficients
nearinc_1978 = results_1978.params["nearinc"]
nearinc_1981 = results_1981.params["nearinc"]

# Calculate DiD
did_estimate = nearinc_1981 - nearinc_1978

print(f"Effect in 1978 (before): ${nearinc_1978:,.2f}")
print(f"Effect in 1981 (after):  ${nearinc_1981:,.2f}")
print(f"\nDiD Estimate: ${did_estimate:,.2f}")
print(
    f"Interpretation: The incinerator caused a ${abs(did_estimate):,.2f} {'decrease' if did_estimate < 0 else 'increase'} in nearby housing prices",
)

# Pooled regression with interaction term
# This automatically computes the DiD estimator
model_joint = smf.ols(formula="rprice ~ nearinc * C(year)", data=kielmc)
results_joint = model_joint.fit()

print("\nPOOLED REGRESSION WITH INTERACTION (Difference-in-Differences)")
print("=" * 70)
table_joint = pd.DataFrame(
    {
        "Coefficient": results_joint.params,
        "Std. Error": results_joint.bse,
        "t-statistic": results_joint.tvalues,
        "p-value": results_joint.pvalues,
    },
)
display(table_joint.round(4))
print(f"\nR-squared: {results_joint.rsquared:.4f}")
print(f"Observations: {int(results_joint.nobs)}")

# Interpretation of the DiD estimate from pooled regression
print("\nINTERPRETATION OF POOLED REGRESSION:")
print("=" * 70)

# The interaction term is the DiD estimator
interaction_coef = results_joint.params["nearinc:C(year)[T.1981]"]
interaction_se = results_joint.bse["nearinc:C(year)[T.1981]"]
interaction_t = results_joint.tvalues["nearinc:C(year)[T.1981]"]
interaction_p = results_joint.pvalues["nearinc:C(year)[T.1981]"]

print(f"\nDiD Estimate (interaction term): ${interaction_coef:,.2f}")
print(f"Standard Error: ${interaction_se:,.2f}")
print(f"t-statistic: {interaction_t:.4f}")
print(f"p-value: {interaction_p:.4f}")

if interaction_p < 0.05:
    print(
        "\nYES The effect is statistically significant at the 5% level",
    )
    print(
        f"  The incinerator caused a ${abs(interaction_coef):,.2f} {'decrease' if interaction_coef < 0 else 'increase'} in nearby housing prices",
    )
else:
    print(
        "\nNO The effect is NOT statistically significant at the 5% level",
    )
    print("  We cannot conclude the incinerator had a causal effect on prices")

crime2 = wool.data("crime2")

# Wide format: one row per city with year-82 and year-87 columns.
crime_wide = (
    crime2.assign(city=(crime2.index // 2).astype(int))
    .pivot(index="city", columns="year", values=["crmrte", "unem"])
    .dropna()
)
crime_wide.columns = [f"{var}_{yr}" for var, yr in crime_wide.columns]

# First differences
crime_wide["dcrmrte"] = crime_wide["crmrte_87"] - crime_wide["crmrte_82"]
crime_wide["dunem"] = crime_wide["unem_87"] - crime_wide["unem_82"]

# Pooled OLS (ignores the city effect) vs first-difference estimator.
pooled = smf.ols("crmrte ~ unem", data=crime2).fit()
fd = smf.ols("dcrmrte ~ dunem", data=crime_wide).fit()

comparison = pd.DataFrame(
    {
        "Estimator": ["Pooled OLS (crmrte on unem)", "First differences (Δcrmrte on Δunem)"],
        "β̂_unem": [pooled.params["unem"], fd.params["dunem"]],
        "SE": [pooled.bse["unem"], fd.bse["dunem"]],
        "t-stat": [pooled.tvalues["unem"], fd.tvalues["dunem"]],
        "R²": [pooled.rsquared, fd.rsquared],
    },
)
display(comparison.round(3))

kielmc = wool.data("kielmc")

# Standard interaction specification on the level rprice.
did_levels = smf.ols("rprice ~ y81 * nearinc", data=kielmc).fit()

# Same in logs, which removes much of the heteroskedasticity in housing prices.
did_logs = smf.ols("np.log(rprice) ~ y81 * nearinc", data=kielmc).fit()

# Report the four sample means that determine the DiD estimate.
cell_means = (
    kielmc.assign(group=np.where(kielmc["nearinc"] == 1, "near", "far"),
                  period=np.where(kielmc["y81"] == 1, "1981", "1978"))
    .groupby(["group", "period"], observed=True)["rprice"]
    .mean()
    .unstack()
)
print("Sample means of real price (1978$):")
display(cell_means.round(0))
manual_did = (cell_means.loc["near", "1981"] - cell_means.loc["near", "1978"]) - (
    cell_means.loc["far", "1981"] - cell_means.loc["far", "1978"]
)
print(f"\nManual DiD (delta near minus delta far) = {manual_did:,.0f} dollars")
print(f"Regression DiD (y81:nearinc coefficient) = {did_levels.params['y81:nearinc']:,.0f} dollars")
print(f"Regression DiD in logs (≈ percentage change) = {did_logs.params['y81:nearinc']:.3f}")

# Load crime panel data
crime4 = wool.data("crime4")

print("Dataset Information:")
print(f"Total observations: {len(crime4)}")
print(f"Number of counties: {crime4['county'].nunique()}")
print(f"Years: {sorted(crime4['year'].unique())}")
print(f"Time periods: {len(crime4['year'].unique())}")

# Display first few observations
print("\nFirst few observations:")
display(crime4.head(10))

# Summary statistics for key variables
key_vars = ["crmrte", "prbarr", "prbconv", "prbpris", "avgsen", "polpc"]
print("\nSummary Statistics for Key Variables:")
display(crime4[key_vars].describe().round(4))

# Set the panel structure (county and year as index)
crime4_panel = crime4.set_index(["county", "year"], drop=False)

# Create year dummies manually
for year in [83, 84, 85, 86, 87]:
    crime4_panel[f"d{year}"] = (crime4_panel["year"] == year).astype(int)

print("Panel data structure set up successfully")
print(f"Panel index levels: {crime4_panel.index.names}")

# Estimate first-difference model using linearmodels
# This package provides specialized panel data estimators

fd_model = FirstDifferenceOLS.from_formula(
    formula="np.log(crmrte) ~ year + d83 + d84 + d85 + d86 + d87 + "
    "np.log(prbarr) + np.log(prbconv) + np.log(prbpris) + "
    "np.log(avgsen) + np.log(polpc)",
    data=crime4_panel,
)

fd_results = fd_model.fit()

print("FIRST-DIFFERENCE ESTIMATION RESULTS")
print("=" * 70)
print(fd_results)

# Interpretation of results
print("\nINTERPRETATION OF FIRST-DIFFERENCE ESTIMATES:")
print("=" * 70)

# Extract coefficients
params = fd_results.params

print("\nElasticity Interpretations (all variables in logs):")
print(
    f"\n1. log(prbarr) - Probability of Arrest: {params['np.log(prbarr)']:.4f}",
)
print(
    "   A 1% increase in arrest probability reduces crime by "
    f"{abs(params['np.log(prbarr)']):.2f}%"
)

print(
    f"\n2. log(prbconv) - Probability of Conviction: {params['np.log(prbconv)']:.4f}",
)
print(
    "   A 1% increase in conviction probability affects crime by "
    f"{params['np.log(prbconv)']:.2f}%"
)

print(
    f"\n3. log(prbpris) - Probability of Prison: {params['np.log(prbpris)']:.4f}",
)
print(
    "   A 1% increase in prison probability affects crime by "
    f"{params['np.log(prbpris)']:.2f}%"
)

print(f"\n4. log(avgsen) - Average Sentence: {params['np.log(avgsen)']:.4f}")
print(
    "   A 1% increase in sentence length affects crime by "
    f"{params['np.log(avgsen)']:.2f}%"
)

print(f"\n5. log(polpc) - Police Per Capita: {params['np.log(polpc)']:.4f}")
print(
    "   A 1% increase in police per capita affects crime by "
    f"{params['np.log(polpc)']:.2f}%"
)

print("\n" + "=" * 70)
print("Note: First-differencing controls for all time-invariant county")
print("characteristics (e.g., poverty, culture, geography) that might")
print("be correlated with both crime and law enforcement variables.")

# Compare FD and FE estimators
fe_model = PanelOLS.from_formula(
    formula="np.log(crmrte) ~ year + d83 + d84 + d85 + d86 + d87 + "
    "np.log(prbarr) + np.log(prbconv) + np.log(prbpris) + "
    "np.log(avgsen) + np.log(polpc) + EntityEffects",
    data=crime4_panel,
)

fe_results = fe_model.fit()

print("FIXED EFFECTS ESTIMATION RESULTS")
print("=" * 70)
print(fe_results)

# Compare FD and FE coefficients
comparison = pd.DataFrame(
    {
        "First-Difference": fd_results.params,
        "Fixed Effects": fe_results.params,
        "Difference": fd_results.params - fe_results.params,
    },
)

print("\nCOMPARISON OF FD AND FE ESTIMATES:")
print("=" * 70)
display(comparison.round(4))

print("\nNote: With many time periods (T=7 here), FD and FE can give")
print("different estimates. The choice depends on the structure of")
print("serial correlation in the errors.")