AKM variance decomposition¶
This is the application the library was written for. Having estimated worker and firm effects, you want to know how much of the variance of log earnings each one explains, and how much comes from the covariance between them: whether good workers sort into good firms.
Fitted values split additively into the covariate part and one term per fixed-effect dimension,
so the variance of \(y\) is the sum of those five variances plus twice every
pairwise covariance between them. Each term is an aggregate over rows, and
resid() hands back a lazy scan of the row-level output, so Polars computes
them in a streaming pass: the numbers come back without the rows ever being
collected.
The fit¶
fit = feols_stream(
"log_earn ~ age_squared + age_cubed | worker_id + firm_id + year",
"data/panel.parquet", workdir="scratch", vcov={"CRV1": "worker_id"},
keep=["region"], # carry region into the row-level file
)
The row-level file holds, per observation, the outcome, the covariate
contribution xb, one column per fixed-effect dimension with that row’s
estimated effect, and the residual, along with the fixed-effect ids and
row_id. keep= carries further columns of the source into it, here region,
for the decomposition by region below.
The decomposition, in one streaming pass¶
import polars as pl
PARTS = ["xb", "fe_worker_id", "fe_firm_id", "fe_year", "resid"]
stats = fit.resid().select(
pl.col("log_earn").var().alias("var_y"),
*[pl.col(c).var().alias(f"var_{c}") for c in PARTS],
pl.cov("fe_worker_id", "fe_firm_id").alias("cov_worker_firm"),
pl.corr("fe_worker_id", "fe_firm_id").alias("corr_worker_firm"),
).collect(engine="streaming").row(0, named=True)
engine="streaming" keeps Polars from materializing the frame. On the
simulated panel:
Variance of log earnings: 0.3953
component variance share
------------------------------------------------------
worker effects 0.1738 44.0%
firm effects 0.0759 19.2%
year effects 0.0002 0.1%
covariates (Xb) 0.0054 1.4%
residual 0.0812 20.5%
2 x covariances between the parts 0.0587 14.8%
------------------------------------------------------
total 0.3953 100.0%
2 x cov(worker, firm) alone: 0.0596 (15.1% of Var(y))
Correlation of worker and firm effects: 0.259
A positive correlation means assortative matching: better workers at better
firms. The same moments by group are a group_by in the same streaming pass;
region is in the file because the fit was asked to keep it. Any column of the
source could be, or the residuals can be joined back onto the source on
row_id instead (see the sample).
Identification¶
Worker and firm effects are only comparable within a connected set of firms,
linked by workers who moved between them. With more than one component, effects
in different components are not comparable. The fit reports how many components
the worker–firm network has (fit.n_components), and fixef("firm_id") gives
each firm’s component, numbered by size so that component 0 is the largest.
Connected sets shows how to find them before fitting, and
how to restrict a fit to the largest one.
The bias in these numbers¶
With few observations per worker, the variance of the worker effects is biased
up and their covariance with the firm effects biased down, by estimation error:
the “limited mobility bias” of Andrews et al. (2008). This example reports the
plug-in decomposition, which is what the estimated effects give you. The number
of observations per worker, 8.5 here, is the one to judge it by, since the bias
grows as it falls. leave_out_kss removes the bias; the
leave-out variance components section explains how.
The full script: examples/akm_variance.py
"""AKM variance decomposition, without ever holding the data in memory.
This is the application the library was written for. Having estimated worker
and firm effects, you want to know how much of the variance of log earnings
each one explains, and how much comes from the covariance between them --
whether good workers sort into good firms.
Fitted values split additively into the covariate part and one term per fixed
effect dimension,
y = Xb + worker + firm + year + residual
so the variance of y is the sum of those five variances plus twice every
pairwise covariance between them. Each term is an aggregate over rows, and
`resid()` hands back a lazy scan of the row-level output, so Polars computes
them in a streaming pass: the numbers come back without the rows ever being
collected.
python examples/akm_variance.py
"""
import polars as pl
from hdfe_stream import feols_stream
from simulated_data import panel_path, workdir
fit = feols_stream(
"log_earn ~ age_squared + age_cubed | worker_id + firm_id + year",
panel_path(),
workdir=workdir("akm"),
vcov={"CRV1": "worker_id"},
# carry region into the row-level file, for the decomposition by region
# below; the fixed-effect id columns are always there
keep=["region"],
verbose=False,
)
# The row-level file holds, per observation: the outcome, the covariate
# contribution `xb`, one column per fixed effect dimension with that row's
# estimated effect, and the residual.
rows = fit.resid()
print("columns available per row:", rows.collect_schema().names())
# One streaming pass over the file computes the whole decomposition.
# engine="streaming" keeps Polars from materializing the frame.
PARTS = ["xb", "fe_worker_id", "fe_firm_id", "fe_year", "resid"]
stats = rows.select(
pl.col("log_earn").var().alias("var_y"),
*[pl.col(c).var().alias(f"var_{c}") for c in PARTS],
pl.cov("fe_worker_id", "fe_firm_id").alias("cov_worker_firm"),
pl.corr("fe_worker_id", "fe_firm_id").alias("corr_worker_firm"),
).collect(engine="streaming").row(0, named=True)
var_y = stats["var_y"]
variances = {c: stats[f"var_{c}"] for c in PARTS}
# everything not in the individual variances is covariance between the parts
covariances = var_y - sum(variances.values())
print(f"\nVariance of log earnings: {var_y:.4f}\n")
print(f"{'component':<34}{'variance':>11}{'share':>9}")
print("-" * 54)
for label, key in [("worker effects", "fe_worker_id"),
("firm effects", "fe_firm_id"),
("year effects", "fe_year"),
("covariates (Xb)", "xb"),
("residual", "resid")]:
value = variances[key]
print(f"{label:<34}{value:>11.4f}{value / var_y:>8.1%}")
print(f"{'2 x covariances between the parts':<34}{covariances:>11.4f}"
f"{covariances / var_y:>8.1%}")
print("-" * 54)
print(f"{'total':<34}{var_y:>11.4f}{1.0:>8.1%}")
# The sorting term on its own: this is the one people quote
sorting = 2 * stats["cov_worker_firm"]
print(f"\n2 x cov(worker, firm) alone: {sorting:.4f} ({sorting / var_y:.1%} of Var(y))")
print(f"\nCorrelation of worker and firm effects: {stats['corr_worker_firm']:.3f}")
print("(positive means assortative matching: better workers at better firms)")
# The same moments by group are a group_by in the same streaming pass. `region`
# is in the file because the fit was asked to keep it; any column of the source
# could be (or join the residuals back onto the source on `row_id` instead;
# see examples/formulas.py).
by_region = (rows.group_by("region")
.agg(pl.len().alias("rows"),
pl.col("log_earn").var().alias("var_y"),
pl.col("fe_firm_id").var().alias("var_firm"),
(2 * pl.cov("fe_worker_id", "fe_firm_id")).alias("sorting"))
.with_columns((pl.col("var_firm") / pl.col("var_y")).alias("firm share"),
(pl.col("sorting") / pl.col("var_y")).alias("sorting share"))
.sort("region")
.collect(engine="streaming"))
print("\nby region:")
print(by_region)
# Identification: worker and firm effects are only comparable within a
# connected set of firms, linked by workers who moved between them. With more
# than one component, effects in different components are not comparable.
print(f"\nconnected components: {fit.n_components}")
print("worker groups that carry identifying information:", f"{fit.n_identifying:,}")
# Firm effects, with the component each firm belongs to
firms = fit.fixef("firm_id").collect()
print("\nfirm effects by component:")
print(firms.group_by("component").agg(
firms=pl.len(),
observations=pl.col("n_obs").sum(),
mean_effect=pl.col("fe_firm_id").mean(),
sd_effect=pl.col("fe_firm_id").std(),
).sort("firms", descending=True))
# A caution worth knowing: with few observations per worker, Var(worker) is
# biased up and Cov(worker, firm) biased down by estimation error -- the
# "limited mobility bias" of Andrews et al. This example reports the raw
# plug-in decomposition, which is what the estimated effects give you.
#
# `leave_out_kss` removes that bias; see examples/leave_out_kss.py. The number
# below is the one to judge it by, since the bias grows as it falls.
print("\nobservations per worker: "
f"{fit.n_obs / fit.n_levels['worker_id']:.1f} on average")