ML for SS: ML for Causality

Dr. Richard M. Crowley

rcrowley@smu.edu.sg

https://rmc.link/

Overview

Papers

Chernozhukov et al. 2017 AER (including web appendix)

• Introduces a ML-based method for causal identification useful in standard DID and IV approaches
  • Focused on calculating ATE and ATTE

Gentzkow, Shapiro and Taddy (2019) Econometrica

• Shows the methodological benefits that can come from careful merging of econometrics and machine learning

Davis and Heller (2020) RES

• Applies causal forests to randomized field experiments
• Focused on heterogeneity of experimental effects

Technical Discussion: DoubleML

Python

• Use the DoubleML library
  • For a basic DoubleML model
  • Also handles clustering, DiD, and IV methods

R

• The doubleML library is available in R as well
  • We’ll try this for a model with clustered standard errors
• The AER paper’s source code is all in R!

Polars incompatability

Due to a lack of Stata .dta support in Polars, we’ll use Pandas today. There is experimental support for it in polars-readstat

Technical Discussion: Causal Random Forest

Python

• Use the EconML library
  • It can fit a generalized random forest (GRF)
  • It can also do a hybrid model with DML

R

• The grf package handles generalized random forest models

We will apply these models to the same data as DML

Double ML: Theory

Background

• There are a number of relevant papers published in economics in recent years developing and using Double ML
• The method is developed largely from:
  • Chernozhukov et al. (2017 AER), “Double/debiased/Neyman machine learning of treatment effects”
  • Chernozhukov et al. (2018 Econometrics J), “Double/debiased machine learning for treatment and structural parameters.”

Impact or overlap with methodological work by Susan Athey, Matthew Gentzkow, Trevor Hastie, Guido Imbens, Matt Taddy, and Stefan Wager

What is Double ML?

1. Split your sample as you would for \(K\)-fold cross validation, into sets \(\{I_k\}_{k\in \{1,...,K\}}\)
   • \(K\) samples of \(N/K\) observations each
   • Let \(I_k^c = \cup \{I_j\}_{j\ne k}\)
2. Construct \(K\) estimators using a machine learning estimator over nuisance parameters (e.g., controls) applied to the data \(I_K^c\)
3. Average the \(K\) estimators to obtain a final estimator
   • This average estimator is approximately unbiased and normally distributed
   • The estimator is also asymptotically efficient

And repeat. Bootstrap this out and take the mean or median of the estimators

Where Double ML excels: Endogenous treatment

• Suppose a policy affects a subset of individuals (people, corporations, etc.)
• Suppose individuals have the ability to alter their treatment status
  • E.g., state laws (move), labor laws, etc.
• Linear controls may be insufficient to claim causality of the treatment on anything

There are a lot of older methods that try to address this, though incompletely

1. Linear controls
2. Propensity score adjustments (e.g., weighting)
3. Matching methods
4. “doubly-robust” estimators

Why is machine learning needed?

• Suppose a true form of a specification is as follows
  • \(T\) is a treatment indicator, \(C\) is a vector of controls

\[ \begin{aligned} Y &= g_0(T, C) + \varepsilon_1\\ T &= m_0(C) + \varepsilon_2 \end{aligned} \]

• We often assume \(g_0\) to be something like \(\alpha + \theta_0\ T + \gamma \cdot C\)
• We often assume \(m_0\) to be a constant (i.e., assume that \(T\) is exogenous)

We know these assumptions aren’t true! (in many cases)

Why is machine learning needed?

How can we estimate a more general form for \(g_0\) and \(m_0\)?

• We could use a more flexible econometric approach, such as including interactions between \(T\) and \(C\)
  • This is still very restrictive: purely linear
• We could include transformations of \(C\) and its interactions
  • This is still restrictive: \(T\) is additively separable
• We could use a nonparametric estimator!
  • This is where machine learning is very useful: efficient and reasonably accurate nonparametric estimation
    • LASSO, random forest, XGBoost, etc.

Model variants

• Interactive regression model (IRM)
  • The model described in the previous slides
• Partially linear regression model (PLR)
  • Use if you can separate your treatment effect from the controls but suspect nonlinear effects of controls
    • Solves \(Y = \theta_0 T + g_0(C) + \varepsilon_0\) and \(T = m_0(C) + \varepsilon_2\)
• There are also instrumental variable variants of both IRM and PLR

