Session 7 · Four Papers, One Toolkit

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
import statsmodels.api as sm
from statsmodels.miscmodels.ordinal_model import OrderedModel

# The data: one CSV per paper on the course site, read over the internet
DATA = "https://warint.github.io/quantitative-methods/data/session-07/"

# Step 1 · Load the data: one row is a country-day
d = pd.read_csv(DATA + "paper1_fraiberger_2021.csv")
d = d.dropna(subset=["F0_index", "L1_sent1_best", "L1_vix"])
d["up"] = (d["F0_index"] > 0).astype(int)            # 1 if the market rose that day
print(f"{len(d):,} country-days, {d['iso'].nunique()} countries, share of up days {d['up'].mean():.3f}")

# Step 2 · Three lines of Python (session 03): today's return on yesterday's news
# The formula reads "F0_index explained by L1_sent1_best". statsmodels
# adds the intercept for you; write ~ 0 + x if you ever want it suppressed,
# and know that R-squared stops being interpretable when you do.
model = smf.ols("F0_index ~ L1_sent1_best", data=d).fit()

# .params holds the two coefficients. The slope is in percentage points of
# return per one standard deviation of sentiment — always say the units.
print(model.params)

# Step 3 · How much did the line explain? (session 03)
# R-squared is the share of variance the line accounts for — and nothing about
# whether the line is the right shape. The rest of this hour is that difference.
print(f"R-squared = {model.rsquared:.3f}")

# Step 4 · Two predictors: what "controlling for" means (session 03)
# Adding a second predictor changes what the first coefficient means: it is now
# the association holding global risk (yesterday's VIX) fixed, not the raw
# association. The two answer different questions, and the number will move.
model2 = smf.ols("F0_index ~ L1_sent1_best + L1_vix", data=d).fit()
print(model2.params)

# Step 5 · The residual plot is the instrument (session 03)
# Residuals vs Fitted: is there structure the straight line missed?
plt.scatter(model2.fittedvalues, model2.resid, s=12, alpha=0.6)
plt.axhline(0, color="black", lw=1)
plt.xlabel("Fitted values"); plt.ylabel("Residuals")
plt.show()

# Step 6 · Binary logit: log-odds (session 04, LR 6) — does the market go up today?
# Fit the logistic regression. The coefficients come out as log-odds.
m2 = smf.logit("up ~ L1_sent1_best", data=d).fit(disp=0)
print(m2.params.round(4))

# Step 7 · Binary logit: odds (session 04, LR 7)
# Log-odds -> odds ratio: take the exponential of each coefficient
print(np.exp(m2.params["L1_sent1_best"]).round(4))

# Step 8 · Binary logit: predicted probability (session 04, LR 8)
# predict() does the same calculation for you: a bad, an ordinary and a good news day
print(m2.predict(pd.DataFrame({"L1_sent1_best": [-1, 0, 1]})).round(3))

# Step 9 · Ordinal logit (session 04, OLR 3) — the day's return in thirds: down < flat < up
d["day"] = pd.qcut(d["F0_index"], 3, labels=["down", "flat", "up"])
# Keep the order: down < flat < up
order = ["down", "flat", "up"]
d["day_c"] = pd.Categorical(d["day"], categories=order, ordered=True)
om = OrderedModel(d["day_c"], d[["L1_sent1_best", "L1_vix"]], distr="logit").fit(method="bfgs", disp=0)
print(om.summary())

# Step 10 · Ordinal logit: log-odds and odds (session 04, OLR 4)
# Log-odds: the slopes as fitted (the cutpoints are left out)
print(om.params[["L1_sent1_best", "L1_vix"]].round(3))
# Proportional-odds ratios
print(np.exp(om.params[["L1_sent1_best", "L1_vix"]]).round(3))

# Step 11 · Ordinal logit: probabilities (session 04, OLR 5)
# A bad and a good news day, with the VIX at 20. Columns: down, flat, up.
new = pd.DataFrame({"L1_sent1_best": [-1, 1], "L1_vix": [20, 20]})
print(pd.DataFrame(np.asarray(om.predict(new)), columns=order, index=["bad news", "good news"]).round(3))

# Step 12 · Multinomial logit (session 04, MLR 3) — the same three days, without their order
# Code 0 (flat) is the baseline: every coefficient reads "relative to a flat day"
levels = ["flat", "down", "up"]
d["day_m"] = pd.Categorical(d["day"], categories=levels).codes
X = d[["L1_sent1_best"]].astype(float)
X["L1_vix"] = d["L1_vix"].astype(float)
X = sm.add_constant(X)
# Two sets of coefficients at once, each against the flat-day baseline
mn = sm.MNLogit(d["day_m"], X).fit(disp=0)
print(mn.summary().tables[1])          # the coefficient table

