Utilities

Reporting

hdfe_stream.etable(models, **kwargs)

pf.etable for streaming results: a regression table of HDFEResult and HDFEMulti models, mixed freely with ordinary pyfixest models.

A multi-model result expands into one column per model. The spacing of the fixed-effect side is normalized, so "a + b" and "a+b" share a row (pf.etable itself matches rows on the literal text). With type="tex", labels are escaped for LaTeX.

Parameters:

• models (HDFEResult, HDFEMulti, or a list of them and pyfixest models) – The models, one column each.

• **kwargs – Passed to pf.etable: type, model_stats, labels, felabels, file_name, …

Returns:

object – What pf.etable returns for the type asked for.

Examples

>>> etable([fit1, fit2], type="df", model_stats=["N", "se_type", "r2"])

Multiple testing

hdfe_stream.bonferroni(models, param)

Bonferroni-adjusted p-values for param = 0 across models.

Each p-value, from HDFEResult.tidy, is multiplied by the number of hypotheses (models times parameters), capped at 1. As pyfixest’s bonferroni, which takes one parameter.

Parameters:

• models (HDFEResult, HDFEMulti, or a list of them) – The models.

• param (str or list of str) – The coefficient(s) tested in every model; each must be in every model.

Returns:

polars.DataFrame – One row per model and parameter: model (its position), fml, the tidy columns, and Bonferroni Pr(>|t|).

hdfe_stream.holm(models, param)

Holm’s step-down adjusted p-values for param = 0 across models.

The p-values are sorted from smallest to largest, the j-th of m is multiplied by m - j + 1, and each adjusted p-value is the running maximum of those, capped at 1. Like Bonferroni’s, the adjustment controls the familywise error rate under any dependence between the tests, and no adjusted p-value is larger than Bonferroni’s, so it rejects at least as often (Holm, 1979).

Parameters:

• models – As in bonferroni.

• param – As in bonferroni.

Returns:

polars.DataFrame – As in bonferroni, with Holm Pr(>|t|).

hdfe_stream.rwolf(models, param, reps=9999, seed=None, weights_type='rademacher')

Romano-Wolf adjusted p-values for param = 0 across models, from the wild cluster bootstrap.

Every model is bootstrapped on the same draws (each cluster gets the same weight in every model), and the hypotheses are stepped down from the largest absolute t statistic: each one’s p-value is the share of draws whose largest absolute bootstrap t statistic, over the hypotheses not yet stepped past, reaches its t statistic, with one added to the count and to the draws, and the p-values are then made monotone (Romano and Wolf 2005; Clarke, Romano and Wolf 2020). This controls the familywise error rate while using the dependence between the tests, so it rejects more than Bonferroni or Holm when they are correlated. As pyfixest’s rwolf, with the same draws for the same seed.

Parameters:

• models (HDFEResult, HDFEMulti, or a list of them) – Fit with the same bootcluster=.

• param (str or list of str) – The coefficient(s) tested in every model.

• reps (int, default 9999) – Bootstrap draws.

• seed (int, optional) – Seed for the draws.

• weights_type ({"rademacher", "mammen", "webb", "norm"}, default "rademacher") – The distribution of the cluster weights.

Returns:

polars.DataFrame – As in bonferroni, with RW Pr(>|t|).

hdfe_stream.wyoung(models, param, reps=9999, seed=None, weights_type='rademacher')

Westfall-Young adjusted p-values for param = 0 across models, from the wild cluster bootstrap.

The free step-down resampling method (Westfall and Young 1993): each test’s t statistic, and each draw’s, is turned into a p-value from the t distribution with G - 1 degrees of freedom; in each draw, the p-values are replaced by their running minimum over the tests with larger original p-values; and a test’s adjusted p-value is the share of draws in which that minimum is at most its own p-value, made monotone. As pyfixest’s wyoung, with the same draws for the same seed.

Parameters:

• models – As in rwolf.

• param – As in rwolf.

• reps – As in rwolf.

• seed – As in rwolf.

• weights_type – As in rwolf.

Returns:

polars.DataFrame – As in bonferroni, with WY Pr(>|t|).

Run directories

hdfe_stream.cleanup(workdir=None, older_than_hours=None, force=False, dry_run=False)