What does this give us?

• Average treatment effect (ATE)
  • How does the treatment effect the outcome across groups?
  • \(\mathbb{E}\left[g_0(1,C) - g_0(0,C)\right]\)
• Average treatment effect of the treated (ATTE)
  • How does the treatment effect only those under the treatment?
  • \(\mathbb{E}\left[g_0(1,C) - g_0(0,C)|T=1\right]\)

Reconciling these slides notation with the paper

• These slides use a somewhat simpler oriented notation.
• Reconciliation from slides to papers:
  • \(T\) is \(D\)
  • \(C\) is \(X\)
  • \(\varepsilon_0\) is \(U\) or \(\zeta\) depending on the paper
  • \(\varepsilon_1\) is \(V\)

Implementing DoubleML

Walking through an implementation of DoubleML

Problem: How does 401k participation impact wealth?

• This problem is walked through in Chernozhukov et al. (2017 AER, Web Appendix)
  • The R code for the AER paper is available from AER as well
    • Quite clean code at that!
• We will implement this in python using the DoubleML library
  • Which Chernozhukov was involved in the development of

Importing the data

• Conveniently, the data is available from the DoubleML package

# Grab the dataset
import doubleml.datasets
df = dml.datasets.fetch_401K('DataFrame')
df

         nifa  net_tfa        tw  age      inc  ...  twoearn  e401  p401  pira  hown
0         0.0      0.0    4500.0   47   6765.0  ...        0     0     0     0     1
1      6215.0   1015.0   22390.0   36  28452.0  ...        0     0     0     0     1
2         0.0  -2000.0   -2000.0   37   3300.0  ...        0     0     0     0     0
3     15000.0  15000.0  155000.0   58  52590.0  ...        1     0     0     0     1
4         0.0      0.0   58000.0   32  21804.0  ...        0     0     0     0     1
...       ...      ...       ...  ...      ...  ...      ...   ...   ...   ...   ...
9910  98498.0  98858.0  157858.0   52  73920.0  ...        0     1     1     0     1
9911    287.0   6230.0   15730.0   41  42927.0  ...        1     1     1     1     1
9912     99.0   6099.0    7406.0   40  23619.0  ...        0     1     0     1     0
9913      0.0    -32.0    2468.0   47  14280.0  ...        0     1     1     0     0
9914   4000.0   5000.0    8857.0   33  11112.0  ...        0     1     1     0     0

[9915 rows x 14 columns]

Using your own data

• We can also do this manually, by importing the Stata file from AER
• We then need to prep the data into the format DoubleML expects
  • This is fairly straightforward, just defining our Y, treatment, and control variables

df = pd.read_stata('../../Data/S3_sipp1991.dta')

y = 'net_tfa'
treat = 'e401'
controls = [x for x in df.columns.tolist() if x not in [y, treat]]

df_dml = dml.DoubleMLData(df, y_col=y, d_cols=treat, x_cols=controls)

What is the data format used by DoubleML?

print(df_dml)

================== DoubleMLData Object ==================

------------------ Data summary      ------------------
Outcome variable: net_tfa
Treatment variable(s): ['e401']
Covariates: ['nifa', 'tw', 'age', 'inc', 'fsize', 'educ', 'db', 'marr', 'twoearn', 'p401', 'pira', 'hown']
Instrument variable(s): None
No. Observations: 9915
------------------ DataFrame info    ------------------
<class 'pandas.DataFrame'>
RangeIndex: 9915 entries, 0 to 9914
Columns: 14 entries, nifa to hown
dtypes: float32(4), int8(10)
memory usage: 251.9 KB

• Pandas dataframe
• A pre-specified outcome variable
• One or more treatment indicators
• One or more controls
• Optional instruments

Set up the Nuisance functions

• Recall that there are two functions, \(m_0\) and \(g_0\) that need to be solved for this method
  • \(m_0\) must match the treatment’s form, \(g_0\) the DV’s form
• We can specify any form for these that we want, if they follow Scikit-learn’s API

\(g_0\): Continuous GBM

g_0 = GradientBoostingRegressor(
  loss='squared_error',
  learning_rate=0.01,
  n_estimators=1000,
  subsample=0.5,
  max_depth=2
  )

\(m_0\): Binary GBM

m_0 = GradientBoostingClassifier(
  loss='exponential',
  learning_rate=0.01,
  n_estimators=1000,
  subsample=0.5,
  max_depth=2
  )

