Low-level interface

StreamingHDFE and StreamingGLM take column names, or (name, Polars expression) pairs, instead of a formula, and need no pyfixest. The formula functions build one of them for each set of fixed effects.

class hdfe_stream.StreamingHDFE

Low-level interface: column names or Polars expressions instead of a formula (see feols_stream for formulas). It needs no pyfixest.

Build the estimator, then call fit on the data.

Parameters:

• y (str, tuple or list) – Dependent variable(s): a column name, or a list of names and/or (name, polars.Expr) pairs. Each y is estimated as its own model (sharing one solve) unless models says otherwise.

• x (str, tuple, list or None) – Covariates, as for y. May be empty (a model of fixed effects only).

• fe (sequence of str) – Fixed-effect dimensions, e.g. ["worker_id", "firm_id", "year"] or ["worker_id", "firm_id^year"] (^ interacts columns). One dimension may carry varying slopes, "worker_id[t]" or "worker_id[t, t2]" (fixest syntax: fixed-effect intercepts plus slopes); it is then the streamed dimension. One dimension is a within regression: the streamed dimension is demeaned group by group and there is nothing left to solve for. None (or []) is ordinary least squares: the rows are read in batches and never grouped, and an "Intercept" column is added to each model’s covariates unless models lists them explicitly.

• workdir (str or Path, optional) – Directory under which each fit creates its own run directory (hdfe_run_<time>_<id>/) for intermediate and result files. Default: the HDFE_STREAM_WORKDIR environment variable if set, else the system temporary directory (honors TMPDIR). Put it on a disk with room for several times the input’s Parquet size (the sorted working copy of the rows is less compressible than typical input; about 4-5x in tests); diagnostics["disk_peak_gb"] reports it.

• stream (str, optional) – Name of the dimension to stream. Default: the dimension with slopes, else the highest approximate cardinality (ties: listed first). The choice is logged and stored in diagnostics["stream"]; pass it explicitly when cardinality is a poor guide (e.g. firm-year effects that outnumber workers in a subsample).

• weights (str or polars.Expr, optional) – Column name or Polars expression with strictly positive weights (weighted least squares; rows with missing weights are dropped).

• weights_type ({"aweights", "fweights"}, default "aweights") – Analytic weights (N = rows) or frequency weights (N = sum of weights), as in pyfixest.

• fixef_rm ({"singleton", "none"}, default "singleton") – "singleton" (as in pyfixest) drops observations whose level of some fixed effect appears in no other observation, repeating until none is left, since dropping one can create another in a different dimension. They contribute nothing to the coefficients but count toward N, the number of clusters and the fixed-effect parameters. "none" keeps them. The search takes a few extra passes over the fixed-effect columns (one per dimension per round, and typically two or three rounds), and holds one row per singleton level in memory. It runs after missing values are removed and does not look at weights; the number dropped is in diagnostics["singletons"]. Nothing is dropped without fixed effects.

• connected ({"largest", "leave_one_out"}, optional) – "largest" keeps only the rows in the largest connected set (by rows) of the pair connected_fe, by default the first two fixed effects as listed (worker_id and firm_id in "worker_id + firm_id + year"): the largest component of the graph in which each row joins its two levels, as connected_components finds it. "leave_one_out" keeps the leave-one-out connected set of the pair instead, the sample KSS estimation needs, with its first dimension as the workers (connected_components(..., leave_one_out=True) describes it); leave_out_kss uses it. Either is found among the rows left after missing values are removed (and min_degree has pruned), and comes before the singletons are dropped (and before the separation check of a GLM), so the fit is that of the data restricted beforehand. With two fixed effects the set stays as it is after them; with more, dropping a singleton of a further dimension can split it, and connectedness of one pair is necessary but not sufficient for identification (see connected_components). Default None keeps every component. HDFEResult.sample reports the rows dropped as "disconnected" or "leave_one_out"; the counts are in diagnostics["connected"].

• min_degree (int or (int, int), default 1) – Drop the rows of levels of the pair connected_fe with fewer distinct partners in the other than this, repeating until none is left, before connected and the singletons; an int, or a pair (first, second). (2, 1) keeps movers only: workers seen at two or more firms, with "worker_id + firm_id". As in connected_components, which documents it. HDFEResult.sample reports the rows as "min_degree"; the count is in diagnostics["connected"]["min_degree"]. The default, 1, drops nothing.

• connected_fe (str or (str, str), optional) – The pair of fixed effects connected and min_degree work on, as "worker_id + firm_id" or ("worker_id", "firm_id"); both must be fixed effects of the model. Default: the first two listed.

• models (list of dict, optional) – {"fml", "y", "x"} dicts selecting which of the variables enter each model; all models share passes 0-1 and the solve. An IV model adds "iv": {"endog": [...], "z": [...]}, where z lists the full instrument set (excluded instruments and exogenous regressors).