# Step 13 · Multinomial logit: log-odds and relative risk ratios (session 04, MLR 4)
# Log-odds and RRR, for down vs flat and up vs flat
print(pd.concat({"log-odds": mn.params, "RRR": np.exp(mn.params)}, axis=1)
      .set_axis(["log-odds: down", "log-odds: up", "RRR: down", "RRR: up"], axis=1).round(3))

# Step 14 · Multinomial logit: probabilities (session 04, MLR 5)
# The VIX at 20, across news days from bad to good
grid = pd.DataFrame({"const": 1.0, "L1_sent1_best": [-2, -1, 0, 1, 2], "L1_vix": 20.0})
out = pd.DataFrame(mn.predict(grid).values, columns=levels, index=grid["L1_sent1_best"])
print(out.round(3))

import warnings
warnings.filterwarnings("ignore")            # convergence notes are not part of the answer

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.miscmodels.ordinal_model import OrderedModel

pd.set_option("display.width", 130)       # as in the session-04 slides, a little wider
pd.set_option("display.max_columns", 10)

# The data: one CSV per paper on the course site, read over the internet
DATA = "https://warint.github.io/quantitative-methods/data/session-07/"

# Step 1 · Load the groups, and what happened to each one next week
FILE = DATA + "paper2_saganowski_2019.csv"     # the authors' file, unchanged
raw = pd.read_csv(FILE)
print(raw.shape)
print(raw["event_type"].value_counts(normalize=True).round(3))   # the base rates

# Step 2 · Three features a person can read
groups = pd.DataFrame({
    "event": raw["event_type"],
    "log2_size": np.log2(raw["g1_size"]),            # +1 = the group is twice as big
    "avg_degree": raw["g1_avg_group_degree_total"],  # ties per member, inside the group
    "aging": raw["g1_IlhanAging"]})                  # how long its members have been around
print(groups[["log2_size", "avg_degree", "aging"]].median().round(2))

# Step 3 · Binary logit: log-odds (session 04, LR 6)
groups["dissolve"] = (groups["event"] == "dissolving").astype(int)
# Fit the logistic regression. The coefficients come out as log-odds.
m2 = smf.logit("dissolve ~ log2_size", data=groups).fit(disp=0)
print(m2.params.round(4))

# Step 4 · Binary logit: odds (session 04, LR 7)
# Log-odds -> odds ratio: take the exponential of each coefficient
print(np.exp(m2.params["log2_size"]).round(4))

# Step 5 · Binary logit: predicted probability (session 04, LR 8)
# predict() does the same calculation for you: groups of 4, 8 and 16 members
print(m2.predict(pd.DataFrame({"log2_size": [2, 3, 4]})).round(3))

# Step 6 · Ordinal logit: the code (session 04, OLR 3)
# Keep the order: shrinking < continuing < growing
order = ["shrinking", "continuing", "growing"]
dat = groups[groups["event"].isin(order)].copy()
dat["event_c"] = pd.Categorical(dat["event"], categories=order, ordered=True)
om = OrderedModel(dat["event_c"], dat[["log2_size", "avg_degree", "aging"]], distr="logit").fit(method="bfgs", disp=0)
print(om.summary())

# Step 7 · Ordinal logit: log-odds -> odds (session 04, OLR 4)
# Log-odds: the slopes as fitted (the cutpoints are left out)
print(om.params[["log2_size", "avg_degree", "aging"]].round(3))
# Proportional-odds ratios
print(np.exp(om.params[["log2_size", "avg_degree", "aging"]]).round(3))

# Step 8 · Ordinal logit: probabilities (session 04, OLR 5)
# A 4-member group, and an otherwise identical 16-member one.
# Columns: shrinking, continuing, growing.
new = pd.DataFrame({"log2_size": [2, 4], "avg_degree": [2.7, 2.7], "aging": [8.5, 8.5]})
print(pd.DataFrame(np.asarray(om.predict(new)), columns=order, index=["4 members", "16 members"]).round(3))

# Step 10 · Multinomial logit: log-odds -> relative risk ratios (session 04, MLR 4)
# Log-odds and RRR, for shrinking, growing and continuing vs dissolving
print(pd.concat({"log-odds": mn.params, "RRR": np.exp(mn.params)}, axis=1)
      .set_axis(["log-odds: shrinking", "log-odds: growing", "log-odds: continuing",
                 "RRR: shrinking", "RRR: growing", "RRR: continuing"], axis=1).round(3))