Run the DML model: Average Treatment Effects

# Fix the random number generator for replicability
np.random.seed(1234)
# Run the model
dml_model_irm = dml.DoubleMLIRM(df_dml, g_0, m_0)
# Output the model's findings
print(dml_model_irm.fit())

================== DoubleMLIRM Object ==================

------------------ Data summary      ------------------
Outcome variable: net_tfa
Treatment variable(s): ['e401']
Covariates: ['nifa', 'tw', 'age', 'inc', 'fsize', 'educ', 'db', 'marr', 'twoearn', 'p401', 'pira', 'hown']
Instrument variable(s): None
No. Observations: 9915

------------------ Score & algorithm ------------------
Score function: ATE

------------------ Machine learner   ------------------
Learner ml_g: GradientBoostingRegressor(learning_rate=0.01, max_depth=2, n_estimators=1000,
                          subsample=0.5)
Learner ml_m: GradientBoostingClassifier(learning_rate=0.01, loss='exponential', max_depth=2,
                           n_estimators=1000, subsample=0.5)
Out-of-sample Performance:
Regression:
Learner ml_g0 RMSE: [[11830.57121056]]
Learner ml_g1 RMSE: [[19910.06662348]]
Classification:
Learner ml_m Log Loss: [[0.28146494]]

------------------ Resampling        ------------------
No. folds: 5
No. repeated sample splits: 1

------------------ Fit summary       ------------------
             coef     std err         t         P>|t|        2.5 %      97.5 %
e401  3355.420045  362.414861  9.258506  2.073132e-20  2645.099969  4065.74012

Run the DML model: ATTE

• ATTE: Average Treatment Effects of the Treated

# Run the model
dml_model_irm_ATTE = dml.DoubleMLIRM(df_dml, g_0, m_0, score='ATTE')
# Output the model's findings
print(dml_model_irm_ATTE.fit())

================== DoubleMLIRM Object ==================

------------------ Data summary      ------------------
Outcome variable: net_tfa
Treatment variable(s): ['e401']
Covariates: ['nifa', 'tw', 'age', 'inc', 'fsize', 'educ', 'db', 'marr', 'twoearn', 'p401', 'pira', 'hown']
Instrument variable(s): None
No. Observations: 9915

------------------ Score & algorithm ------------------
Score function: ATTE

------------------ Machine learner   ------------------
Learner ml_g: GradientBoostingRegressor(learning_rate=0.01, max_depth=2, n_estimators=1000,
                          subsample=0.5)
Learner ml_m: GradientBoostingClassifier(learning_rate=0.01, loss='exponential', max_depth=2,
                           n_estimators=1000, subsample=0.5)
Out-of-sample Performance:
Regression:
Learner ml_g0 RMSE: [[11009.26883587]]
Learner ml_g1 RMSE: [[20395.48301999]]
Classification:
Learner ml_m Log Loss: [[0.28201358]]

------------------ Resampling        ------------------
No. folds: 5
No. repeated sample splits: 1

------------------ Fit summary       ------------------
              coef     std err          t          P>|t|        2.5 %  e401  10103.925289  400.757188  25.212087  2.952209e-140  9318.455633   

Other twists on the model

1. Change the machine learning backend
   • Our models used dml2
   • You can switch to dml1 using dml_procedure='dml1'
     • Note that this requires using an older package version, as dml1 is deprecated
   • dml1 follows the math in these slides
     • Solve for a condition equal to zero for each model, and then average the estimators
     • dml2 solves the for the average of the condition being equal to zero overall
2. Run multiple iterations of the model
   • The paper uses 100 iterations, emulate this by adding n_rep=100
3. Change the machine learning models fed to the DoubleML model
   • An example of using “Histogram-based Gradient Boosting” is in the Jupyter notebook
     • This is a much faster GBM-like model

DoubleML DiD

An example of DML DiD

• On the next slides, I present the code for the DML DiD test from Crowley et al. 2023
• The code uses an estimator that is doubly robust following Sant’Anna and Zhao (2020)
• The code uses the DiD approach from Chang (2020)

Pros

• DiD with DoubleML!

Cons

• Hard to use many Fixed Effects
• No model for clustered standard errors (yet)
• Only works with basic DiD estimation

Initializing the model in Python

