Reproducing Machine Learning in Economics with Python

“

Author: Shengsheng Li (A Little Ant in Management)Email:[email protected]

“

Source: Jiang Jinhe, Huang Shan. Research on the Impact of New Trade Formats on Green Technology Innovation: Evidence from the Policy of Cross-Border E-Commerce Comprehensive Experimental Zone [J]. Quantitative Economic and Technical Economic Research, 2024, 41(12): 133-154. DOI:10.13653/j.cnki.jqte.20240926.001.

No one supports your ambition, you must tread through the snow to the mountain peak

• (1) Parameter Description
• (2) Code Reproduction
• (3) One-time Generation of Fixed Effects Code Reference

Stata runs machine learning very slowly, not just slow, but extremely slow; running a method takes about several hours, while using Python can complete it in just a few seconds. For reference, Stata’s machine learning can refer to a dual machine learning example for DID. If there are issues running the code, please refer to the previous post on the solution for “Dual Machine Learning ddml not running”.

I do not recommend using Stata for machine learning. Therefore, I referred to the “Python” section to reproduce Jiang Jinhe and Huang Shan’s dual machine learning article in this post, which uses Python to reproduce the machine learning part of Jiang Jinhe, Huang Shan’s research on the impact of new trade formats on green technology innovation: evidence from the policy of cross-border e-commerce comprehensive experimental zone [J]. Quantitative Economic and Technical Economic Research, 2024, 41(12): 133-154. I spent an afternoon exploring it, so there may be some errors. Since I reinstalled Python today (the old laptop was scrapped), I used the latest version of Python, which has avoided some pitfalls, for example: cross_val_predict in sklearn 1.3 and later has stricter checks on the estimator passed in: it requires the estimator to directly implement fit and predict, while StackingRegressor first performs cross-validation predictions on the first layer model and then uses the second layer model (default is RidgeCV) for stacking. At the same time, I provided a one-time generation of fixed effects (entity-centered) to avoid one-hot dimension explosion in Python code.

If you are interested in the code and data in this post, reply with the keyword Machine Learning 1.

Note: The following content is sourced from the official website of Quantitative Economic and Technical Economic Research and is for communication and learning purposes only. If there is any infringement, please contact the above email.

(1) Parameter Description

Parameter Name	Meaning
cv=5	Number of cross-validation folds. The training data is divided into 5 parts (5-fold CV), using 4 parts for training and 1 part for validation, repeating this 5 times to average the validation error, which is used to select the optimal regularization strength α.
random_state=42	Random seed. Cross-validation will shuffle the order when dividing the data, setting a fixed value allows the results to be reproducible (the same division every time you run it). 42 is just a conventional writing, it can be changed to 0 or 2024.
max_iter=5000	Maximum number of iterations. Lasso uses coordinate descent for iterative solving, defaulting to 1000 iterations. If the data volume is large or there are many features, 1000 iterations may not converge before stopping; raising the limit to 5000 can avoid “not converged” warnings or biases.

(2) Code Reproduction

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 14 15:22:40 2025