# Step 11 · Multinomial logit: probabilities (session 04, MLR 5)
# Typical groups (median degree and aging), from 4 to 32 members
grid = pd.DataFrame({"const": 1.0, "log2_size": [2, 3, 4, 5],
                     "avg_degree": 2.7, "aging": 8.5})
out = pd.DataFrame(mn.predict(grid).values, columns=levels, index=grid["log2_size"])
print(out.round(3))

import warnings
warnings.filterwarnings("ignore")       # convergence notes are not part of the answer
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

# The data: one CSV per paper on the course site, read over the internet
DATA = "https://warint.github.io/quantitative-methods/data/session-07/"

# Step 1 · Load the data: one row per country pair, log FDI, 56 candidates
fdi = pd.read_csv(DATA + "paper3_blonigen_piger_2014.csv")
YCOL = "log_fdi"
XCOLS = [c for c in fdi.columns if c not in ("parent", "host", "year", "fdi_stock_musd")]
d = fdi[fdi.fdi_stock_musd > 0].dropna(subset=XCOLS).copy()   # positive FDI, complete rows
d["log_fdi"] = np.log(d["fdi_stock_musd"])
for v in ["parent_gdp", "host_gdp", "distance"]:                # gravity variables in logs
    d[v] = np.log(d[v])
print(f"country pairs: {len(d)}   candidate predictors: {len(XCOLS)}")

# Step 2 · OLS: three lines of Python (session 3)
# The formula reads "log_fdi explained by host_gdp" (both in logs).
model = smf.ols("log_fdi ~ host_gdp", data=d).fit()

# .params holds the two coefficients. The slope is an elasticity: % more FDI
# per 1% more host GDP — always say the units.
print(model.params)

# Step 3 · How much did the line explain?
# R-squared is the share of variance the line accounts for — and nothing about
# whether the line is the right shape.
print(f"R-squared = {model.rsquared:.3f}")

# Step 4 · Two predictors: what "controlling for" means
# Adding distance changes what the host-GDP coefficient means: it is now the
# association holding distance fixed, not the raw association.
model2 = smf.ols("log_fdi ~ host_gdp + distance", data=d).fit()
print(model2.params)
print(f"R-squared = {model2.rsquared:.3f}")

# Step 5 · Ridge (session 5)
from sklearn.linear_model import RidgeCV

Xc, yc = d[XCOLS].astype(float), d[YCOL].astype(float)
ridge = make_pipeline(StandardScaler(),                        # same scale for all
                      RidgeCV(alphas=np.logspace(-3, 3, 60))   # try 60 values of λ
                      ).fit(Xc, yc)
print(f"λ = {ridge[-1].alpha_:.2f}")
print(f"ridge kept {np.sum(ridge[-1].coef_ != 0)} of {len(XCOLS)} — it never drops a predictor")

# Step 6 · Lasso
from sklearn.linear_model import LassoCV

lasso = make_pipeline(StandardScaler(),
                      LassoCV(cv=5, random_state=0, max_iter=50_000)   # λ by 5-fold CV
                      ).fit(Xc, yc)
kept = [c for c, b in zip(XCOLS, lasso[-1].coef_) if b != 0]
print(f"λ = {lasso[-1].alpha_:.4f}")
print(f"lasso kept {len(kept)} of {len(XCOLS)}")

# Step 7 · Elastic net
from sklearn.linear_model import ElasticNetCV

enet = make_pipeline(StandardScaler(),
                     ElasticNetCV(l1_ratio=0.5,            # half lasso, half ridge
                                  cv=5, random_state=0, max_iter=50_000)
                     ).fit(Xc, yc)
kept = [c for c, b in zip(XCOLS, enet[-1].coef_) if b != 0]
print(f"λ = {enet[-1].alpha_:.4f}")
print(f"elastic net kept {len(kept)} of {len(XCOLS)}")

# Step 8 · Which one predicts best? Fit on 70% of the pairs, score on 30%
Xtr, Xte, ytr, yte = train_test_split(d[XCOLS], d[YCOL], test_size=0.3, random_state=7)

def test_rmse(model):                  # fit on 70% of rows, score on the other 30%
    pipe = make_pipeline(StandardScaler(), model).fit(Xtr, ytr)   # standardise first
    kept = int(np.sum(np.abs(pipe[-1].coef_) > 1e-8))             # non-zero coefficients
    return kept, mean_squared_error(yte, pipe.predict(Xte)) ** 0.5