Here we define \(m_0\) in ml_m and \(g_0\) in ml_g. This paper uses XGBoost as the underlying algorithm, so the model is fully nonparametric.

# Set up the ML functions to use within the DoubleML routine
# Initial
ml_g = ml_l_xgb = XGBRegressor(objective = "reg:squarederror", eta = 0.1,
                        n_estimators = 100)
ml_m = XGBClassifier(use_label_encoder = False,
                        objective = "binary:logistic",
                        eval_metric = "logloss",
                        eta = 0.1, n_estimators = 100)

# Set a random seed for replicability
np.random.seed(65474)

Setting up the data and estimator

# Build the data object
data = df[['Avg_Fake', 'treat', 'post_reg', 'Size_w', 'ROA_w', 'Market_to_Book_w',
           'Leverage_w', 'SOE_indicator', 'dailyret_w', 'retvol_qtr_w',
           'disclosure_rating', 'avg_daily_overall_tone_w', 'log_avg_total_words',
           'log_daily_post_count']]
obj_dml_data = dml.DoubleMLData(data, 'Avg_Fake', 'treat', t_col='post_reg')

# Compile the DMLDiD estimator
dml_did_obj = dml.DoubleMLDIDCS(obj_dml_data, ml_g, ml_m)

# Run and output the results
output = dml_did_obj.fit()  # Takes 7 minutes using 16 cores

print(output)  # shown on the next slide

• The above code is very similar to our base DoubleML model
  • Just note the use of dml.DoubleMLDIDCS() instead of dml.DoubleMLIRM()

Model output

================== DoubleMLDIDCS Object ==================

------------------ Data summary      ------------------
Outcome variable: Avg_Fake
Treatment variable(s): ['treat']
Covariates: ['Size_w', 'ROA_w', 'Market_to_Book_w', 'Leverage_w', 'SOE_indicator', 'dailyret_w', 'retvol_qtr_w', 'disclosure_rating', 'avg_daily_overall_tone_w', 'log_avg_total_words', 'log_daily_post_count']
Instrument variable(s): None
Time variable: post_reg
No. Observations: 2633704

------------------ Score & algorithm ------------------
Score function: observational
DML algorithm: dml2

------------------ Machine learner   ------------------
Learner ml_g: XGBRegressor(base_score=None, booster=None, callbacks=None,
             colsample_bylevel=None, colsample_bynode=None,
             colsample_bytree=None, early_stopping_rounds=None,
             enable_categorical=False, eta=0.1, eval_metric=None,
             feature_types=None, gamma=None, gpu_id=None, grow_policy=None,
             importance_type=None, interaction_constraints=None,
             learning_rate=None, max_bin=None, max_cat_threshold=None,
             max_cat_to_onehot=None, max_delta_step=None, max_depth=None,
             max_leaves=None, min_child_weight=None, missing=nan,
             monotone_constraints=None, n_estimators=100, n_jobs=None,
             num_parallel_tree=None, predictor=None, ...)
Learner ml_m: XGBClassifier(base_score=None, booster=None, callbacks=None,
              colsample_bylevel=None, colsample_bynode=None,
              colsample_bytree=None, early_stopping_rounds=None,
              enable_categorical=False, eta=0.1, eval_metric='logloss',
              feature_types=None, gamma=None, gpu_id=None, grow_policy=None,
              importance_type=None, interaction_constraints=None,
              learning_rate=None, max_bin=None, max_cat_threshold=None,
              max_cat_to_onehot=None, max_delta_step=None, max_depth=None,
              max_leaves=None, min_child_weight=None, missing=nan,
              monotone_constraints=None, n_estimators=100, n_jobs=None,
              num_parallel_tree=None, predictor=None, ...)
Out-of-sample Performance:
Learner ml_g_d0_t0 RMSE: [[0.09086587]]
Learner ml_g_d0_t1 RMSE: [[0.08300815]]
Learner ml_g_d1_t0 RMSE: [[0.09469256]]
Learner ml_g_d1_t1 RMSE: [[0.08113207]]
Learner ml_m RMSE: [[0.47046926]]

------------------ Resampling        ------------------
No. folds: 5
No. repeated sample splits: 1
Apply cross-fitting: True

------------------ Fit summary       ------------------
           coef   std err          t          P>|t|     2.5 %    97.5 %