@author: Shengsheng Li
"""

pip install scikit-learn
pip install matplotlib
pip install EconML

import pandas as pd
import numpy as np
import statsmodels.api as sm
from sklearn.preprocessing import OneHotEncoder, PolynomialFeatures
from sklearn.linear_model import LinearRegression, LassoCV, Lasso, Ridge, RidgeCV, ElasticNet, ElasticNetCV
from sklearn.model_selection import cross_val_predict
from sklearn.ensemble import StackingRegressor, RandomForestRegressor, GradientBoostingRegressor
from sklearn.utils import resample
from sklearn.svm import SVR, LinearSVR
from econml.dml import DML, LinearDML, SparseLinearDML
from matplotlib import pyplot as plt
from sklearn.linear_model import LassoCV
from sklearn.model_selection import cross_val_predict


# 1. Read data
data = pd.read_stata(r'F:/Statalearning/do/2025/did3/Data/data.dta')

# 2. Variable splitting
## Outcome variable
Y = data[['gti']]

## Treatment variable
T = data[['cbecpolicy']]

## Control variables
W = data[['pgdp', 'uis', 'fin', 'inet', 'hca', 'urban', 'ino', 'enr', 'upm', 'gov', 'ure', 'ums']]


## Record individual and time fixed effects, and sum to W
id = data[['id']]
year = data[['year']]
id_oh = OneHotEncoder(sparse_output=False).fit_transform(id)
year_oh = OneHotEncoder(sparse_output=False).fit_transform(year)
W = np.column_stack((W, year_oh, id_oh))


# 3. Estimate the residuals of T

# Estimate the residuals of T
reg_t = LassoCV(cv=5, random_state=42, max_iter=5000)
T_hat = cross_val_predict(reg_t, W, T.values.ravel(), cv=5)
T_residual = T.values.ravel() - T_hat


# 4. Estimate the residuals of Y
# Estimate the residuals of Y
reg_y = LassoCV(cv=5, random_state=42, max_iter=5000)
Y_hat = cross_val_predict(reg_y, W, Y.values.ravel(), cv=5)
Y_residual = Y.values.ravel() - Y_hat


# 5. Perform linear regression with T_residual and Y_residual
X = sm.add_constant(T_residual)  # Add constant term
model_sm = sm.OLS(Y_residual, X).fit()

tau_hat = np.linalg.lstsq(T_residual.reshape(-1, 1), Y_residual, rcond=None)[0][0]
print(f'Hand-calculated PLM estimate τ = {tau_hat:.4f}')


# 6. statsmodels OLS residual regression (with standard errors)
ols = sm.OLS(Y_residual, sm.add_constant(T_residual)).fit()
print(ols.summary().tables[1])

# Reference
# Step 4: Output coefficients and standard errors
coef = model_sm.params[1]  # Coefficient of T_residual
std_err = model_sm.bse[1]  # Standard error of T_residual

# Step 5: Calculate 95% confidence interval
conf_int = model_sm.conf_int(alpha=0.05)[1]  # Confidence interval of T_residual

print(f"Coefficient for T: {coef}")
print(f"Standard Error: {std_err}")
print(f"95% Confidence Interval for T: [{conf_int[0]}, {conf_int[1]}]")


# 7. Compare with econml's LinearDML (DoubleML)
dml = LinearDML(model_y=LassoCV(cv=5, random_state=42, max_iter=5000),
                model_t=LassoCV(cv=5, random_state=42, max_iter=5000),
                cv=5, random_state=42)
dml.fit(y, T, W=W)
print('\nLinearDML Average Treatment Effect')
print(dml.summary())


# Code reproduction

# 1. Read data
data = pd.read_stata(r'F:/Statalearning/do/2025/did3/Data/data.dta')

# 2. Variable splitting
## Outcome variable
Y = data[['gti']]

## Treatment variable
T = data[['cbecpolicy']]

## Control variables
W = data[['pgdp', 'uis', 'fin', 'inet', 'hca', 'urban', 'ino', 'enr', 'upm', 'gov', 'ure', 'ums']]


## Record individual and time fixed effects, and sum to W
id = data[['id']]
year = data[['year']]
id_oh = OneHotEncoder(sparse_output=False).fit_transform(id)
year_oh = OneHotEncoder(sparse_output=False).fit_transform(year)
W = np.column_stack((W, year_oh, id_oh))

# 3. lassocv output results, corresponding to Table 1, Column 1
reg_t = LassoCV(cv=5, random_state=42, max_iter=5000)
T_hat = cross_val_predict(reg_t, W, T.values.ravel(), cv=5)
T_residual = T.values.ravel() - T_hat

# Estimate the residuals of Y
reg_y = LassoCV(cv=5, random_state=42, max_iter=5000)
Y_hat = cross_val_predict(reg_y, W, Y.values.ravel(), cv=5)
Y_residual = Y.values.ravel() - Y_hat

# Perform linear regression with T_residual and Y_residual
X = sm.add_constant(T_residual)  # Add constant term
model_sm = sm.OLS(Y_residual, X).fit()

tau_hat = np.linalg.lstsq(T_residual.reshape(-1, 1), Y_residual, rcond=None)[0][0]
print(f'Hand-calculated PLM estimate τ = {tau_hat:.4f}')

# statsmodels OLS residual regression (with standard errors)
ols = sm.OLS(Y_residual, sm.add_constant(T_residual)).fit()
print(ols.summary().tables[1])

# 4. elasticcv output results, corresponding to Table 1, Column 2
reg_t = ElasticNetCV(cv=5, random_state=42, max_iter=5000)
T_hat = cross_val_predict(reg_t, W, T.values.ravel(), cv=5)
T_residual = T.values.ravel() - T_hat

# Estimate the residuals of Y
reg_y = ElasticNetCV(cv=5, random_state=42, max_iter=5000)
Y_hat = cross_val_predict(reg_y, W, Y.values.ravel(), cv=5)
Y_residual = Y.values.ravel() - Y_hat

# Perform linear regression with T_residual and Y_residual
X = sm.add_constant(T_residual)  # Add constant term
model_sm = sm.OLS(Y_residual, X).fit()

tau_hat = np.linalg.lstsq(T_residual.reshape(-1, 1), Y_residual, rcond=None)[0][0]
print(f'Hand-calculated PLM estimate τ = {tau_hat:.4f}')

# statsmodels OLS residual regression (with standard errors)
ols = sm.OLS(Y_residual, sm.add_constant(T_residual)).fit()
print(ols.summary().tables[1])

# 4. elasticcv output results, corresponding to Table 1, Column 3
# Estimate the residuals of T
reg_t = SVR(kernel='rbf', C=1.0, gamma='scale')
T_hat = cross_val_predict(reg_t, W, T.values.ravel(), cv=5)
T_residual = T.values.ravel() - T_hat

# Estimate the residuals of Y
reg_y = SVR(kernel='rbf', C=1.0, gamma='scale')
Y_hat = cross_val_predict(reg_y, W, Y.values.ravel(), cv=5)
Y_residual = Y.values.ravel() - Y_hat

# Perform linear regression with T_residual and Y_residual
X = sm.add_constant(T_residual)  # Add constant term
model_sm = sm.OLS(Y_residual, X).fit()

tau_hat = np.linalg.lstsq(T_residual.reshape(-1, 1), Y_residual, rcond=None)[0][0]
print(f'Hand-calculated PLM estimate τ = {tau_hat:.4f}')

# statsmodels OLS residual regression (with standard errors)
ols = sm.OLS(Y_residual, sm.add_constant(T_residual)).fit()
print(ols.summary().tables[1])


# Corresponding to Table 1, Column 3
Y = data[['gti1']]

reg_t = LassoCV(cv=5, random_state=42, max_iter=5000)
T_hat = cross_val_predict(reg_t, W, T.values.ravel(), cv=5)
T_residual = T.values.ravel() - T_hat

# Estimate the residuals of Y
reg_y = LassoCV(cv=5, random_state=42, max_iter=5000)
Y_hat = cross_val_predict(reg_y, W, Y.values.ravel(), cv=5)
Y_residual = Y.values.ravel() - Y_hat

# Perform linear regression with T_residual and Y_residual
X = sm.add_constant(T_residual)  # Add constant term
model_sm = sm.OLS(Y_residual, X).fit()

tau_hat = np.linalg.lstsq(T_residual.reshape(-1, 1), Y_residual, rcond=None)[0][0]
print(f'Hand-calculated PLM estimate τ = {tau_hat:.4f}')

# statsmodels OLS residual regression (with standard errors)
ols = sm.OLS(Y_residual, sm.add_constant(T_residual)).fit()
print(ols.summary().tables[1])

# Corresponding to Table 1, Column 4
Y = data[['gti2']]

reg_t = LassoCV(cv=5, random_state=42, max_iter=5000)
T_hat = cross_val_predict(reg_t, W, T.values.ravel(), cv=5)
T_residual = T.values.ravel() - T_hat

# Estimate the residuals of Y
reg_y = LassoCV(cv=5, random_state=42, max_iter=5000)
Y_hat = cross_val_predict(reg_y, W, Y.values.ravel(), cv=5)
Y_residual = Y.values.ravel() - Y_hat

# Perform linear regression with T_residual and Y_residual
X = sm.add_constant(T_residual)  # Add constant term
model_sm = sm.OLS(Y_residual, X).fit()

tau_hat = np.linalg.lstsq(T_residual.reshape(-1, 1), Y_residual, rcond=None)[0][0]
print(f'Hand-calculated PLM estimate τ = {tau_hat:.4f}')

# statsmodels OLS residual regression (with standard errors)
ols = sm.OLS(Y_residual, sm.add_constant(T_residual)).fit()
print(ols.summary().tables[1])

(3) One-time Generation of Fixed Effects Code Reference

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 14 16:36:15 2025

@author: Shengsheng Li
"""