results = {
    "OLS, all predictors": test_rmse(LinearRegression()),
    "ridge":       test_rmse(RidgeCV(alphas=np.logspace(-3, 3, 60))),   # λ by CV
    "lasso":       test_rmse(LassoCV(cv=5, random_state=0, max_iter=50_000)),
    "elastic net": test_rmse(ElasticNetCV(l1_ratio=0.5, cv=5, random_state=0, max_iter=50_000)),
}
res = pd.DataFrame(results, index=["predictors kept", "test RMSE"]).T.astype({"predictors kept": int})
print(res.round(3))

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
import statsmodels.api as sm
from statsmodels.miscmodels.ordinal_model import OrderedModel
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

# The data: one CSV per paper on the course site, read over the internet
DATA = "https://warint.github.io/quantitative-methods/data/session-07/"

# Step 1 · Load the firms: one row is one firm in one year, 1989-1996
firms = pd.read_csv(DATA + "paper4_topalova_khandelwal_2011.csv")
cols = ["tfp1000", "lagtariff", "laginputariff", "prgroup", "govown",
        "foreign", "medium", "small", "age"]
d = firms[firms["yr"] < 1997].dropna(subset=cols)       # no gaps in what we use
print(f"{len(d):,} firm-years, {d['companyname'].nunique():,} firms")
print(d[["tfp1000", "lagtariff", "laginputariff"]].describe().round(3))

# Step 2 · Three lines of Python (session 3): productivity on the output tariff
# The formula reads "tfp1000 explained by lagtariff". statsmodels adds the
# intercept for you. The slope is in log points of productivity per one unit
# of tariff (1.0 = 100 percentage points) — always say the units.
model = smf.ols("tfp1000 ~ lagtariff", data=d).fit()
print(model.params)

# Step 3 · How much did the line explain? (session 3)
# R-squared is the share of variance the line accounts for — and nothing about
# whether the line is the right shape.
print(f"R-squared = {model.rsquared:.3f}")

# Step 4 · Two predictors: what "controlling for" means (session 3)
# Adding a second predictor changes what the first coefficient means: it is now
# the association holding the input tariff fixed, not the raw association.
model2 = smf.ols("tfp1000 ~ lagtariff + laginputariff", data=d).fit()
print(model2.params)
print(f"R-squared = {model2.rsquared:.3f}")

# Step 5 · Residuals vs fitted (session 3): structure the straight line missed?
fit = model2
plt.scatter(fit.fittedvalues, fit.resid, s=12, alpha=0.6)
plt.axhline(0, color="black", lw=1)
plt.xlabel("Fitted values"); plt.ylabel("Residuals")
plt.savefig("p4_resid_vs_fitted.png", dpi=110, bbox_inches="tight")   # saved in your folder
plt.show()

# Step 6 · Binary logit: log-odds (session 4, LR 6) — is the firm more productive than the median firm?
d["high_tfp"] = (d["tfp1000"] > d["tfp1000"].median()).astype(int)   # 1 = above the median
# Fit the logistic regression. The coefficients come out as log-odds.
m2 = smf.logit("high_tfp ~ lagtariff", data=d).fit(disp=0)
print(m2.params.round(4))

# Step 7 · Binary logit: odds (session 4, LR 7)
# Log-odds -> odds ratio: take the exponential of each coefficient
print(np.exp(m2.params["lagtariff"]).round(4))

# Step 8 · Binary logit: predicted probability (session 4, LR 8)
# predict() does the same calculation for you
print(m2.predict(pd.DataFrame({"lagtariff": [0.4, 0.7, 1.0]})).round(3))   # inside the data: 0.1 to 3.5

# Step 9 · Ordinal logit (session 4, OLR 3) — productivity in thirds: low < middle < high
d["tfp_third"] = pd.qcut(d["tfp1000"], 3, labels=["low", "middle", "high"])
# Keep the order: low < middle < high
order = ["low", "middle", "high"]
d["tfp_c"] = pd.Categorical(d["tfp_third"], categories=order, ordered=True)
om = OrderedModel(d["tfp_c"], d[["lagtariff", "laginputariff", "age"]], distr="logit").fit(method="bfgs", disp=0)
print(om.summary())

# Step 10 · Ordinal logit: log-odds and odds (session 4, OLR 4)
# Log-odds: the slopes as fitted (the cutpoints are left out)
print(om.params[["lagtariff", "laginputariff", "age"]].round(3))
# Proportional-odds ratios
print(np.exp(om.params[["lagtariff", "laginputariff", "age"]]).round(3))

