Python: Estimating Effects

A quick reference for the methods used in class, building on PP5000. We will add methods as we introduce them. See the weekly notes for the reasoning behind the estimates.

Week 1: Differences in means and regression

These syntax templates use a data frame called df, an outcome Y, and a treatment indicator D coded 1 for treatment and 0 for control. Replace these names with the variables in your dataset.

import pandas as pd
import statsmodels.formula.api as smf

Task	Python	What to read
Group means	`sample.groupby("D")["Y"].mean()`	Mean outcome in each group
Difference in means	`means.loc[1] - means.loc[0]`	Treatment mean minus control mean
Regression	`smf.ols("Y ~ D", data=sample)`	Fit the model below; read the coefficient on `D`
Region indicators	`C(region)` in the formula	Adjust for categorical region controls

Difference in means

Use observations with both treatment and outcome observed, so the means and regression refer to the same sample.

sample = df.dropna(subset=["Y", "D"])
means = sample.groupby("D")["Y"].mean()
difference = means.loc[1] - means.loc[0]
print(means)
print(difference)

The difference is in the units of Y. With random assignment and complete outcome data, it estimates the effect of assignment to treatment. Without a credible research design, a difference in means need not be a causal effect.

Regression and standard errors

Use heteroskedasticity-consistent standard errors

For the regressions in Week 1, use heteroskedasticity-consistent standard errors, also called robust standard errors. Specify .fit(cov_type="HC1", use_t=True) when fitting the model. Calling .fit() alone uses the conventional standard errors that assume homoskedastic disturbances. Later exercises may require clustered standard errors instead; follow the instructions for the research design.

The complete specification for a regression of the outcome on an intercept and the treatment indicator is:

import statsmodels.formula.api as smf

sample = df.dropna(subset=["Y", "D"])
model = smf.ols("Y ~ D", data=sample).fit(
    cov_type="HC1", use_t=True
)
print(model.summary())

The formula includes an intercept by default. Its estimate is the control-group mean; the coefficient on D equals the treatment-group mean minus the control-group mean.

cov_type="HC1" requests heteroskedasticity-consistent standard errors. It changes the estimated uncertainty, not the estimated coefficients. use_t=True requests t-based inference, as in the Week 1 exercise.

estimate = model.params["D"]
standard_error = model.bse["D"]
confidence_interval = model.conf_int().loc["D"]
n_observations = model.nobs

Report the estimate, its units, its standard error or confidence interval, and the number of observations. HC1 allows unequal disturbance variances; it does not account for dependence within groups. We will revisit standard errors for other research designs.

Region indicators

Here is the complete specification with region fixed effects and heteroskedasticity-consistent standard errors:

import statsmodels.formula.api as smf

sample_region = df.dropna(subset=["Y", "D", "region"])
model_region = smf.ols(
    formula="Y ~ D + C(region)",
    data=sample_region
).fit(cov_type="HC1", use_t=True)
print(model_region.summary())

C(region) treats region codes as categories and includes indicators, omitting a reference category. Writing region alone treats the codes as a quantitative variable. The coefficient on D is now the regression-adjusted difference; it need not equal the raw difference in means. Missing values in added controls can also change the estimation sample.

Binary outcomes: the linear probability model

When Y is coded 0 or 1, this OLS regression is a linear probability model. A coefficient of 0.04 means a difference of 4 percentage points, not a 4 per cent increase.

If the conditional probability is $p$, the conditional variance of the binary outcome—and hence of its disturbance around the conditional mean—is $p\left(1-p\right)$. This generally varies across observations, motivating heteroskedasticity-consistent standard errors.

To express an estimate in percentage points, multiply the estimate in probability units by 100. Scale its standard error and confidence-interval endpoints by the same factor.

Week 2: Propensity scores and weighting

Use an outcome Y, a treatment indicator D coded 1 for treatment and 0 for control, and pre-treatment characteristics X1 and X2. Replace these names with your variables. The weights below estimate the average treatment effect on the treated (ATT) under the assumptions discussed in the Week 2 notes.

Task	Python	What to read
Estimate a logit	`smf.logit("D ~ X1 + X2", data=sample_ipw).fit()`	Model the probability of treatment given the characteristics
Predict propensity scores	`score_model.predict(sample_ipw)`	Estimated treatment probability for each observation
ATT weights	`1` for treated; `pscore / (1 - pscore)` for controls	Give more weight to controls whose characteristics are common among treated observations
Weighted regression	`smf.wls("Y ~ D", data=sample_ipw, weights=sample_ipw["att_weight"])`	Fit the model below; read the coefficient on `D`

Estimate the propensity score

Keep a common sample for the propensity-score model and outcome regression.

import statsmodels.formula.api as smf

sample_ipw = df.dropna(subset=["Y", "D", "X1", "X2"]).copy()

score_model = smf.logit(
    "D ~ X1 + X2", data=sample_ipw
).fit()

sample_ipw["pscore"] = score_model.predict(sample_ipw)

smf.logit estimates the logit model. .predict() returns the estimated probability of treatment for each observation: its propensity score. For this step, we use those probabilities to construct the weights.

Construct ATT weights

Treated observations receive weight 1. Controls receive $\hat p_i / \left(1-\hat p_i\right)$.

# Start with a weight of 1 for every observation.
sample_ipw["att_weight"] = 1.0

# Replace the weights for the control observations.
controls = sample_ipw["D"] == 0
sample_ipw.loc[controls, "att_weight"] = (
    sample_ipw.loc[controls, "pscore"]
    / (1 - sample_ipw.loc[controls, "pscore"])
)

controls identifies the rows where D equals zero. .loc[controls, ...] selects those rows. Check overlap in propensity scores and covariate balance, and inspect large weights, before interpreting the estimate.

Estimate the weighted regression

The complete outcome regression includes an intercept and the treatment indicator:

weighted = smf.wls(
    "Y ~ D",
    data=sample_ipw,
    weights=sample_ipw["att_weight"],
).fit(cov_type="HC1", use_t=True)

estimate = weighted.params["D"]
standard_error = weighted.bse["D"]
confidence_interval = weighted.conf_int().loc["D"]
n_observations = weighted.nobs

The coefficient on D is the treated-group mean minus the weighted control-group mean, in the units of Y. cov_type="HC1" requests heteroskedasticity-consistent standard errors. These standard errors treat the estimated weights as fixed; a full calculation of uncertainty also accounts for estimating the propensity scores.