Remove leftover run directories, e.g. from killed processes.

Only directories created by hdfe_stream (they carry a marker file) are touched, and empty ones named like them (whose files were deleted but which could not be removed while a file in them was still open, as on NFS); every new fit removes those too.

Runs of this process are removed unless a fit is still running in them (this includes result files of finished fits; their resid() and fixef() stop working). Skipped unless force=True: runs of other processes still alive on this host, runs from other hosts (their liveness can’t be checked), and runs younger than older_than_hours if given.

Parameters:

• workdir (str or Path, optional) – The directory to sweep. Default: the HDFE_STREAM_WORKDIR environment variable if set, else the system temporary directory.

• older_than_hours (float, optional) – Skip runs younger than this.

• force (bool, default False) – Remove runs of live processes and of other hosts too.

• dry_run (bool, default False) – List what would be removed, and remove nothing.

Returns:

list of (str, int) – (path, bytes) for each run that was (or, with dry_run, would be) removed.

Simulated data

Simulated matched employer-employee panels with AKM structure.

Used by the test suite and the examples, and useful for trying the library out before pointing it at real data:

from hdfe_stream import feols_stream
from hdfe_stream.simulate import simulate_akm

simulate_akm(n_workers=50_000).write_parquet("sim.parquet")
feols_stream("log_earn ~ age_squared | worker_id + firm_id + year",
             "sim.parquet").summary()

simulate_akm is the base panel: workers move between firms over a run of years, and log earnings are a worker effect plus a firm effect plus an age profile plus noise. Because the worker and firm effects are drawn explicitly, the true variance decomposition is known, and there is enough worker mobility to connect the firms into one component.

simulate_rich adds the columns the wider feature set needs – categoricals, coarser cluster variables, weights, instruments, a second outcome, a column with missing values – and simulate_trends adds worker-specific time trends for varying-slopes (FEIS) models. simulate_bottleneck builds a panel whose variance components are weakly identified, for leave-out estimation. All of them take a seed and are deterministic.

Only numpy and Polars are used, so nothing here needs the optional dependencies.

hdfe_stream.simulate.simulate_akm(n_workers=50000, n_firms=None, years=range(2005, 2015), p_move=0.08, p_obs=0.85, sd_worker=0.4, sd_firm=0.25, sd_noise=0.3, seed=0, keep_effects=False)

An AKM panel: one row per worker-year actually observed.

Worker and firm effects are positively assorted – better workers tend to be at better firms – so the firm effect is not orthogonal to the worker effect and the fixed effects genuinely have to be estimated jointly.

Parameters:

• n_workers (int, default 50_000) – The number of workers.

• n_firms (int, optional) – The number of firms. Default: n_workers / 15 (at least 50), which keeps firms large enough to be identified.

• years (iterable of int) – The years observed (default 2005-2014).

• p_move (float, default 0.08) – Probability a worker changes firm between consecutive years. Mobility is what connects firms, so lowering it a lot will split the data into several components.

• p_obs (float, default 0.85) – Probability a worker is observed in a given year, so the panel is unbalanced.

• sd_worker (float) – Standard deviations of the worker effect (0.4), the firm effect (0.25) and the residual (0.3).

• sd_firm (float) – Standard deviations of the worker effect (0.4), the firm effect (0.25) and the residual (0.3).

• sd_noise (float) – Standard deviations of the worker effect (0.4), the firm effect (0.25) and the residual (0.3).

• seed (int or anything numpy.random.default_rng accepts, default 0) – The seed; the panel is deterministic given it.

• keep_effects (bool, default False) – Also return the worker and firm effects that generated the data, as true_worker_effect and true_firm_effect. Off by default so they cannot be used as regressors by accident. Useful for checking an estimator against a known answer – in particular for bias corrections to variance components, where the quantity being estimated is a moment of exactly these columns.

Returns:

polars.DataFrame – Shuffled so it is not sorted on disk, with columns worker_id (int64, deliberately not dense), firm_id (string), year (int32), age (float), age_squared, age_cubed (age polynomials scaled into a sane range) and log_earn (the outcome).

Examples

>>> simulate_akm(n_workers=1_000_000).write_parquet("sim.parquet")

hdfe_stream.simulate.simulate_rich(n_workers=50000, n_firms=None, seed=0, extra_seed=9, **kwargs)

