Python vs Stata / R / Excel / SQL / SAS: trade-offs and task-based selection

import pandas as pd
import sqlite3

# in-memory demo: build a tiny "firm panel", then aggregate IN the database
con = sqlite3.connect(":memory:")
pd.DataFrame({
    "province": ["AH", "AH", "AH", "ZJ", "ZJ"],
    "year":     [2010, 2010, 2011, 2010, 2011],
    "employment": [50, 66, 71, 90, 88],
}).to_sql("firm_panel", con, index=False)

query = """
SELECT province, year,
       AVG(employment) AS mean_emp,
       COUNT(*)        AS n_firms
FROM firm_panel
WHERE year BETWEEN 2010 AND 2020
GROUP BY province, year
"""
panel = pd.read_sql(query, con)
print(panel)

Expected output

  province  year  mean_emp  n_firms
0       AH  2010      58.0        2
1       AH  2011      71.0        1
2       ZJ  2010      90.0        1
3       ZJ  2011      88.0        1

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

rng = np.random.default_rng(0)
n = 300
df = pd.DataFrame({
    "firm":  rng.integers(0, 30, n),
    "year":  rng.integers(2010, 2020, n),
    "treat": rng.integers(0, 2, n),
    "post":  rng.integers(0, 2, n),
})
df["y"] = 1 + 1.2 * (df.treat * df.post) + 0.3 * (df.firm % 5) + rng.normal(size=n)

# Stata (one line):    reghdfe y i.treat##i.post, absorb(firm year) cluster(firm)
# R/fixest (one line): feols(y ~ treat*post | firm + year, cluster = ~firm, df)
# Python: spell out fixed effects as factors + cluster-robust SE
m = smf.ols("y ~ treat*post + C(firm) + C(year)", data=df).fit(
    cov_type="cluster", cov_kwds={"groups": df["firm"]})
print(round(float(m.params["treat:post"]), 3))

Expected output

1.214