For Problem Set 1, use re78 as the outcome and treat as the treatment indicator. Estimate a separate propensity-score model for each observational comparison sample, using the characteristics specified in the assignment.

Week 3: Instrumental variables

Use linearmodels for two-stage least squares. With outcome Y, endogenous treatment D, instrument Z, and a categorical control group:

import linearmodels.iv as iv

sample_iv = df.dropna(subset=["Y", "D", "Z", "group"])
model_iv = iv.IV2SLS.from_formula(
    "Y ~ 1 + C(group) + [D ~ Z]", data=sample_iv
)
result_iv = model_iv.fit(cov_type="robust", debiased=True)

estimate = result_iv.params["D"]
standard_error = result_iv.std_errors["D"]
confidence_interval = result_iv.conf_int().loc["D"]

[D ~ Z] tells the estimator to instrument D with Z. C(group) includes category indicators in both stages. cov_type="robust" requests heteroskedasticity-consistent IV standard errors. In this package, the standard-error attribute is .std_errors.

The Week 3 exercise uses the Oregon survey weights and clusters by household:

model_oregon = iv.IV2SLS.from_formula(
    "outpatient ~ 1 + C(stratum) + [medicaid ~ selected]",
    data=analysis, weights=analysis["weight_12m"],
)
result_oregon = model_oregon.fit(
    cov_type="clustered",
    clusters=analysis["household_id"],
    debiased=True,
)

Here stratum identifies the observed combinations of household-list size and survey wave. The setup file constructs it and selects the common analysis sample. Household clustering permits heteroskedasticity and correlation within households.

First stage, reduced form and the IV ratio

With one excluded instrument, estimate the first stage and reduced form using the same observations, controls and weights as IV. For the Oregon example:

first_oregon = smf.wls(
    "medicaid ~ selected + C(stratum)",
    data=analysis, weights=analysis["weight_12m"],
).fit(
    cov_type="cluster",
    cov_kwds={"groups": analysis["household_id"]},
    use_t=True,
)

reduced_oregon = smf.wls(
    "outpatient ~ selected + C(stratum)",
    data=analysis, weights=analysis["weight_12m"],
).fit(
    cov_type="cluster",
    cov_kwds={"groups": analysis["household_id"]},
    use_t=True,
)

first_stage_effect = first_oregon.params["selected"]
itt_effect = reduced_oregon.params["selected"]
iv_ratio = itt_effect / first_stage_effect
print(iv_ratio)
print(result_oregon.params["medicaid"])

The first stage measures the effect of lottery selection on Medicaid coverage. The reduced form measures its effect on the outcome, also called the intention-to-treat effect here. Their coefficient ratio equals the IV coefficient in this exactly identified model. Obtain its standard error from the fitted IV model.

Display the first stage from an IV model

print(result_iv.first_stage)
print(result_iv.first_stage.diagnostics)

The first command displays the first-stage coefficients and diagnostics. The second gives a data frame of diagnostics, including partial $R^2$ and a test of the excluded instruments. Read f.dist alongside f.stat: with robust or clustered covariance, linearmodels reports a chi-squared Wald statistic for this test. The examples below show how to obtain the conventional first-stage F used in PS2.

Week 4: Leniency instruments and simulations

Construct leave-one-out leniency

In the Week 4 simulation, caseworker identifies the decision-maker and training is 1 for admission and 0 otherwise. The data contain one observation per applicant and no missing admission decisions.

Task	Python	Result
Count admitted applicants by caseworker	`data.groupby("caseworker")["training"].transform("sum")`	Each applicant receives their caseworker’s total admissions
Count applicants by caseworker	`data.groupby("caseworker")["training"].transform("count")`	Each applicant receives their caseworker’s caseload
Leave out the applicant	`(admissions - data["training"]) / (caseload - 1)`	Admission rate among the caseworker’s other applicants

admissions = data.groupby("caseworker")["training"].transform("sum")
caseload = data.groupby("caseworker")["training"].transform("count")

data["leniency"] = (
    (admissions - data["training"]) / (caseload - 1)
)

transform() returns a value for each original row, making it convenient for constructing a new variable. Subtracting the applicant’s own admission decision removes its direct contribution to their instrument. This calculation requires at least two applicants per caseworker.

Estimate IV with leniency

The simulation assigns applicants across caseworkers without additional assignment groups:

import linearmodels.iv as iv

leniency_iv = iv.IV2SLS.from_formula(
    "earnings_after ~ 1 + [training ~ leniency]",
    data=data,
).fit(cov_type="robust", debiased=True)
print(leniency_iv.summary)

In an empirical application, assignment may be comparable only within an office, court or time period. If assignment_group identifies those groups, include their indicators in both stages:

leniency_group_iv = iv.IV2SLS.from_formula(
    "Y ~ 1 + C(assignment_group) + [D ~ leniency]",
    data=assignment_sample,
).fit(cov_type="robust", debiased=True)

Here assignment_sample is the application’s common analysis sample, with outcome Y, treatment D and its constructed leniency instrument. The definition of leniency and the choice of clustered standard errors should follow the assignment process in that application.

Keep a simulated sample reproducible

import numpy as np

rng = np.random.default_rng(12345)

Use this generator for the random draws in the simulation. Re-running the same code with the same seed reproduces the sample. Choose your own integer. np.random.default_rng() with empty parentheses produces a new sample each time.

Week 5: Multiple instruments and first-stage evidence

The examples below use the AK variable names from Problem Set 2. The setup file constructs q1, q2, q3, age_years and age_squared. Quarter 4 is the omitted category.

Estimate IV with several instruments

import statsmodels.formula.api as smf
import linearmodels.iv as iv

controls = "black + married + smsa + C(division) + age_years + age_squared"

ak_iv = iv.IV2SLS.from_formula(
    "lnwkwage ~ 1 + " + controls + " + [educ ~ q1 + q2 + q3]",
    data=data,
).fit(cov_type="unadjusted", debiased=True)
print(ak_iv.summary)

The expression [educ ~ q1 + q2 + q3] specifies schooling as endogenous and the three quarter indicators as excluded instruments. Included controls enter both stages.

Conventional standard errors for replication

PS2 uses conventional standard errors to compare with AK and BJB. These assume homoskedastic disturbances.