simulate_akm plus the columns that exercise the wider feature set.

On top of the base panel, it adds:

• state, region: coarse firm-level groupings, for clustering on something that is not a fixed effect (region is coarser still);

• cohort: year-of-birth decade, another non-fixed-effect cluster variable;

• treat, occ, cat: a binary, an integer-coded and a string categorical, for C(), i() and interaction terms;

• y2: a second outcome, for multiple-estimation formulas;

• x_na: a covariate that is null in about 5% of rows;

• wa, wf: analytic and frequency weights;

• z1, z2, x2: two instruments and the endogenous regressor they predict;

• y_iv: an outcome that depends on x2, for 2SLS.

log_earn also picks up an occupation effect and a post-2010 treatment effect, so i(year, treat) event-study terms have something to find.

Parameters:

• n_workers – As in simulate_akm.

• n_firms – As in simulate_akm.

• seed – As in simulate_akm.

• extra_seed (default 9) – The seed of the added columns.

• **kwargs – Passed to simulate_akm.

Returns:

polars.DataFrame

hdfe_stream.simulate.simulate_trends(n_workers=600, n_firms=60, seed=3, **kwargs)

simulate_akm plus worker-specific time trends, for FEIS models.

Each worker gets their own linear and quadratic trend in t (the year as a float), so y ~ x | worker_id[t] + firm_id is the correctly specified model and a plain worker intercept is not. It adds:

• t, t2: the year, and its centered square;

• x, v: an exogenous covariate and an unobservable;

• xe, z: an endogenous regressor and its instrument;

• wa: analytic weights;

• y: the outcome, generated with the worker trends in it.

Defaults are small because the brute-force comparison for varying slopes needs one dummy per worker per slope.

Parameters:

• n_workers – As in simulate_akm, with smaller defaults (600, 60 and 3).

• n_firms – As in simulate_akm, with smaller defaults (600, 60 and 3).

• seed – As in simulate_akm, with smaller defaults (600, 60 and 3).

• **kwargs – Passed to simulate_akm.

Returns:

polars.DataFrame

hdfe_stream.simulate.simulate_bottleneck(n_bridge=6, firms_per_block=8, movers_per_block=150, stayers_per_firm=8, block_gap=0.5, sd_worker=0.4, sd_firm=0.3, sd_noise=0.15, seed=0, keep_effects=False)

A panel whose AKM variance components are weakly identified.

Two blocks of firms with plenty of mobility inside each and only n_bridge workers moving between them: a bottleneck in the mobility network. The difference between the blocks’ average firm effects is then one linear combination of the effects that is estimated far worse than the rest, so one eigenvalue dominates the estimator’s variance – KSS’s weakly identified case, where the normal interval undercovers and the q = 1 interval is the one to use. leave_out_kss(..., se=True) reports both.

Parameters:

• n_bridge (int, default 6) – Workers moving between the blocks. Fewer is weaker; with only one or two, the leading direction rests on so few rows that even the q = 1 theory’s conditions fail, and six is a comfortable middle.

• firms_per_block (int) – The size of each block (8, 150 and 8). Movers are observed four years across two firms of their block, stayers three years at one.

• movers_per_block (int) – The size of each block (8, 150 and 8). Movers are observed four years across two firms of their block, stayers three years at one.

• stayers_per_firm (int) – The size of each block (8, 150 and 8). Movers are observed four years across two firms of their block, stayers three years at one.

• block_gap (float, default 0.5) – Added to every firm effect in the second block. It sets how far the badly estimated direction’s true value is from zero, which is what decides whether the ordinary interval merely widens or badly undercovers.

• sd_worker (float) – Standard deviations of the effects (0.4 and 0.3) and the residual (0.15).

• sd_firm (float) – Standard deviations of the effects (0.4 and 0.3) and the residual (0.15).

• sd_noise (float) – Standard deviations of the effects (0.4 and 0.3) and the residual (0.15).

• seed (int or anything numpy.random.default_rng accepts, default 0) – The seed.

• keep_effects (bool, default False) – Also return the effects that generated the data, as in simulate_akm.

Returns:

polars.DataFrame – With worker_id, firm_id and log_earn (and, with keep_effects, true_worker_effect and true_firm_effect).