• solver ({"auto", "explicit", "stream_cg", "within"}, default "auto") – "auto": "explicit" unless S would exceed max_s_gb, then "stream_cg". "explicit": build S = D_o’ M_0 D_o once as a sparse matrix and iterate in memory; its memory is about nnz(S), which grows with distinct level co-occurrences (firm-year pairs, firm pairs sharing a mover), not with rows. "stream_cg": never form S; each iteration streams the identifying cells; the fallback when S is too big. "within": identifying cells in memory, solved by the within package (no varying slopes).

• precond ({"jacobi", "amg"}, default "jacobi") – "amg" is pyamg smoothed aggregation, for the explicit solver only.

• tol (float, default 1e-10) – Conjugate-gradient convergence tolerance, on each column’s residual relative to its right-hand side.

• maxiter (int, default 5000) – The most conjugate-gradient iterations; a solve that reaches it warns.

• max_s_gb (float, optional) – Memory budget, in GB, for the reduced system: the cap on the explicit matrix S, above which solver="auto" uses "stream_cg", and on the per-thread accumulators that "stream_cg" and the leave-out passes use, which together take at most half of it (with fewer threads on those passes when it binds). Default: 25% of physical memory.

• rhs_block (int, default 8) – Variables solved at a time; bounds solver memory at about 6 x (total non-streamed levels) x rhs_block doubles.

• batch_rows (int, default 2_000_000) – Rows per chunk when streaming Parquet, and per slice of whole groups of the streamed dimension when building the cell table.

• rows_per_bucket (int, default 20_000_000) – Target rows per bucket: rows are hash-partitioned by the streamed dimension into buckets so the sort is bucket-sized.

• n_buckets (int, optional) – The number of buckets, instead of rows_per_bucket.

• row_group_size (int, default 500_000) – Rows per row group of the working and result Parquet files.

• assembly ({"auto", "cells", "rows"}, default "auto") – "cells": within-cell cross-products accumulated in pass 1; "rows": an extra row pass after the solve, with no m^2 aggregation columns; "auto": cells when there are at most 8 variables.

• collapse_cells ({"auto", True, False}, default "auto") – Whether the rows are collapsed into cells, one per combination of fixed-effect levels within each group of the streamed dimension. Collapsing makes the later passes cheaper when combinations repeat (a worker’s years at one firm, with worker_id + firm_id). It costs a wider sort key in pass 0, and buys nothing when they do not repeat (worker_id + firm_id^year on a worker-year panel, or worker_id + firm_id + year), where every row is its own cell. "auto" counts the distinct combinations approximately in pass 0’s scan, at little cost, and collapses unless they are at least 90% of the rows; True always collapses, False never does (rows that happen to be adjacent and identical still share a cell). The estimates are the same either way. GLMs and fits with keep_intermediates always collapse. diagnostics["cells_collapsed"] says which happened.

• collin_tol (float, default 1e-10) – Tolerance for dropping collinear covariates (pyfixest’s algorithm applied to the residualized X’X).

• collin_tol_rel (float, default 1e-6) – Also drop a covariate when its variation left after removing the fixed effects is below this share of its raw variation (the default means residual SD < 0.1% of raw SD).

• cells_in_memory (bool, default False) – Load the identifying cells into memory instead of memory-mapping them.

• n_threads (int, optional) – Numba threads. Default: all available.

• scratch_mb (int, default 32) – Budget, in MB, for the row-blocked scratch the leave-out passes hold. Those passes want about a dozen float64 columns per row of the chunk, so the chunk is shrunk below batch_rows as needed to stay inside this. It bounds memory only: every random vector involved is keyed on the absolute row ordinal, so results do not depend on it. Raising it trades memory for speed and the default sits at the knee: on a 1.3M-row panel, 32 MB runs the correction in 30s holding 67 MB, while 512 MB gets to 20s but holds 669 MB.

• keep (sequence of str, optional) – Extra columns to carry into the residual output.

• row_id (str, default "row_id") – Name of the column that gives each row’s position in the source (counted before any row is dropped). It is written to the residual file, so HDFEResult.resid joins back onto the source on it, and HDFEResult.sample adds it to the source. The source must not already have a column of this name.

• outputs ({"auto", "keep"}, default "auto") – "auto": result files (residuals, fixed effects) are deleted when all result objects of the fit are garbage-collected or the interpreter exits; keep the result alive while using resid() and fixef(), or collect or sink what you need. "keep": result files stay until HDFEResult.cleanup or hdfe_stream.cleanup. Intermediates are always removed.

• save_resid (bool, default True) – Write the row-level residual file (the largest output). With False, resid() is unavailable.

• keep_intermediates (bool, default False) – Keep the intermediate files, and on errors the run directory. Needed for leave-out (Kline-Saggio-Solvsten) estimation when leaving out an observation, which reads the sorted rows and cell tables after the fit; see HDFEResult.leave_out_kss. Otherwise for debugging. It also widens the pass-0 sort key so the stored row order is a function of the data, which is what makes the leave-out draws reproducible.

• verbose (bool, default True) – Emit progress messages.

• logger (logging.Logger, optional) – If given, progress messages go to it (at log_level) as they happen, per-iteration CG residuals go to it at DEBUG (and are not printed without a logger), and warnings (dropped collinear variables, non-convergence, negative variances) go to logger.warning. Default: print to stdout and warn through the warnings module.