Covariance choice	`statsmodels` OLS	`linearmodels` IV
Conventional	`.fit()`	`.fit(cov_type="unadjusted", debiased=True)`
Heteroskedasticity-consistent	`.fit(cov_type="HC1", use_t=True)`	`.fit(cov_type="robust", debiased=True)`

For example, OLS with the same controls as the IV model is:

ak_ols = smf.ols(
    "lnwkwage ~ educ + " + controls, data=data
).fit()

Changing the covariance estimator changes the standard errors and associated tests. It leaves the coefficient estimates unchanged. State the choice in the table notes.

Test the excluded instruments jointly

first_stage = smf.ols(
    "educ ~ q1 + q2 + q3 + " + controls, data=data
).fit()

print(first_stage.summary())
first_stage_test = first_stage.f_test("q1 = q2 = q3 = 0")
print(first_stage_test)
first_stage_f = float(first_stage_test.fvalue)

This tests whether the three quarter-of-birth coefficients are jointly zero, conditional on the controls. The F-statistic at the top of the full regression summary tests all slopes together, including controls. Report the excluded-instrument test when assessing instrument relevance.

Calculate partial $R^2$

Estimate a restricted first stage that includes only the controls, using the same observations:

restricted = smf.ols("educ ~ " + controls, data=data).fit()
partial_r2 = (restricted.ssr - first_stage.ssr) / restricted.ssr
print(partial_r2)

.ssr is the sum of squared residuals. Partial $R^2$ measures how much of the residual variation in schooling is explained by adding the excluded instruments. It is reported here as a proportion. Multiply by 100 for the scale used in BJB’s tables.

Construct an interaction

Multiply indicators to identify observations meeting both conditions. For example:

data["q1_y1930"] = data["q1"] * (data["yearborn"] == 1930)

This variable equals 1 for first-quarter births in 1930 and 0 otherwise. Quarter-by-year interactions let the relationship between birth quarter and schooling differ across cohorts. Part 6 of PS2 asks you to construct the full instrument set with Claude and check which columns are independent once the controls are included.

Extract results for your own tables

For a coefficient named D, the main attributes are:

Result	`statsmodels` fitted result `ols_result`	`linearmodels` fitted result `iv_result`
Coefficient	`ols_result.params["D"]`	`iv_result.params["D"]`
Standard error	`ols_result.bse["D"]`	`iv_result.std_errors["D"]`
95% confidence interval	`ols_result.conf_int().loc["D"]`	`iv_result.conf_int().loc["D"]`
Sample size	`int(ols_result.nobs)`	`int(iv_result.nobs)`

Replace D with the coefficient name in your model, such as educ or medicaid. Report sample sizes as integers. Four decimal places are usually sufficient for coefficients and standard errors. Retain enough precision to display small partial $R^2$ values clearly.

Week 7: Regression discontinuity

Week 7 notes. In these examples, X is the running variable, Y is the outcome, and the cutoff is zero. Use a dataset with these variables observed on the same rows.

Local linear regression

Centre the running variable at the cutoff and choose a bandwidth in its units. Here the illustrative bandwidth is 0.10. Fit separate slopes on the two sides, with triangular weights:

import statsmodels.formula.api as smf

local = data.loc[data["X"].abs() < 0.10].copy()
local["eligible"] = (local["X"] >= 0).astype(int)
local["weight"] = 1 - local["X"].abs() / 0.10
sharp_rd = smf.wls(
    "Y ~ eligible * X", data=local, weights=local["weight"]
).fit(cov_type="HC1")
print(sharp_rd.summary())

The coefficient on eligible estimates the jump at zero. eligible * X includes both variables and their interaction. HC1 accounts for heteroskedasticity. For bias-corrected RD inference and bandwidth selection, use rdrobust below.

Bandwidth selection and robust RD inference

import rdrobust as rd

sharp = rd.rdrobust(
    y=data["Y"], x=data["X"], c=0,
    p=1, q=2, kernel="tri", bwselect="mserd"
)
print(sharp)

p=1 fits local straight lines, q=2 estimates their bias using a quadratic, and kernel="tri" gives greater weight to observations closer to the cutoff. Report the conventional point estimate, bias-corrected estimate, robust confidence interval, bandwidths and observations used on each side. For Lee’s data, add cluster=data["statedisdec"] to retain the clustering used in the exercise. Optional baseline controls enter through covs=data[["X1", "X2"]].

Fuzzy RD

Add actual treatment receipt, D, to estimate the ratio of the outcome jump to the treatment-probability jump:

fuzzy = rd.rdrobust(
    y=data["Y"], x=data["X"], fuzzy=data["D"], c=0,
    p=1, q=2, kernel="tri", bwselect="mserd"
)
print(fuzzy)

Density test at the cutoff

import rddensity as density

density_test = density.rddensity(X=data["X"], c=0)
print(density_test)

The null hypothesis is that the running-variable density is continuous at the cutoff. Interpret the test alongside a histogram and evidence about the assignment process. A large p-value leaves uncertainty about manipulation. Package documentation.

Week 8: Basic difference-in-differences

Week 8 notes. Here T identifies the treatment group, post identifies the post-policy period, and unit identifies the unit observed in both periods.

Four means

means = data.groupby(["T", "post"])["Y"].mean()
did_means = (means.loc[1, 1] - means.loc[1, 0]) - (
    means.loc[0, 1] - means.loc[0, 0]
)
print(did_means)

Regression in levels

did = smf.ols("Y ~ T * post", data=data).fit(
    cov_type="cluster", cov_kwds={"groups": data["unit"]}
)
print(did.summary())

T:post is the difference-in-differences estimate. Clustering by unit allows the two disturbances for that unit to be correlated. When policy is assigned at a broader level, consider dependence at that level too. Two states alone provide very limited information about state-level uncertainty.

Baseline covariates and first differences

Let X1 and X2 be measured before treatment and repeated on each unit’s two rows. Interact them with post to allow their associations with the outcome to change:

adjusted = smf.ols(
    "Y ~ T * post + X1 + X2 + post:X1 + post:X2", data=data
).fit(cov_type="cluster", cov_kwds={"groups": data["unit"]})