treat -0.007622  0.000211 -36.155816  3.014536e-286 -0.008035 -0.007209

DoubleML with clustered standard errors

An example of DML with clustered SEs

• On the next slides, I present the code for the clustered standard error DoubleML model from Crowley, Lou, Tan, and Zhang (2025)
• The model is from Chiang et al. (2022)
  • Bootstrapping and sampling over the clustering variables

Pros

• Clustering with DoubleML!
  • One-way
  • Two-way

Cons

• Hard to use many Fixed Effects
• Only works with basic treatment effect estimation
• High computational cost and data requirements

Initializing the model in R

# Required imports
library('hdm')
library('DoubleML')
library('mlr3')
library('mlr3learners')

# Note that all the data is in `df`

# Data setup (building a data.frame for it)
formula = formula(~ -1 + mcap + roa + dt_at + log_followers + log_friends +
                    log_total_tweets + as.factor(year) + event + as.factor(sic1))
data_transf = data.frame(model.matrix(formula, df))

# Set up the learners using glmnet under the hood
lasso_l = lrn("regr.cv_glmnet", nfolds = 10, s = "lambda.min")
lasso_m = lrn("classif.cv_glmnet", nfolds = 10, s = "lambda.min")

• Here we use a LASSO model, so this is only semi-parametric
  • lasso_l is \(g_0\), lasso_m is \(m_0\)

Setting up the data and estimator

# Constructing the DML data object
y_col = 'log_tweets_num'
d_col = 'misinformation'
cluster_cols = c('region', 'year')
dml_df = cbind(df[c(y_col, d_col, cluster_cols)], data_transf)
dml_data = DoubleMLClusterData$new(dml_df,
                                   y_col=y_col,
                                   d_cols=d_col,
                                   cluster_cols=cluster_cols,
                                   x_cols=names(data_transf))

set.seed(1111)
#dml_data$z_cols = z_col
#dml_data$cluster_cols = c('model.id', 'cdid')
n_rep = 10
dml_pl = DoubleMLPLR$new(dml_data, lasso_l, lasso_m,
                         n_folds=2, n_rep=n_rep)
dml_pl$fit()

• Note that the R doubleml package allows for specifying clustering via the cluster_cols argument in the data construction function

Model output

dml_pl$summary()

Estimates and significance testing of the effect of target variables
               Estimate. Std. Error t value Pr(>|t|)  
misinformation   -0.6509     0.2956  -2.202   0.0277 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• This model does not give coefficient estimates

Implementing causal random forests

Some details

1. These models shift their focus to CATE
   • More granular CATEs than DML
   • It can still compute ATE though
2. The models can work on honest trees
   • These use 2 splits of training data: 1 for tree splits and 1 for leaves
   • This reduces overfitting and is mathematically required to compute confidence intervals

Python setup and a simple model

np.random.seed(1234)

# CausalForest expects treatments with shape (n_samples, n_treatments).
# For a binary treatment, encode T numerically as 0/1.
T_grf = np.asarray(T).reshape(-1, 1)

cf = CausalForest(n_estimators=400, min_samples_leaf=20,
    max_depth=12, honest=True, n_jobs=-1, random_state=1234)

# Fit the model
cf.fit(X, T_grf, Y)

# predict() returns one treatment-effect coefficient per treatment.
tau_hat = cf.predict(X).ravel()
ate_hat = tau_hat.mean()

print("Estimated ATE from causal forest:", round(float(ate_hat), 3))

Estimated ATE from causal forest: 7497.59

Treatment effect heterogeneity: By observation

results = df[[y_col, t_col] + x_cols].copy()
results["tau_hat"] = np.asarray(tau_hat).ravel()
results[["tau_hat"]].describe()

       tau_hat
count  9915.00
mean   7497.59
std    8163.20
min    -1254.29
25%    1160.81
50%    4486.01
75%    11754.16
max    37915.73

Treatment effect heterogeneity: By group

subgroup_cols = {"homeowner": "hown", "db_pension": "db", "married": "marr"}
subgroup_results = {}

for subgroup_name, subgroup in subgroup_cols.items():
    summary = (
        results.groupby(subgroup, dropna=False)["tau_hat"]
        .agg(n="size", mean_tau="mean", median_tau="median", std_tau="std")
        .reset_index()
    )
    subgroup_results[subgroup] = summary

    print(f"\nTreatment-effect heterogeneity by {subgroup_name}:")
    print(summary.round(3))

    # Difference between subgroup 1 and subgroup 0
    means = results.groupby(subgroup)["tau_hat"].mean()
    if 0 in means.index and 1 in means.index:
        difference = means.loc[1] - means.loc[0]
        print(f"Difference in mean effects (1 minus 0): {difference:.3f}")

