FastHDFE

This repository provides FastHDFE, a Python package that can be used to estimate linear and multiplicative models with high-dimensional fixed effects.

The core of the package lies in the two functions called reghdfe and ppmlhdfe. The reghdfe-command contains an estimator for linear models, which is optimized for speed and will produce the common output typically expected after a linear regression. Similarly, the ppmlhdfe-command contains an estimator for multiplicative models with high-dimensional fixed effects. As we envision the commands to be useful for different fields, the commands allow the choice of the structure of fixed effects as well as the type of standard errors.

This page covers:

• Requirements
• Test dataset for code examples
• The reghdfe command
• The ppmlhdfe command
• Citation details

Requirements

• Numpy
• Scipy
• Pandas
• tabulate

Test Data

The project uses a small dataset useful for testing. This dataset contains data for 62 exporters, 62 importers and 27 years. Use it like this to execute the example calls:

import pandas as pd
df = pd.read_csv('https://www.usitc.gov/data/gravity/example_trade_and_grav_data_small.csv')
df['lntrade_value'] = np.log(cls._df['trade_value'])

reghdfe

reghdfe is a command available via our package to estimate linear models with high-dimensional fixed effects that has the following structure:

1. Inputs:

   1. dependent variable,

   2. a set of regressors,

   3. a flexible way to specify the type of fixed effects one wants to include,

   4. options to choose convergence criteria,

   5. options to choose from various types of standard errors (homoskedastic standard errors, heteroskedasticity-robust standard errors, clustered standard errors),

   6. option to enable/disable conjugate gradient acceleration,

   7. option to save fixed effects.

2. Outputs:

   1. coefficient estimates,

   2. standard errors,

   3. variance-covariance matrix,

   4. t-values,

   5. number of observations,

   6. model diagnostics (e.g., Aikake Information Criterion, Bayesian Information Criterion),

   7. fitted values, which can be merged with the input data,

   8. residuals, which can be merged with the input data,

   9. optional: fixed effects, which can be merged with the input data.