# Step 11 · Ordinal logit: probabilities (session 4, OLR 5)
# The same firm (20 years old, input tariff 0.3) under a high and a low output tariff.
# Columns: low, middle, high productivity.
new = pd.DataFrame({"lagtariff": [1.0, 0.4], "laginputariff": [0.3, 0.3], "age": [20, 20]})
print(pd.DataFrame(np.asarray(om.predict(new)), columns=order, index=["tariff 1.0", "tariff 0.4"]).round(3))

# Step 13 · Multinomial logit: log-odds and relative risk ratios (session 4, MLR 4)
pd.set_option("display.width", 130)       # six columns: print them side by side
pd.set_option("display.max_columns", 10)
# Log-odds and RRR, for each owner type vs stand-alone
print(pd.concat({"log-odds": mn.params, "RRR": np.exp(mn.params)}, axis=1)
      .set_axis(["log-odds: group", "log-odds: state", "log-odds: foreign",
                 "RRR: group", "RRR: state", "RRR: foreign"], axis=1).round(3))

# Step 14 · Multinomial logit: probabilities (session 4, MLR 5)
# Large firms: the median productivity, then across the range
grid = pd.DataFrame({"const": 1.0, "medium": 0.0, "small": 0.0,
                     "tfp1000": [d["tfp1000"].median(), -0.5, 0.0, 0.5]})
out = pd.DataFrame(mn.predict(grid).values, columns=levels, index=grid["tfp1000"].round(2))
print(out.round(3))

# Step 15 · The candidates for ridge, lasso and elastic net (session 5)
d = d.copy()
d["age2"] = d["age"] ** 2
years = pd.get_dummies(d["yr"], prefix="yr", drop_first=True).astype(float)   # 1989 is the base
d = pd.concat([d, years], axis=1)
YCOL = "tfp1000"
XCOLS = cols[1:] + ["age2"] + list(years.columns)
Xc, yc = d[XCOLS].astype(float), d[YCOL].astype(float)
print(f"{len(XCOLS)} candidates:", XCOLS)

# Step 16 · Ridge (session 5)
from sklearn.linear_model import RidgeCV

ridge = make_pipeline(StandardScaler(),                        # same scale for all
                      RidgeCV(alphas=np.logspace(-3, 3, 60))   # try 60 values of λ
                      ).fit(Xc, yc)
print(f"λ = {ridge[-1].alpha_:.2f}")
print(f"ridge kept {np.sum(ridge[-1].coef_ != 0)} of {len(XCOLS)} — it never drops a predictor")

# Step 17 · Lasso (session 5)
from sklearn.linear_model import LassoCV

lasso = make_pipeline(StandardScaler(),
                      LassoCV(cv=5, random_state=0, max_iter=50_000)   # λ by 5-fold CV
                      ).fit(Xc, yc)
kept = [c for c, b in zip(XCOLS, lasso[-1].coef_) if b != 0]
print(f"λ = {lasso[-1].alpha_:.4f}")
print(f"lasso kept {len(kept)} of {len(XCOLS)}:", kept)

# Step 18 · Elastic net (session 5)
from sklearn.linear_model import ElasticNetCV

enet = make_pipeline(StandardScaler(),
                     ElasticNetCV(l1_ratio=0.5,            # half lasso, half ridge
                                  cv=5, random_state=0, max_iter=50_000)
                     ).fit(Xc, yc)
kept = [c for c, b in zip(XCOLS, enet[-1].coef_) if b != 0]
print(f"λ = {enet[-1].alpha_:.4f}")
print(f"elastic net kept {len(kept)} of {len(XCOLS)}:", kept)

# Step 19 · The four, on rows they never saw (session 5, as reviewed in session 6)
Xtr, Xte, ytr, yte = train_test_split(Xc, yc, test_size=0.3, random_state=7)

def test_rmse(model):                  # fit on 70% of rows, score on the other 30%
    pipe = make_pipeline(StandardScaler(), model).fit(Xtr, ytr)   # standardise first
    kept = int(np.sum(np.abs(pipe[-1].coef_) > 1e-8))             # non-zero coefficients
    return kept, mean_squared_error(yte, pipe.predict(Xte)) ** 0.5

results = {
    "OLS, all predictors": test_rmse(LinearRegression()),
    "ridge":       test_rmse(RidgeCV(alphas=np.logspace(-3, 3, 60))),   # λ by CV
    "lasso":       test_rmse(LassoCV(cv=5, random_state=0, max_iter=50_000)),
    "elastic net": test_rmse(ElasticNetCV(l1_ratio=0.5, cv=5, random_state=0, max_iter=50_000)),
}
res = pd.DataFrame(results, index=["predictors kept", "test RMSE"]).T.astype({"predictors kept": int})
print(res.round(3))