wide = data.pivot(index="unit", columns="post", values="Y")
wide["change"] = wide[1] - wide[0]
baseline = data.loc[data["post"] == 0, ["unit", "T", "X1", "X2"]]
changes = baseline.merge(wide[["change"]], on="unit").dropna()
first_difference = smf.ols("change ~ T + X1 + X2", data=changes).fit(
    cov_type="HC1"
)
print(first_difference.summary())

For a balanced two-period panel using the same complete observations, the coefficient on T in differences equals that on T:post in levels. The covariate coefficients in differences correspond to the covariate-by-post coefficients in levels.

Week 9: TWFE and event studies

Week 9 notes. D is treatment status in each unit-year. year always denotes calendar time.

twfe = smf.ols("Y ~ D + C(unit) + C(year)", data=data).fit(
    cov_type="cluster", cov_kwds={"groups": data["unit"]}
)

For an event study, first record each treated unit’s adoption year as adoption_year. Leave it missing for never-treated units. Create indicators for each event year in the observed window, omitting year −1. This example assumes all treated observations lie between −3 and 3:

data["event_time"] = data["year"] - data["adoption_year"]
data["lead3"] = (data["event_time"] == -3).astype(int)
data["lead2"] = (data["event_time"] == -2).astype(int)
data["event0"] = (data["event_time"] == 0).astype(int)
data["lag1"] = (data["event_time"] == 1).astype(int)
data["lag2"] = (data["event_time"] == 2).astype(int)
data["lag3"] = (data["event_time"] == 3).astype(int)
event_study = smf.ols(
    "Y ~ lead3 + lead2 + event0 + lag1 + lag2 + lag3 + C(unit) + C(year)",
    data=data
).fit(cov_type="cluster", cov_kwds={"groups": data["unit"]})

Extend the indicators to cover the observed event times or explicitly bin the endpoints. The Week 9 exercise uses the authors’ bins and reference periods. With staggered adoption, treatment-effect heterogeneity can complicate TWFE event-study coefficients. We return to this in Week 10.

Week 10: Modern difference-in-differences

Week 10 notes · Exercise and data.

County-clustered standard errors

Observations from the same county in different years may share common shocks. Clustered standard errors allow for this dependence when measuring uncertainty. This changes standard errors and confidence intervals. The estimated coefficients stay the same.

Using the Week 10 data, construct treatment status and estimate TWFE as follows:

import pandas as pd
import statsmodels.formula.api as smf

data = pd.read_csv("mpdta.csv")
data = data.rename(columns={"first.treat": "first_treat"})
data["D"] = ((data["first_treat"] > 0) &
             (data["year"] >= data["first_treat"])).astype(int)

twfe = smf.ols(
    "lemp ~ D + C(countyreal) + C(year)", data=data
).fit(
    cov_type="cluster",                       # Cluster the standard errors.
    cov_kwds={"groups": data["countyreal"]}    # County identifies each cluster.
)
print(twfe.summary())

A cohort-time comparison

To estimate the effect in 2007 for counties first treated in 2004, compare their change since 2003 with that of never-treated counties. Here, zero in first_treat denotes never treated.

sample = data.loc[data["first_treat"].isin([0, 2004])].copy()
wide = sample.pivot(index="countyreal", columns="year", values="lemp")
cohorts = sample.drop_duplicates("countyreal").set_index("countyreal")

changes = pd.DataFrame(index=wide.index)
changes["change"] = wide[2007] - wide[2003]
changes["cohort2004"] = (cohorts["first_treat"] == 2004).astype(int)
changes = changes.dropna(subset=["change"])

means = changes.groupby("cohort2004")["change"].mean()
att_2004_2007 = means.loc[1] - means.loc[0]
print(att_2004_2007)

did_2004_2007 = smf.ols(
    "change ~ cohort2004", data=changes
).fit(cov_type="HC1")
print(did_2004_2007.summary())

There is one observation per county in this change regression, so we use HC1 standard errors. The coefficient on cohort2004 equals the difference in mean changes. Under parallel trends and no anticipation, it estimates $ATT(2004,2007)$.

Week 11: Synthetic controls

Week 11 notes · Worked exercise and setup notebook.

Use pysyncon for synthetic control estimation. Install it once in the Python environment used by Quarto:

python -m pip install pysyncon==1.7.0

Specify and fit the model

This template uses a panel data frame df, with outcome Y, predictors X1 and X2, a unit identifier unit, and calendar year year. For illustration, unit 1 receives treatment in 2010 and units 2–5 form the donor pool. Replace these names, identifiers and years with those appropriate for your application.

import pysyncon as sc

dataprep = sc.Dataprep(
    foo=df,
    predictors=["X1", "X2"],
    predictors_op="mean",
    time_predictors_prior=range(2000, 2010),
    dependent="Y",
    unit_variable="unit",
    time_variable="year",
    treatment_identifier=1,
    controls_identifier=[2, 3, 4, 5],
    time_optimize_ssr=range(2000, 2010),
)

synth = sc.Synth()
synth.fit(dataprep=dataprep)

time_predictors_prior specifies the years over which predictors are averaged. time_optimize_ssr specifies the outcome years used to choose predictor weights. Here both use 2000–2009: Python’s range excludes its upper endpoint.

Inspect the fit and plot the results

Task	Python
Donor weights	`synth.weights()`
Predictor weights	`synth.V`
Predictor balance	`synth.summary()`

The donor weights should be nonnegative and sum to one. Compare the treated and synthetic predictor values and inspect the pre-treatment outcome fit.

synth.path_plot(
    time_period=range(2000, 2021), treatment_time=2010
)

Plot the gap, defined as observed minus synthetic outcome:

synth.gaps_plot(
    time_period=range(2000, 2021), treatment_time=2010
)

The Week 11 exercise provides the German application, including separate training and fitting periods, donor sensitivity checks, and placebo comparisons.