Treatment-effect heterogeneity by homeowner:
   hown     n   mean_tau  median_tau   std_tau
0     0  3617   2375.767     919.405  3805.574
1     1  6298  10439.101    7142.042  8536.744
Difference in mean effects (1 minus 0): 8063.334

Treatment-effect heterogeneity by db_pension:
   db     n  mean_tau  median_tau   std_tau
0   0  7228  7058.466    3999.788  8014.210
1   1  2687  8678.828    5614.549  8439.952
Difference in mean effects (1 minus 0): 1620.362

Treatment-effect heterogeneity by married:
   marr     n  mean_tau  median_tau  std_tau
0     0  3918  5234.574    2567.912  6704.58
1     1  5997  8976.079    5747.977  8676.59
Difference in mean effects (1 minus 0): 3741.505

Pros and Cons of each model type

Feature	DoubleML	Causal.Random.Forest
Architecture	Any ML estimator with both regression and classification	Tree-based, especially honest forests
Output statistic	ATE or group-based CATE	CATE, including individual-level
Applicable to DiD?	Yes, with a dedicated estimator	Yes, via causalfe (Python) or EconML (Python) through the DML hybrid model
Applicable to two-way FEs?	Yes, via xtdml (R) or through de-meaning by group and time	Yes, via causalfe (Python)
Supports clustering?	Yes	Yes via grf (R)

Conclusion

Wrap-up

DoubleML can help in cleanly identifying treatment effects

• Easy to implement
• Many models now available, making it more general-purpose

Causal Random Forest can help in identifying heterogeneous treatment effects

• Easy to implement
• Can be combined with DoubleML

ML and Econometrics are not at odds with one another

• You can use ML to strengthen an econometric framework

ML and econometrics

ML is essentially just another tool in the econometrics toolbox!

Packages used for these slides

Python

• doubleML
• econml
• numpy
• pandas
• sklearn

R

• {doubleml}
• grf
• glmnet
• {hdm}
• kableExtra
• knitr
• mlr3
• mlr3learners
• quarto
• reticulate
• revealjs

References

• Chalfin, Aaron, Oren Danieli, Andrew Hillis, Zubin Jelveh, Michael Luca, Jens Ludwig, and Sendhil Mullainathan. “Productivity and selection of human capital with machine learning.” American Economic Review 106, no. 5 (2016): 124-27.
• Chang, Neng-Chieh. “Double/debiased machine learning for difference-in-differences models.” The Econometrics Journal 23, no. 2 (2020): 177-191.
• Chernozhukov, Victor, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, and Whitney Newey. “Double/debiased/neyman machine learning of treatment effects.” American Economic Review 107, no. 5 (2017): 261-65.
• Chiang, Harold D., Kengo Kato, Yukun Ma, and Yuya Sasaki. “Multiway cluster robust double/debiased machine learning.” Journal of Business & Economic Statistics 40, no. 3 (2022): 1046-1056.
• Crowley, Richard M., Yun Lou, Samuel T. Tan, and Liandong Zhang. “Does Misinformation Regulation Reduce Fake News in Financial Markets? Evidence from East Guba.” Woking paper, Singapore Management University (2023).
• Crowley, Richard M., Yun Lou, Samuel T. Tan, and Liandong Zhang. “Misinformation regulations: Early evidence on corporate social media strategy.” Review of Accounting Studies (2025): 1
• Davis, Jonathan MV, and Sara B. Heller. “Rethinking the benefits of youth employment programs: The heterogeneous effects of summer jobs.” Review of economics and statistics 102, no. 4 (2020): 664-677.
• Deryugina, Tatyana, Garth Heutel, Nolan H. Miller, David Molitor, and Julian Reif. “The mortality and medical costs of air pollution: Evidence from changes in wind direction.” American Economic Review 109, no. 12 (2019): 4178-4219.
• Sant’Anna, Pedro HC, and Jun Zhao. “Doubly robust difference-in-differences estimators.” Journal of Econometrics 219, no. 1 (2020): 101-122.