• log_level (int, default logging.INFO) – The level progress messages are logged at.

Examples

>>> import polars as pl
>>> est = StreamingHDFE(
...     y="y",
...     x=["x", ("late", (pl.col("year") >= 2010).cast(pl.Float64))],
...     fe=["worker_id[t]", "firm_id"],
...     workdir="scratch")
>>> fit = est.fit("data/*.parquet", vcov={"CRV1": "firm_id"})

__init__(y, x, fe, workdir=None, solver='auto', precond='jacobi', keep=(), tol=1e-10, maxiter=5000, batch_rows=2000000, row_group_size=500000, n_buckets=None, rows_per_bucket=20000000, cells_in_memory=False, rhs_block=8, assembly='auto', max_s_gb=None, collin_tol=1e-10, collin_tol_rel=1e-06, n_threads=None, scratch_mb=32, fixef_rm='singleton', connected=None, min_degree=1, connected_fe=None, row_id='row_id', weights=None, weights_type='aweights', stream=None, models=None, verbose=True, logger=None, log_level=20, outputs='auto', save_resid=True, keep_intermediates=False, collapse_cells='auto')

classmethod __new__(*args, **kwargs)

fit(source, vcov='iid', cluster=(), fe_dof='exact', bootcluster=None)

Fit the model(s) to the data.

Parameters:

• source (str, Path, polars.LazyFrame or polars.DataFrame) – A Parquet path or glob, a Polars LazyFrame, or a DataFrame.

• vcov (str or dict, default "iid") – The default vcov: "iid", "hetero" (or "HC1"), {"CRV1": var} or "CRV1:var"; multi-way clustering as {"CRV1": "a+b"} (any number of ways). {"CRV3": var} is the cluster jackknife (one-way, OLS), available when every fixed effect is nested within the clusters or there are none. iid and hetero are always computed; switch with HDFEResult.with_vcov.

• cluster (str or sequence of str, optional) – Further cluster specifications to compute CRV1 for, e.g. ["firm_id", "worker_id+firm_id"]. Any column or "a^b" combination works; the streamed dimension and the other fixed effects reuse their codes.

• fe_dof ({"exact", "pyfixest"}, default "exact") – "exact" counts one redundant level per connected component of the pair of fixed-effect dimensions with the most components, plus one per further dimension: exact with two dimensions, and with more a conservative count (never more redundant levels than there are) that does not depend on their order. For CRV1, a pair of fixed effects both nested in the clusters is not used. "pyfixest" counts one per dimension beyond the first.

• bootcluster (str, optional) – A clustering (one column, or "a^b") to keep the per-cluster sums for that HDFEResult.wildboottest, hdfe_stream.rwolf and hdfe_stream.wyoung need: about G x k^2 numbers per model, for G clusters and k coefficients, collected in the residual pass that the fit makes anyway. CRV1 by it is computed too. OLS only (not IV), and every fixed effect must be nested within the clusters, as for CRV3: firm effects clustered by state are, if no firm is in two states; year effects never are, so put them in as covariates instead ("y ~ treat + i(year) | county" estimates the same treat as "y ~ treat | county + year").

Returns:

HDFEResult or HDFEMulti – One result, or several when there are several models.

class hdfe_stream.StreamingGLM

Low-level GLM interface: column names or Polars expressions instead of a formula (see fepois_stream and feglm_stream for formulas). One dependent variable, any number of models of it.

Takes StreamingHDFE’s parameters, and these. assembly does not apply (the cross-products come from each step’s row pass), and varying slopes, IV and CRV3 are not available. With solver="auto" the first step builds S and later steps use "stream_cg".

Parameters:

• y (str or tuple) – The dependent variable: a column name, or a (name, polars.Expr) pair.

• x (str, tuple, list or None) – Covariates, as in StreamingHDFE.

• fe (sequence of str) – Fixed-effect dimensions, as in StreamingHDFE, without varying slopes.

• family ({"poisson", "logit", "probit"}, default "poisson") – Poisson (log link), logit or probit.

• offset (str or polars.Expr, optional) – Column name or Polars expression added to the linear predictor with a coefficient of one (e.g. log exposure).

• iwls_tol (float, default 1e-8) – IRLS convergence tolerance on the relative change in deviance (pyfixest’s default).

• iwls_maxiter (int, default 25) – The most IRLS steps (pyfixest’s default).

• separation_check (bool, default True) – Drop fixed-effect levels whose outcome is always zero (Poisson) or never varies (logit, probit), for logit and probit repeating until none is left. Covariates that separate the outcome are not detected.

• **options – StreamingHDFE’s other parameters.

Examples

>>> import polars as pl
>>> fit = StreamingGLM("visits", ["x"], ["worker_id", "firm_id", "year"],
...                    family="poisson", offset=pl.col("log_exposure"),
...                    workdir="scratch").fit("data/*.parquet")

__init__(y, x, fe, family='poisson', offset=None, iwls_tol=1e-08, iwls_maxiter=25, separation_check=True, **options)

classmethod __new__(*args, **kwargs)

StreamingGLM has the same fit method as StreamingHDFE.