--- title: "Python: Estimating Effects" --- A quick reference for the methods used in class, building on PP5000. We will add methods as we introduce them. See the [weekly notes](index.qmd) for the reasoning behind the estimates. ## Week 1: Differences in means and regression These syntax templates use a data frame called `df`, an outcome `Y`, and a treatment indicator `D` coded 1 for treatment and 0 for control. Replace these names with the variables in your dataset. ```python import pandas as pd import statsmodels.formula.api as smf ``` | Task | Python | What to read | |:---|:---|:---| | Group means | `sample.groupby("D")["Y"].mean()` | Mean outcome in each group | | Difference in means | `means.loc[1] - means.loc[0]` | Treatment mean minus control mean | | Regression | `smf.ols("Y ~ D", data=sample)` | Fit the model below; read the coefficient on `D` | | Region indicators | `C(region)` in the formula | Adjust for categorical region controls | ### Difference in means Use observations with both treatment and outcome observed, so the means and regression refer to the same sample. ```python sample = df.dropna(subset=["Y", "D"]) means = sample.groupby("D")["Y"].mean() difference = means.loc[1] - means.loc[0] print(means) print(difference) ``` The difference is in the units of `Y`. With random assignment and complete outcome data, it estimates the effect of assignment to treatment. Without a credible research design, a difference in means need not be a causal effect. ### Regression and standard errors ::: {.callout-important} ## Use heteroskedasticity-consistent standard errors For the regressions in Week 1, **use heteroskedasticity-consistent standard errors**, also called robust standard errors. Specify `.fit(cov_type="HC1", use_t=True)` when fitting the model. Calling `.fit()` alone uses the conventional standard errors that assume homoskedastic disturbances. Later exercises may require clustered standard errors instead; follow the instructions for the research design. ::: The complete specification for a regression of the outcome on an intercept and the treatment indicator is: ```python import statsmodels.formula.api as smf sample = df.dropna(subset=["Y", "D"]) model = smf.ols("Y ~ D", data=sample).fit( cov_type="HC1", use_t=True ) print(model.summary()) ``` The formula includes an intercept by default. Its estimate is the control-group mean; the coefficient on `D` equals the treatment-group mean minus the control-group mean. `cov_type="HC1"` requests **heteroskedasticity-consistent standard errors**. It changes the estimated uncertainty, not the estimated coefficients. `use_t=True` requests t-based inference, as in the Week 1 exercise. ```python estimate = model.params["D"] standard_error = model.bse["D"] confidence_interval = model.conf_int().loc["D"] n_observations = model.nobs ``` Report the estimate, its units, its standard error or confidence interval, and the number of observations. HC1 allows unequal disturbance variances; it does not account for dependence within groups. We will revisit standard errors for other research designs. ### Region indicators Here is the complete specification with region fixed effects and heteroskedasticity-consistent standard errors: ```python import statsmodels.formula.api as smf sample_region = df.dropna(subset=["Y", "D", "region"]) model_region = smf.ols( formula="Y ~ D + C(region)", data=sample_region ).fit(cov_type="HC1", use_t=True) print(model_region.summary()) ``` `C(region)` treats region codes as categories and includes indicators, omitting a reference category. Writing `region` alone treats the codes as a quantitative variable. The coefficient on `D` is now the regression-adjusted difference; it need not equal the raw difference in means. Missing values in added controls can also change the estimation sample. ### Binary outcomes: the linear probability model When `Y` is coded 0 or 1, this OLS regression is a **linear probability model**. A coefficient of 0.04 means a difference of **4 percentage points**, not a 4 per cent increase. If the conditional probability is $p$, the conditional variance of the binary outcome—and hence of its disturbance around the conditional mean—is $p\left(1-p\right)$. This generally varies across observations, motivating heteroskedasticity-consistent standard errors. To express an estimate in percentage points, multiply the estimate in probability units by 100. Scale its standard error and confidence-interval endpoints by the same factor. ## Week 2: Propensity scores and weighting Use an outcome `Y`, a treatment indicator `D` coded 1 for treatment and 0 for control, and pre-treatment characteristics `X1` and `X2`. Replace these names with your variables. The weights below estimate the **average treatment effect on the treated (ATT)** under the assumptions discussed in the [Week 2 notes](sessions/02-matching.qmd). | Task | Python | What to read | |:---|:---|:---| | Estimate a logit | `smf.logit("D ~ X1 + X2", data=sample_ipw).fit()` | Model the probability of treatment given the characteristics | | Predict propensity scores | `score_model.predict(sample_ipw)` | Estimated treatment probability for each observation | | ATT weights | `1` for treated; `pscore / (1 - pscore)` for controls | Give more weight to controls whose characteristics are common among treated observations | | Weighted regression | `smf.wls("Y ~ D", data=sample_ipw, weights=sample_ipw["att_weight"])` | Fit the model below; read the coefficient on `D` | ### Estimate the propensity score Keep a common sample for the propensity-score model and outcome regression. ```python import statsmodels.formula.api as smf sample_ipw = df.dropna(subset=["Y", "D", "X1", "X2"]).copy() score_model = smf.logit( "D ~ X1 + X2", data=sample_ipw ).fit() sample_ipw["pscore"] = score_model.predict(sample_ipw) ``` `smf.logit` estimates the logit model. `.predict()` returns the estimated probability of treatment for each observation: its propensity score. For this step, we use those probabilities to construct the weights. ### Construct ATT weights Treated observations receive weight 1. Controls receive $\hat p_i / \left(1-\hat p_i\right)$. ```python # Start with a weight of 1 for every observation. sample_ipw["att_weight"] = 1.0 # Replace the weights for the control observations. controls = sample_ipw["D"] == 0 sample_ipw.loc[controls, "att_weight"] = ( sample_ipw.loc[controls, "pscore"] / (1 - sample_ipw.loc[controls, "pscore"]) ) ``` `controls` identifies the rows where `D` equals zero. `.loc[controls, ...]` selects those rows. Check overlap in propensity scores and covariate balance, and inspect large weights, before interpreting the estimate. ### Estimate the weighted regression The complete outcome regression includes an intercept and the treatment indicator: ```python weighted = smf.wls( "Y ~ D", data=sample_ipw, weights=sample_ipw["att_weight"], ).fit(cov_type="HC1", use_t=True) estimate = weighted.params["D"] standard_error = weighted.bse["D"] confidence_interval = weighted.conf_int().loc["D"] n_observations = weighted.nobs ``` The coefficient on `D` is the treated-group mean minus the weighted control-group mean, in the units of `Y`. `cov_type="HC1"` requests **heteroskedasticity-consistent standard errors**. These standard errors treat the estimated weights as fixed; a full calculation of uncertainty also accounts for estimating the propensity scores. For [Problem Set 1](exercises/02-matching.qmd), use `re78` as the outcome and `treat` as the treatment indicator. Estimate a separate propensity-score model for each observational comparison sample, using the characteristics specified in the assignment. ## Week 3: Instrumental variables Use `linearmodels` for two-stage least squares. With outcome `Y`, endogenous treatment `D`, instrument `Z`, and a categorical control `group`: ```python import linearmodels.iv as iv sample_iv = df.dropna(subset=["Y", "D", "Z", "group"]) model_iv = iv.IV2SLS.from_formula( "Y ~ 1 + C(group) + [D ~ Z]", data=sample_iv ) result_iv = model_iv.fit(cov_type="robust", debiased=True) estimate = result_iv.params["D"] standard_error = result_iv.std_errors["D"] confidence_interval = result_iv.conf_int().loc["D"] ``` `[D ~ Z]` tells the estimator to instrument `D` with `Z`. `C(group)` includes category indicators in both stages. `cov_type="robust"` requests heteroskedasticity-consistent IV standard errors. In this package, the standard-error attribute is `.std_errors`. The [Week 3 exercise](exercises/03-iv-late.qmd) uses the Oregon survey weights and clusters by household: ```python model_oregon = iv.IV2SLS.from_formula( "outpatient ~ 1 + C(stratum) + [medicaid ~ selected]", data=analysis, weights=analysis["weight_12m"], ) result_oregon = model_oregon.fit( cov_type="clustered", clusters=analysis["household_id"], debiased=True, ) ``` Here `stratum` identifies the observed combinations of household-list size and survey wave. The setup file constructs it and selects the common analysis sample. Household clustering permits heteroskedasticity and correlation within households. ### First stage, reduced form and the IV ratio With one excluded instrument, estimate the first stage and reduced form using the same observations, controls and weights as IV. For the Oregon example: ```python first_oregon = smf.wls( "medicaid ~ selected + C(stratum)", data=analysis, weights=analysis["weight_12m"], ).fit( cov_type="cluster", cov_kwds={"groups": analysis["household_id"]}, use_t=True, ) reduced_oregon = smf.wls( "outpatient ~ selected + C(stratum)", data=analysis, weights=analysis["weight_12m"], ).fit( cov_type="cluster", cov_kwds={"groups": analysis["household_id"]}, use_t=True, ) first_stage_effect = first_oregon.params["selected"] itt_effect = reduced_oregon.params["selected"] iv_ratio = itt_effect / first_stage_effect print(iv_ratio) print(result_oregon.params["medicaid"]) ``` The first stage measures the effect of lottery selection on Medicaid coverage. The reduced form measures its effect on the outcome, also called the intention-to-treat effect here. Their coefficient ratio equals the IV coefficient in this exactly identified model. Obtain its standard error from the fitted IV model. ### Display the first stage from an IV model ```python print(result_iv.first_stage) print(result_iv.first_stage.diagnostics) ``` The first command displays the first-stage coefficients and diagnostics. The second gives a data frame of diagnostics, including partial $R^2$ and a test of the excluded instruments. Read `f.dist` alongside `f.stat`: with robust or clustered covariance, `linearmodels` reports a chi-squared Wald statistic for this test. The examples below show how to obtain the conventional first-stage F used in PS2. ## Week 4: Leniency instruments and simulations ### Construct leave-one-out leniency In the [Week 4 simulation](exercises/04-judges-examiners.qmd), `caseworker` identifies the decision-maker and `training` is 1 for admission and 0 otherwise. The data contain one observation per applicant and no missing admission decisions. | Task | Python | Result | |:---|:---|:---| | Count admitted applicants by caseworker | `data.groupby("caseworker")["training"].transform("sum")` | Each applicant receives their caseworker's total admissions | | Count applicants by caseworker | `data.groupby("caseworker")["training"].transform("count")` | Each applicant receives their caseworker's caseload | | Leave out the applicant | `(admissions - data["training"]) / (caseload - 1)` | Admission rate among the caseworker's other applicants | ```python admissions = data.groupby("caseworker")["training"].transform("sum") caseload = data.groupby("caseworker")["training"].transform("count") data["leniency"] = ( (admissions - data["training"]) / (caseload - 1) ) ``` `transform()` returns a value for each original row, making it convenient for constructing a new variable. Subtracting the applicant's own admission decision removes its direct contribution to their instrument. This calculation requires at least two applicants per caseworker. ### Estimate IV with leniency The simulation assigns applicants across caseworkers without additional assignment groups: ```python import linearmodels.iv as iv leniency_iv = iv.IV2SLS.from_formula( "earnings_after ~ 1 + [training ~ leniency]", data=data, ).fit(cov_type="robust", debiased=True) print(leniency_iv.summary) ``` In an empirical application, assignment may be comparable only within an office, court or time period. If `assignment_group` identifies those groups, include their indicators in both stages: ```python leniency_group_iv = iv.IV2SLS.from_formula( "Y ~ 1 + C(assignment_group) + [D ~ leniency]", data=assignment_sample, ).fit(cov_type="robust", debiased=True) ``` Here `assignment_sample` is the application's common analysis sample, with outcome `Y`, treatment `D` and its constructed leniency instrument. The definition of leniency and the choice of clustered standard errors should follow the assignment process in that application. ### Keep a simulated sample reproducible ```python import numpy as np rng = np.random.default_rng(12345) ``` Use this generator for the random draws in the simulation. Re-running the same code with the same seed reproduces the sample. Choose your own integer. `np.random.default_rng()` with empty parentheses produces a new sample each time. ## Week 5: Multiple instruments and first-stage evidence The examples below use the AK variable names from [Problem Set 2](exercises/05-weak-instruments.qmd). The setup file constructs `q1`, `q2`, `q3`, `age_years` and `age_squared`. Quarter 4 is the omitted category. ### Estimate IV with several instruments ```python import statsmodels.formula.api as smf import linearmodels.iv as iv controls = "black + married + smsa + C(division) + age_years + age_squared" ak_iv = iv.IV2SLS.from_formula( "lnwkwage ~ 1 + " + controls + " + [educ ~ q1 + q2 + q3]", data=data, ).fit(cov_type="unadjusted", debiased=True) print(ak_iv.summary) ``` The expression `[educ ~ q1 + q2 + q3]` specifies schooling as endogenous and the three quarter indicators as excluded instruments. Included controls enter both stages. ### Conventional standard errors for replication PS2 uses conventional standard errors to compare with AK and BJB. These assume homoskedastic disturbances. | Covariance choice | `statsmodels` OLS | `linearmodels` IV | |:---|:---|:---| | Conventional | `.fit()` | `.fit(cov_type="unadjusted", debiased=True)` | | Heteroskedasticity-consistent | `.fit(cov_type="HC1", use_t=True)` | `.fit(cov_type="robust", debiased=True)` | For example, OLS with the same controls as the IV model is: ```python ak_ols = smf.ols( "lnwkwage ~ educ + " + controls, data=data ).fit() ``` Changing the covariance estimator changes the standard errors and associated tests. It leaves the coefficient estimates unchanged. State the choice in the table notes. ### Test the excluded instruments jointly ```python first_stage = smf.ols( "educ ~ q1 + q2 + q3 + " + controls, data=data ).fit() print(first_stage.summary()) first_stage_test = first_stage.f_test("q1 = q2 = q3 = 0") print(first_stage_test) first_stage_f = float(first_stage_test.fvalue) ``` This tests whether the three quarter-of-birth coefficients are jointly zero, conditional on the controls. The F-statistic at the top of the full regression summary tests all slopes together, including controls. Report the excluded-instrument test when assessing instrument relevance. ### Calculate partial $R^2$ Estimate a restricted first stage that includes only the controls, using the same observations: ```python restricted = smf.ols("educ ~ " + controls, data=data).fit() partial_r2 = (restricted.ssr - first_stage.ssr) / restricted.ssr print(partial_r2) ``` `.ssr` is the sum of squared residuals. Partial $R^2$ measures how much of the residual variation in schooling is explained by adding the excluded instruments. It is reported here as a proportion. Multiply by 100 for the scale used in BJB's tables. ### Construct an interaction Multiply indicators to identify observations meeting both conditions. For example: ```python data["q1_y1930"] = data["q1"] * (data["yearborn"] == 1930) ``` This variable equals 1 for first-quarter births in 1930 and 0 otherwise. Quarter-by-year interactions let the relationship between birth quarter and schooling differ across cohorts. Part 6 of PS2 asks you to construct the full instrument set with Claude and check which columns are independent once the controls are included. ## Extract results for your own tables For a coefficient named `D`, the main attributes are: | Result | `statsmodels` fitted result `ols_result` | `linearmodels` fitted result `iv_result` | |:---|:---|:---| | Coefficient | `ols_result.params["D"]` | `iv_result.params["D"]` | | Standard error | `ols_result.bse["D"]` | `iv_result.std_errors["D"]` | | 95% confidence interval | `ols_result.conf_int().loc["D"]` | `iv_result.conf_int().loc["D"]` | | Sample size | `int(ols_result.nobs)` | `int(iv_result.nobs)` | Replace `D` with the coefficient name in your model, such as `educ` or `medicaid`. Report sample sizes as integers. Four decimal places are usually sufficient for coefficients and standard errors. Retain enough precision to display small partial $R^2$ values clearly. ## Week 7: Regression discontinuity [Week 7 notes](sessions/07-rdd.qmd). In these examples, `X` is the running variable, `Y` is the outcome, and the cutoff is zero. Use a dataset with these variables observed on the same rows. ### Local linear regression Centre the running variable at the cutoff and choose a bandwidth in its units. Here the illustrative bandwidth is 0.10. Fit separate slopes on the two sides, with triangular weights: ```python import statsmodels.formula.api as smf local = data.loc[data["X"].abs() < 0.10].copy() local["eligible"] = (local["X"] >= 0).astype(int) local["weight"] = 1 - local["X"].abs() / 0.10 sharp_rd = smf.wls( "Y ~ eligible * X", data=local, weights=local["weight"] ).fit(cov_type="HC1") print(sharp_rd.summary()) ``` The coefficient on `eligible` estimates the jump at zero. `eligible * X` includes both variables and their interaction. HC1 accounts for heteroskedasticity. For bias-corrected RD inference and bandwidth selection, use `rdrobust` below. ### Bandwidth selection and robust RD inference ```python import rdrobust as rd sharp = rd.rdrobust( y=data["Y"], x=data["X"], c=0, p=1, q=2, kernel="tri", bwselect="mserd" ) print(sharp) ``` `p=1` fits local straight lines, `q=2` estimates their bias using a quadratic, and `kernel="tri"` gives greater weight to observations closer to the cutoff. Report the conventional point estimate, bias-corrected estimate, robust confidence interval, bandwidths and observations used on each side. For Lee's data, add `cluster=data["statedisdec"]` to retain the clustering used in the exercise. Optional baseline controls enter through `covs=data[["X1", "X2"]]`. ### Fuzzy RD Add actual treatment receipt, `D`, to estimate the ratio of the outcome jump to the treatment-probability jump: ```python fuzzy = rd.rdrobust( y=data["Y"], x=data["X"], fuzzy=data["D"], c=0, p=1, q=2, kernel="tri", bwselect="mserd" ) print(fuzzy) ``` ### Density test at the cutoff ```python import rddensity as density density_test = density.rddensity(X=data["X"], c=0) print(density_test) ``` The null hypothesis is that the running-variable density is continuous at the cutoff. Interpret the test alongside a histogram and evidence about the assignment process. A large p-value leaves uncertainty about manipulation. [Package documentation](https://rdpackages.github.io/). ## Week 8: Basic difference-in-differences [Week 8 notes](sessions/08-diff-in-diff.qmd). Here `T` identifies the treatment group, `post` identifies the post-policy period, and `unit` identifies the unit observed in both periods. ### Four means ```python means = data.groupby(["T", "post"])["Y"].mean() did_means = (means.loc[1, 1] - means.loc[1, 0]) - ( means.loc[0, 1] - means.loc[0, 0] ) print(did_means) ``` ### Regression in levels ```python did = smf.ols("Y ~ T * post", data=data).fit( cov_type="cluster", cov_kwds={"groups": data["unit"]} ) print(did.summary()) ``` `T:post` is the difference-in-differences estimate. Clustering by unit allows the two disturbances for that unit to be correlated. When policy is assigned at a broader level, consider dependence at that level too. Two states alone provide very limited information about state-level uncertainty. ### Baseline covariates and first differences Let `X1` and `X2` be measured before treatment and repeated on each unit's two rows. Interact them with `post` to allow their associations with the outcome to change: ```python adjusted = smf.ols( "Y ~ T * post + X1 + X2 + post:X1 + post:X2", data=data ).fit(cov_type="cluster", cov_kwds={"groups": data["unit"]}) wide = data.pivot(index="unit", columns="post", values="Y") wide["change"] = wide[1] - wide[0] baseline = data.loc[data["post"] == 0, ["unit", "T", "X1", "X2"]] changes = baseline.merge(wide[["change"]], on="unit").dropna() first_difference = smf.ols("change ~ T + X1 + X2", data=changes).fit( cov_type="HC1" ) print(first_difference.summary()) ``` For a balanced two-period panel using the same complete observations, the coefficient on `T` in differences equals that on `T:post` in levels. The covariate coefficients in differences correspond to the covariate-by-post coefficients in levels. ## Week 9: TWFE and event studies [Week 9 notes](sessions/09-twfe-event-studies.qmd). `D` is treatment status in each unit-year. `year` always denotes calendar time. ```python twfe = smf.ols("Y ~ D + C(unit) + C(year)", data=data).fit( cov_type="cluster", cov_kwds={"groups": data["unit"]} ) ``` For an event study, first record each treated unit's adoption year as `adoption_year`. Leave it missing for never-treated units. Create indicators for each event year in the observed window, omitting year −1. This example assumes all treated observations lie between −3 and 3: ```python data["event_time"] = data["year"] - data["adoption_year"] data["lead3"] = (data["event_time"] == -3).astype(int) data["lead2"] = (data["event_time"] == -2).astype(int) data["event0"] = (data["event_time"] == 0).astype(int) data["lag1"] = (data["event_time"] == 1).astype(int) data["lag2"] = (data["event_time"] == 2).astype(int) data["lag3"] = (data["event_time"] == 3).astype(int) event_study = smf.ols( "Y ~ lead3 + lead2 + event0 + lag1 + lag2 + lag3 + C(unit) + C(year)", data=data ).fit(cov_type="cluster", cov_kwds={"groups": data["unit"]}) ``` Extend the indicators to cover the observed event times or explicitly bin the endpoints. The [Week 9 exercise](exercises/09-twfe-event-studies.qmd) uses the authors' bins and reference periods. With staggered adoption, treatment-effect heterogeneity can complicate TWFE event-study coefficients. We return to this in Week 10. ## Week 10: Modern difference-in-differences [Week 10 notes](sessions/10-modern-did.qmd) · [Exercise and data](exercises/10-modern-did.qmd). ### County-clustered standard errors Observations from the same county in different years may share common shocks. Clustered standard errors allow for this dependence when measuring uncertainty. This changes standard errors and confidence intervals. The estimated coefficients stay the same. Using the Week 10 data, construct treatment status and estimate TWFE as follows: ```python import pandas as pd import statsmodels.formula.api as smf data = pd.read_csv("mpdta.csv") data = data.rename(columns={"first.treat": "first_treat"}) data["D"] = ((data["first_treat"] > 0) & (data["year"] >= data["first_treat"])).astype(int) twfe = smf.ols( "lemp ~ D + C(countyreal) + C(year)", data=data ).fit( cov_type="cluster", # Cluster the standard errors. cov_kwds={"groups": data["countyreal"]} # County identifies each cluster. ) print(twfe.summary()) ``` ### A cohort-time comparison To estimate the effect in 2007 for counties first treated in 2004, compare their change since 2003 with that of never-treated counties. Here, zero in `first_treat` denotes never treated. ```python sample = data.loc[data["first_treat"].isin([0, 2004])].copy() wide = sample.pivot(index="countyreal", columns="year", values="lemp") cohorts = sample.drop_duplicates("countyreal").set_index("countyreal") changes = pd.DataFrame(index=wide.index) changes["change"] = wide[2007] - wide[2003] changes["cohort2004"] = (cohorts["first_treat"] == 2004).astype(int) changes = changes.dropna(subset=["change"]) means = changes.groupby("cohort2004")["change"].mean() att_2004_2007 = means.loc[1] - means.loc[0] print(att_2004_2007) did_2004_2007 = smf.ols( "change ~ cohort2004", data=changes ).fit(cov_type="HC1") print(did_2004_2007.summary()) ``` There is one observation per county in this change regression, so we use HC1 standard errors. The coefficient on `cohort2004` equals the difference in mean changes. Under parallel trends and no anticipation, it estimates $ATT(2004,2007)$. ## Week 11: Synthetic controls [Week 11 notes](sessions/11-synthetic-controls.qmd) · [Worked exercise and setup notebook](exercises/11-synthetic-controls.qmd). Use `pysyncon` for synthetic control estimation. Install it once in the Python environment used by Quarto: ```bash python -m pip install pysyncon==1.7.0 ``` ### Specify and fit the model This template uses a panel data frame `df`, with outcome `Y`, predictors `X1` and `X2`, a unit identifier `unit`, and calendar year `year`. For illustration, unit 1 receives treatment in 2010 and units 2–5 form the donor pool. Replace these names, identifiers and years with those appropriate for your application. ```python import pysyncon as sc dataprep = sc.Dataprep( foo=df, predictors=["X1", "X2"], predictors_op="mean", time_predictors_prior=range(2000, 2010), dependent="Y", unit_variable="unit", time_variable="year", treatment_identifier=1, controls_identifier=[2, 3, 4, 5], time_optimize_ssr=range(2000, 2010), ) synth = sc.Synth() synth.fit(dataprep=dataprep) ``` `time_predictors_prior` specifies the years over which predictors are averaged. `time_optimize_ssr` specifies the outcome years used to choose predictor weights. Here both use 2000–2009: Python's `range` excludes its upper endpoint. ### Inspect the fit and plot the results | Task | Python | |:---|:---| | Donor weights | `synth.weights()` | | Predictor weights | `synth.V` | | Predictor balance | `synth.summary()` | The donor weights should be nonnegative and sum to one. Compare the treated and synthetic predictor values and inspect the pre-treatment outcome fit. ```python synth.path_plot( time_period=range(2000, 2021), treatment_time=2010 ) ``` Plot the gap, defined as observed minus synthetic outcome: ```python synth.gaps_plot( time_period=range(2000, 2021), treatment_time=2010 ) ``` The [Week 11 exercise](exercises/11-synthetic-controls.qmd) provides the German application, including separate training and fitting periods, donor sensitivity checks, and placebo comparisons.