import pandas as pd
import numpy as np
import statsmodels.api as sm
from sklearn.linear_model import LassoCV
from sklearn.model_selection import cross_val_predict
from econml.dml import LinearDML

# 1. Read data
data = pd.read_stata(r'F:/Statalearning/do/2025/did3/Data/data.dta')

# 2. Variable splitting
y = data['gti'].values
T = data['cbecpolicy'].values
X = data[['pgdp', 'uis', 'fin', 'inet', 'hca',
          'urban', 'ino', 'enr', 'upm', 'gov', 'ure', 'ums']]

# 3. Entity-time centering (i.e., removing one-hot to solve dimension explosion)
id_mean = data.groupby('id')[['gti', 'cbecpolicy'] + X.columns.tolist()].transform('mean')
year_mean = data.groupby('year')[['gti', 'cbecpolicy'] + X.columns.tolist()].transform('mean')
overall_mean = data[['gti', 'cbecpolicy'] + X.columns.tolist()].mean()

y_demean = y - id_mean['gti'].values - year_mean['gti'].values + overall_mean['gti']
T_demean = T - id_mean['cbecpolicy'].values - year_mean['cbecpolicy'].values + overall_mean['cbecpolicy']
X_demean = X - id_mean[X.columns] - year_mean[X.columns] + overall_mean[X.columns]

W = X_demean.values          # Already includes entity/time fixed effects

# 4. Manual Partially Linear Model (PLM)
reg_t = LassoCV(cv=5, random_state=42, max_iter=5000)
reg_y = LassoCV(cv=5, random_state=42, max_iter=5000)

T_hat   = cross_val_predict(reg_t, W, T_demean, cv=5)
T_resid = T_demean - T_hat

Y_hat   = cross_val_predict(reg_y, W, y_demean, cv=5)
Y_resid = y_demean - Y_hat

tau_hat = np.linalg.lstsq(T_resid.reshape(-1, 1), Y_resid, rcond=None)[0][0]
print(f'Hand-calculated PLM estimate τ = {tau_hat:.4f}')

# 5. statsmodels OLS residual regression (with standard errors)
ols = sm.OLS(Y_resid, sm.add_constant(T_resid)).fit()
print(ols.summary().tables[1])

# 6. Compare with econml's LinearDML (DoubleML)
dml = LinearDML(model_y=LassoCV(cv=5, random_state=42, max_iter=5000),
                model_t=LassoCV(cv=5, random_state=42, max_iter=5000),
                cv=5, random_state=42)
dml.fit(y, T, W=W)
print('\nLinearDML Average Treatment Effect')
print(dml.summary())