3. Syntax:

   reghdfe(y,
           x, 
           fixedeffects              = None,
           data                      = None,
           setype                    = "none",
           clustvars                 = None,
           FEgen                     = False,
           verb                      = 0,
           eps_MAP                   = 1e-9,
           minsubiter                = 0,
           maxsubiter                = 1000,
           accel                     = True,
           nt                        = 0,
           print_stata_style_summary = True)
   

   1. Parameters:

      • y: single entry list containing the name of the dependent variable as a string.

      • x: string list of the names of the explanatory variables.

      • fixedeffects: string list of the names of the fixed effects. Variables can be combined by "#", e.g. "state#year" creates state-year fixed effects. Default is None, which causes estimation with a constant.

      • data: data frame containing all relevant variables.

      • setype: string, type of standard errors. Allowed values are "none", "homoskedastic", "heteroskedastic" ("r", "robust" and "HC1" are synonyms), "cluster" ("multiway" as synonym). Default is None.

      • clustvars: string list of cluster variables. Max length 4, ignored if setype is not "cluster". Default is None.

      • FEgen: bool, if True, the observation-wise sum of the fixed effects, as well as the different types of fixed effects, are returned.

   2. Additional Parameters:

      • accel: bool, if True, MAP employs conjugate gradient acceleration technique following Hernández-Ramos et al. ( 2011 ). Default is True.

      • nt: integer, sets the number of threads to use. Default is 0, which uses "system settings". Has no effect if the backend is Python.

      • verb: int $\in {0, 1, 2}$, level of "in-between" output displayed (0 is none). Default is 0.

      • eps_MAP: float $\in (0, \infty)$, convergence criterion of the alternating projections algorithm given. Default is $\mathtt{10^{-9}}$.

      • minsubiter: minimum number of iterations within the demeaning process. Default is 0.

      • maxsubiter: maximum number of iterations within the demeaning process. Default is 1000.

      • print_stata_style_summary: bool, regression output table is printed. Default is True.

   3. Returns a class object with the following attributes:

      • params: point estimates.

      • bse: standard errors if calculated, "NULL" otherwise.

      • covmat: variance-covariance matrix if calculated, "NULL" otherwise.

      • covmat_nofix: variance-covariance matrix prior to the "eigenvalue fix". Only present if eigenvalue_fix is True.

      • eigenvalue_fix: True only if the covariance matrix was not positive semi-definite and the eigenvalue fix had to be applied (negative eigenvalues replaced by zero); otherwise the attribute is not set.

      • cov_type: string, type of standard errors.

      • clustvars: list, string list with cluster dimensions if cov_type is "cluster", "NULL" otherwise.

      • tvalues: $t$-values for $H_0 = 0$.

      • pvalues: $p$-values for a two-sided $t$-test with $H_0 = 0$.

      • nobs: int, number of observations used for estimation. Some observations may have been excluded and do not count in here.

      • df_model: number of parameters in $X$ plus the theoretical upper bound of fixed effects (collinear fixed effects count as estimated parameters here; this may lead to a slight over-estimation of small-sample corrected standard errors).

      • df_resid: nobs $-$ df_model.

      • aic: Aikake Information Criterion.

      • bic: Bayesian Information Criterion.

      • llf: value of the log-likelihood.

      • deviance: sum of the squared deviance residuals, equal to the sum of squared residuals.

      • pearson_chi2: Pearson's $\chi^2$-value, equal to the sum of squared residuals.

      • family_name: "Gaussian distribution".

      • family_link: "Identity function", ( $g(\cdot)$ ).

      • scale: ML estimate of $\sigma^2$. To obtain the OLS estimates, multiply by nobs/(nobs-1).

      • method: "OLS".

      • fit_history: dict with key 'iteration', equal to 1, as analytical OLS solution is computed.

      • fittedvalues: array of the fitted values.

      • fe_rowsum: array of the observation-wise sum of the fixed effects $\boldsymbol{\eta}-\mathbf{x}'\boldsymbol{\beta}$.

      • fixed_effects: dictionary of all estimates of all fixed effects specified.

      • resid: array of residuals.

      • msg: list, information/warning messages during the estimation process.

      • formula_R: string, formula as you would type it in R.

      • formula_Stata: string, formula as you would type it in Stata.

      • model: GLM class, for compatibility with the gme-package.

   4. Example:

      from fasthdfe import reghdfe
      result = reghdfe(y            = ['lntrade_value'],
                       x            = ['LN_DIST','agree_pta','contiguity','common_language'],
                       fixedeffects = ['exp_ind#industry_id#year','imp_ind#industry_id#year','exp_ind#imp_ind#industry_id'],
                       data         = df,
                       setype       = 'cluster',
                       clustvars    = ['exp_ind','imp_ind','year','industry_id'])
      

ppmlhdfe

ppmlhdfe is a user-friendly Python command to estimate multiplicative models with high-dimensional fixed effects with PPML that has the following structure:

1. Inputs:

   1. dependent variable,

   2. a set of regressors,

   3. a flexible way to specify the type of fixed effects one wants to include,

   4. options to choose convergence criteria,

   5. options to choose from various types of standard errors (homoskedastic standard errors, heteroskedasticity-robust standard errors, clustered standard errors),

   6. options for initial values of the linear predictor (simple, OLS-based, or user-supplied starting values for selected coefficients),

   7. options to detect and drop separated observations (fixed-effect-style separation, iterative rectifier, or none, mirroring Stata's separation() option),

   8. option to keep or drop singleton observations,

   9. an offset option for constrained coefficients,

   10. option to enable/disable conjugate gradient acceleration,

   11. option to standardize regressors before estimation for improved numerical stability (mirroring Stata's standardize_data(1) default),

   12. option to save fixed effects.

2. Outputs:

   1. coefficient estimates,

   2. standard errors,

   3. variance-covariance matrix,

   4. t-values,

   5. number of observations,

   6. diagnostics on the number of observations kept and dropped (total, invalid values, singletons/fe-separation, and separation detected by the iterative rectifier),

   7. model diagnostics (e.g., Aikake Information Criterion, Bayesian Information Criterion, deviance and pseudo-likelihood),

   8. fitted values, which can be merged with the input data,

   9. residuals, which can be merged with the input data,

   10. optional: fixed effects, which can be merged with the input data.

3. Syntax:

   ppmlhdfe(y,
            x,
            fixedeffects              = None, 
            data                      = None, 
            setype                    = "none", 
            clustvars                 = None,
            off                       = None, 
            eta_0                     = None, 
            FEgen                     = False,
            start                     = "simple", 
            olsguess                  = False,   # olsguess kept for back-compat
            separation                = "ir",
            sep_tol                   = 1e-5, 
            sep_maxiter               = 100, 
            sep_K_penalty             = 1e6,
            keepsingletons            = False,
            nocheck                   = False, 
            verb                      = 0,
            eps_beta                  = 1e-8, 
            miniter                   = 1, 
            maxiter                   = 1000,
            eps_MAP                   = 1e-9, 
            minsubiter                = 0, 
            maxsubiter                = 16000,
            start_inner_tol           = 1e-4,
            standardize_data          = True,
            accel                     = True, 
            nt                        = 0, 
            print_stata_style_summary = True)
   

   1. Parameters:

      • y: single entry list containing the name of the dependent variable as string.

      • x: string list of the names of the explanatory variables.

      • fixedeffects: string list of the names of the fixed effects. Variables can be combined by "#", e.g. "state#year" creates state-year fixed effects. Default is None, which causes estimation with a constant.

      • data: data frame containing all relevant variables.

      • setype: string, type of standard errors. Allowed values are "none", "homoskedastic", "heteroskedastic" ("r", "robust" and "HC1" are synonyms), "cluster" ("multiway" as synonym). Default is None.

      • clustvars: string list of cluster variables. Max length 4, ignored if setype is not "cluster". Default is None.

      • FEgen: bool, if True, the observation-wise sum of the fixed effects (fe_rowsum), as well as the different types of fixed effects (fixed_effects), can be accessed after estimation (when the object is saved with "results", then with result.fe_rowsum and result.fixed_effects, respectively).

      • eta_0: single entry list containing the string of the initial values for the linear predictor $\eta$. If set, overrides start. Default is None.

      • start: strategy for the initial values of $\eta$. Accepts one of the following:

        • “simple” (default): $\eta = \ln\big((y + \overline{y})/2\big)$, the Stata ppmlhdfe-style default.

        • “ols”: use OLS on the fixed-effect-demeaned data, then take $\eta = \ln(\text{fitted values})$.

        • dict: user-supplied starting values for selected regression coefficients, e.g. {“log_distance”:-1.0, “contiguity”: 0.4}. Unspecified variables start at zero.

      • olsguess: bool, deprecated. Kept for backward compatibility. If True, equivalent to start=“ols”. Default is False.

      • off: single entry list containing the string of the offset variable.

      • separation: string, list of strings, bool, or None, default “ir”. Controls detection and removal of separated observations, mirroring the separation() option of Stata's ppmlhdfe. Allowed values:

        • “fe”: drop observations whose fixed-effect group has $y$ constant across all observations (this is already performed unconditionally by the singleton/constant-$y$ loop; listing it is informational).

        • “ir” / “relu”: iterative rectifier; detects separation arising from both fixed effects and regressors (Correia et al. 2021).

        • “none” / “off” / False / None: disable additional separation checks (the implicit "fe"-style check still runs).

        • “default” / “auto” / True: equivalent to [“ir”] (the default, matching Stata).

        • list: apply all listed methods, e.g. [“fe”,“ir”].

      • sep_tol: float, threshold above which an observation is flagged as a separation certificate by the iterative rectifier. Default is $10^{-5}$.

      • sep_maxiter: int, maximum number of outer iterations of the iterative rectifier. Default is $100$.

      • sep_K_penalty: float, penalty weight on $y>0$ observations in the IR weighted regression. Too large can cause numerical issues in the inner weighted MAP; too small allows leakage onto $y>0$ observations. Default is $10^{6}$.

      • keepsingletons: bool, default False (matching Stata's default). If False, drop any fixed-effect group with no within-group variation in $y$ (singletons, constant-$y$ groups, and classical fe-style separation). If True, only drop groups where $y$ is identically zero---the minimum required for Poisson identification---and keep singletons and constant-positive-$y$ groups. Useful for replicating analyses that report the full sample, but singletons carry no identifying information and can distort cluster-robust standard errors.

      • standardize_data: bool, default True (matching Stata's standardize_data(1)). If True, each regressor is divided by its sample standard deviation before estimation, and the resulting coefficients and variance-covariance matrix are rescaled to original units at the end. This substantially improves numerical stability when regressors have very different magnitudes (a common case in trade gravity, where dummies sit next to log distances or large monetary values), and prevents the cluster-robust variance computation from becoming near-singular when a regressor's identifying variation is concentrated in a small subset of observations after fixed-effect absorption.

   2. Additional Parameters:

      • accel: bool, if True, MAP employs conjugate gradient acceleration technique following Hernández-Ramos et al. ( 2011 ). Default is True.

      • nt: integer, sets the number of threads to use. Default is 0, which uses "system settings".

      • nocheck: bool, if True, no consistency check is skipped. Default is False.

      • verb: int $\in {0, 1, 2}$, level of output displayed (0 is none). Default is $0$.

      • eps_beta: float $\in (0, \infty)$, convergence criterion of the IRLS loop, $||\boldsymbol{\beta}^{(r+1)}-\boldsymbol{\beta}^{(r)}||_2 < \mathtt{eps_beta}$. Convergence is declared only once this criterion is satisfied in two consecutive iterations, mirroring Stata's convention. Default is $\mathtt{10^{-8}}$, matching Stata ppmlhdfe's tolerance(1e-8).

      • eps_MAP: float $\in (0, \infty)$, target convergence criterion of the alternating projections algorithm. The effective inner tolerance applied at convergence is $\max(\min(\mathtt{eps_MAP},,0.1\cdot\mathtt{eps_beta}),,10^{-12})$, following Stata reghdfe's rule that the inner solve must be at least ten times tighter than the outer IRLS criterion but never tighter than machine precision. Default is $\mathtt{10^{-9}}$, matching Stata reghdfe's itol(1e-9).

      • start_inner_tol: float $\in (0, \infty)$, starting tolerance for the inner partialling-out loop on the first IRLS iteration. The inner tolerance is loosened to this value early in IRLS (when the outer $\boldsymbol{\beta}$-change is large) and tightens proportionally to the outer convergence measure across iterations, with the floor described under eps_MAP. This adaptive schedule avoids spending inner iterations on precision that the outer loop cannot yet exploit, while guaranteeing the required tightness at convergence. Default is $\mathtt{10^{-4}}$, matching Stata reghdfe's start_inner_tol(1e-4).

      • miniter: minimum number of iterations of the IRLS loop. Default is $1$.

      • maxiter: maximum number of iterations of the IRLS loop. Default is $1000$, matching Stata ppmlhdfe's maxiter(1000).

      • minsubiter: minimum number of iterations within the demeaning process. Default is $0$.

      • maxsubiter: maximum number of iterations within the demeaning process. Default is $16000$, matching Stata reghdfe's iterations(16000).

      • print_stata_style_summary: bool, regression output table and Stata-style drop/separation notices are printed. Default is True.

   3. Returns a class object with the following attributes:

      • params: point estimates.

      • bse: standard errors if calculated, "NULL" otherwise.

      • covmat: variance-covariance matrix if calculated, "NULL" otherwise.

      • covmat_nofix: variance-covariance matrix prior to the "eigenvalue fix". Only present if eigenvalue_fix is True.

      • eigenvalue_fix: True only if the covariance matrix was not positive semi-definite and the eigenvalue fix had to be applied (negative eigenvalues replaced by zero); otherwise the attribute is not set.

      • cov_type: string, type of standard errors.

      • clustvars: list, string list with cluster dimensions if cov_type is "cluster", "NULL" otherwise.

      • tvalues: $t$-values for $H_0 = 0$.

      • pvalues: $p$-values for a two-sided $t$-test with $H_0 = 0$.

      • nobs: int, number of observations used for estimation. Some observations may have been excluded and do not count in here.

      • df_model: number of parameters in $X$ plus the theoretical upper bound of fixed effects (collinear fixed effects count as estimated parameters here; this may lead to a slight over-estimation of small-sample corrected standard errors).

      • df_resid: nobs $-$ df_model.

      • aic: Aikake Information Criterion.

      • bic: Bayesian Information Criterion.

      • llf: value of the log-likelihood.

      • deviance: sum of the squared deviance residuals. For $y=0$, the deviance is approximated using L'Hôpital's rule.

      • pearson_chi2: Pearson's $\chi^2$-value.

      • family_name: "Poisson distribution".

      • family_link: "Natural logarithm function", ( $g(\cdot)$ ).

      • scale: 1, per definition of the Poisson distribution.

      • method: "IRLS", Iteratively re-weighted least squares.

      • fit_history: dict with key 'iteration', which is the number of iterations until convergence.

      • fittedvalues: array of the fitted values $\hat{y}=\mu = g^{-1}(\eta) = \exp(\eta)$.

      • fe_rowsum: array of the observation-wise sum of the fixed effects $\boldsymbol{\eta} - \mathbf{x}'\boldsymbol{\beta}$.

      • fixed_effects: dictionary of all estimates of all fixed effects specified.

      • resid_additive: array of additive residuals ( $y-\hat{y}$ ).

      • resid_multiplicative: array of multiplicative residuals ( $y/\hat{y}$ ).

      • msg: list, information/warning messages during the estimation process.

      • formula_R: string, formula as you would type it in R.

      • formula_Stata: string, formula as you would type it in Stata.

      • model: GLM class, for compatibility with the gme-package.

      • separation_methods: list of strings, the separation methods that were applied (e.g. [“ir”]).

      • num_separated: int, number of observations dropped by the iterative rectifier. Equal to N_dropped_ir_sep.

      • keepsingletons: bool, the value of the keepsingletons argument passed in.

      • N_full_initial: int, number of observations in the input data frame before any drops.

      • N_dropped_invalid: int, number of observations dropped due to missing, infinite, or negative values.

      • N_dropped_fe_sep: int, number of observations dropped by the singleton / constant-$y$ / fe-style separation loop.

      • N_dropped_ir_sep: int, number of observations dropped by the iterative rectifier.

   4. Example:

      from fasthdfe import ppmlhdfe
      result = ppmlhdfe(y            = ['trade'],
                        x            = ['LN_DIST','agree_pta','contiguity','common_language'],
                        fixedeffects = ['exp_ind#industry_id#year','imp_ind#industry_id#year','exp_ind#imp_ind#industry_id'],
                        data         = df,
                        eps_beta     = 1e-5,
                        setype       = 'cluster',
                        clustvars    = ['exp_ind', 'imp_ind','year','industry_id'])
      

Citation

@misc{ubt_epub9520,
           title = {Fast Estimation of Linear and Poisson Models with High-Dimensional Fixed Effects in Python : The FastHDFE package},
            note = {This is a technical report corresponding to the Python package https://pypi.org/project/fasthdfe},
         address = {Bayreuth, Germany},
           month = {July},
            year = {2026},
        keywords = {high-dimensional fixed effects, Poisson pseudo-maximum-likelihood, gravity models, method of alternating projections, multi-way clustering, Python, econometric software},
          author = {Larch, Mario and Schoenfeld, Mirco and Shikher, Serge},
             url = {https://epub.uni-bayreuth.de/id/eprint/9520/},
             doi = {https://doi.org/10.15495/EPub_UBT_00009520}
}