*! version 0.6.2 16feb2024

Title

reghdfejl -- Linear regression with multiple fixed effects, accelerated by Julia

Syntax

Fixed effects regression:

reghdfejl
depvar [indepvars]
[if] [in] [weight] , [absorb(absvars, [savefe])] [reghdfejl options]

Instrumented fixed effects regression:

reghdfejl
depvar [indepvars] (endogvars = instruments)
[if] [in] [weight] , [absorb(absvars, [savefe])] [reghdfejl options]

reghdfejl mask

reghdfejl unmask

Partial fixed effects out of variables:

partialhdfejl
depvars
[if] [in] [weight] , absorb(reghdfejl##absorb:absvars) [partialhdfejl options]

reghdfejl option

absorb(absvars, [savefe]) categorical variables representing the fixed effects to be absorbed

vce(reghdfejl##model_opts:vcetype) vcetype may be unadjusted (default), robust, bs/bootstrap, cluster fvvarlist (allowing multi-way clustering)

residuals(newvar) save regression residuals; required for postestimation "predict <varname>, d"

residuals(newvar) save regression residuals; required for postestimation "predict <varname>, d"

ivreg2 call ivreg2 for IV estimation

iterate(#) maximum number of iterations; default is 16,000

keepsingletons do not drop fixed effect singletons

nosample do not create e(sample), saving some space and speed

compact temporarily saves all data to disk in order to free memory

threads(#) number of CPU threads Julia should use

gpu use NVIDIA or Apple Silicon GPU

level(#) set confidence level; default is normally 95

display options standard options governing results display

verbose show and leave behind Julia regression command & data

depvar and indepvars may contain factor variables and time-series operators. depvar
cannot be of the form i.y though, only #.y (where # is a number)

partialhdfejl option

absorb(absvars) categorical variables representing the fixed effects to be absorbed

tolerance(#) criterion for convergence. default is 1e-8

iterate(#) maximum number of iterations; default is 16,000

compact temporarily saves all data to disk in order to free memory

gpu use NVIDIA or Apple Silicon GPU

generate(varlist) names for partialled-out results

prefix(string) prefix stub for partialled-out results

replace overwrite any existing variables identified by generate() or prefix()

Exactly one of generate() and prefix() must be specified.

vce(bootstrap, bsoptions) accepts these standard bootstrap options: robust, cluster(varname), seed(#), reps(#), mse,
size(#), saving(filename [, double replace]). In addition it accepts a procs() option to accelerate the bootstrap through multitasking.

Description

reghdfejl is designed as a slot-in replacement for reghdfe. It is missing some features of reghdfe. And it is
not guarantee to exactly match reghdfe's results. But it
can run ~10 times faster because it relies on the Julia program FixedEffectsModel.jl,
by Matthieu Gomez, which implements similar methods (Correia 2016).
reghdfejl also fits instrumental variables models
with two-stage least squares. In this capacity it is not by default as full-featured as ivreghdfe, because it does not work as a wrapper for
ivreg2. The ivreg2 option of reghdfejl overrides that default, but then reghdfejl is no faster than ivreghdfe.

To run, reghdfejl requires that the Stata command jl be installed; "ssc install julia" should suffice. It also needs
Julia 1.9.4 or later, which is free. See these installation instructions.

Because Julia performs just-in-time compilation, and because reghdfejl may need to install Julia packages, there can be long lags on first
use. It can also take longer the first time in a session that you use a feature, such as multiway clustering, that forces compilation
of another function in the Julia package.

If reghdfejl appears to be failing to install the needed packages,
you can try intervening manually: start Julia outside of Stata, hit the "]" key to enter the package manager, and type
add <pkgname> for each package. The needed packages are Vcov, FixedEffectModels, DataFrames, and Metal (for
Macs) or CUDA (otherwise).

reghdfejl lacks some reghdfe features that are typically secondary for users:

* It does not correct the estimates of the degrees of freedom consumed by absorbed fixed effects for collinearity
and redundance among fixed-effect dummies. reghdfe
displays these corrections in a table after the main results. reghdfejl does not.

* It does not offer the Group FE features.

* It does not accept frequency weights ("fweights").

* It does not allow control over whether the constant term is reported. The constant is always absorbed.

* It does not offer options such as technique() that give finer control over the algorithm. But these are largely obviated
by reghdfejl's speed.

reghdfejl starts by copying the data needed for estimation into a Julia DataFrame. Duplicating the data takes a bit of time and potentially
a lot of RAM. In
extreme cases, it will more than double the storage demand because even variables stored in Stata in small types, such as byte, will be stored
double-precision--8 bytes per value--in Julia. If the memory demand is too great, performance will plummet. reghdfejl therefore is most
useful when you have plenty of RAM, when the number of non-absorbed regressors is low, and when the number of absorbed terms is high
(for then the computational efficiency of Julia shines).

reghdfejl offers several novel features that can increase speed. The first is access to multithreading in Julia, even in
flavors of Stata that do not offer multiprocessing. The threads(#)
option pertains to this feature. But it can only reduce the number of CPU threads Julia uses. The default number--and the
maximum--is set when the Julia session inside Stata is started. It is
possible for the default to be too high as well as too low. If you set it high, then you can experiment using threads(#). See
help jl for more on determining and controlling the number of threads.

In a similar vein, reghdfejl offers accelerated bootstrapping for computing standard errors, via the bs/bootstrap
suboption of the vce() option. The results should be the same, asymptotically, as if one prefixes a reghdfejl command line with
Stata's bootstrap command. But they should come much faster because copying of data between Stata and Julia is minimized, and multiple
copies of Julia are launched for parallelization.

A final new feature is access to GPU-based computation. The gpu specifies the use of NVIDIA or Apple Silicon
GPUs for computation. Typically this modestly increases speed. On non-Apple computers, reghdfejl installs the latest
version of the Julia package CUDA.jl, which currently requires CUDA drivers 11.0 or later. (CUDA is NVIDIA's programming
interface for GPUs.) You can visit the
CUDA download site for the latest drivers.

The command reghdfejl mask redirects all reghdfe calls to reghdfejl. reghdfejl unmask stops the redirection. This is useful is when using other Stata packages that call reghdfe, such as
eventdd. It can speed up those commands as well. Since reghdfe and reghdfejl do not accept exactly the same options,
nor return exactly the same result sets, no guarantee can be given that this hack will work in any particular case.

partialhdfejl partials a given set of fixed effects out of several variables at once, restricting to the sample defined by the missingness pattern of all the
variables together. The results are stored in variables whose names are listed in a generate() option, by implied by a prefix() option. For example,
partialhdfejl x y, a(z) prefix(_) stores the partialled-out versions of x and y in _x and _y. partialhdfe throws an error
if the destination variables already exist, unless replace is also specified.

absorb() syntax

absvar

varname categorical variable to be absorbed

i.varname categorical variable to be absorbed (same as above; the i. prefix is always implicit)

i.var1#i.var2 absorb the interactions of multiple categorical variables

i.var1#c.var2 absorb heterogeneous slopes, where var2 has a different slope estimate depending on var1. Use carefully (see below)!

var1##c.var2 absorb heterogenous intercepts and slopes. Equivalent to "i.var1 i.var1#c.var2", but much faster

var1##c.(var2 var3) multiple heterogeneous slopes are allowed together. Alternative syntax: var1##(c.var2 c.var3)

v1#v2#v3##c.(v4 v5) factor operators can be combined

- To save the estimates of specific absvars, write newvar=absvar. This works with reghdfejl, not partialhdfejl.

- However, be aware that estimates for the fixed effects are generally inconsistent and not econometrically identified.

- Using categorical interactions (e.g. x#z) is easier and faster than running egen group(...) beforehand.

- Singleton observations are dropped iteratively until no more singletons are found (see the linked article for details).

- Slope-only absvars ("state#c.time") have poor numerical stability and slow convergence. If you need those, either i) increase tolerance or
ii) use slope-and-intercept absvars ("state##c.time"), even if the intercept is redundant.
For instance if absvar is "i.zipcode i.state##c.time" then i.state is redundant given i.zipcode, but
convergence will still be much faster.

Options

absorb(absvars, [savefe]) list of categorical variables (or interactions) representing the fixed effects to be absorbed.
This is equivalent to including an indicator/dummy variable for each category of each absvar. absorb() is required.

To save a fixed effect, prefix the absvar with "newvar=".
For instance, the option absorb(firm_id worker_id year_coefs=year_id) will include firm, worker, and year fixed effects,
but will only save the estimates for the year fixed effects (in the new variable year_coefs).

If you want to run predict afterward but don't particularly care about the names of each fixed effect, use the savefe suboption.
This will delete all preexisting variables matching __hdfe*__ and create new ones as required.
Example: reghdfejl price weight, absorb(turn trunk, savefe).

vce(reghdfejl##model_opts:vcetype) specifies the type of standard error reported.

unadjusted|ols estimates conventional standard errors, valid under the assumptions of homoscedasticity and no correlation between observations even in small samples.

robust estimates heteroscedasticity-consistent standard errors (Huber/White/sandwich estimators), which still assume independence between observations.

cluster clustervars estimates consistent standard errors even when the observations
are correlated within groups. Multi-way-clustering is allowed.

vce(bootstrap, bsoptions)} and vce(bs, bsoptions) are synonyms. They request estimation of standard errors using the non-parametric or "pairs"
bootstrap. In expectation, reghdfejl ..., ... vce(bs, bsoptions) should return the same results as bs, bsoptions: reghdfejl ..., .... But the
first is faster because it avoids copying data from Stata to Julia on every replication and can exploit mulitasking. bsoptions
may include any of the following suboptions: robust, cluster(varname), seed(#), reps(#), mse, size(#),
saving(filename [, double replace]), and
procs(#). All but the last are standard bootstrap options. The last instructs reghdfejl to launch several
copies of Julia in order to run the bootstrap in parallel. The procs(#) suboption is semantically distinct from reghdfejl's threads()
option. The latter triggers low-level multitasking: the Julia package FixedEffectModels.jl spreads certain loops across multiple threads within one Julia
instance. The former instead runs
multiple copies of Julia, each of which loads and runs FixedEffectModels.jl, on just one thread each. While
the options are implemented in different ways, the principles governing the optimal number of threads and processes are the
same. See help jl.

Note that while setting the seed(#) suboption allows for exact reproducibility of results, even with the same seed, changing the procs(#)
suboption will (slightly) change results. The latter affects how the bootstrap is distributed across the Julia processes, each with its own
psuedorandom number stream.

residuals[(newvar)]} saves the regression residuals in a new variable. residuals without parenthesis saves them
in the variable _reghdfejl_resid, overwriting it if it already exists.

This option carries a small time cost but is required for subsequent calls to predict, d.

tolerance(#) specifies the tolerance criterion for convergence. The default is 1e-8.
In general, low tolerances (1e-8 to 1e-14) return more accurate results, more slowly.

ivreg2 causes reghdfejl to work more like ivreghdfe when performing instrumental variables estimation. It calls Julia to partial
fixed effects out of all other variables, and then passes the results to ivreg2. This gives access to a much wider
array of estimation options and weak identification diagnostics. However, because ivreg2 dominates the run time, this is no faster than using
ivreghdfe.

iterations(#)
specifies the maximum number of iterations; the default is 16,000.

nosample will not create e(sample), saving some space and speed.

keepsingletons prevents dropping of observations that constitute singleton fixed-effect groups (Correia 2015).

compact temporarily saves all data to disk in orer to free memory for estimation--at the cost of a bit of time.

level(#) sets the confidence level for reported confidence intervals. The default is controlled by set level and is usually 95.

noheader suppresses the display of the table of summary
statistics at the top of the output; only the coefficient table is displayed.
This option is often used in programs and ado-files.

notable suppresses display of the coefficient table.

nofootnote suppresses display of the footnote table that lists the absorbed fixed effects, including the number of categories/levels of each fixed effect,
redundant categories (collinear or otherwise not counted when computing degrees-of-freedom), and the difference between both.

verbose causes reghdfejl to show more of its work--to display the Julia copy of the data set, the formula for the regression model, and the regression command. The data set and formula
are left behind for the user to work with through jl, not erased as they normally are.

Postestimation syntax

Only estat summarize, predict, and test are currently supported.

The syntax of estat summarize and predict is:

estat summarize

Summarizes depvar and the variables described in _b (i.e. not the excluded instruments)

predict
newvar
[if] [in]
[, statistic]

May require you to previously save the fixed effects (except for option xb).
To see how, see the details of the absorb option
Equation: y = xb + d_absorbvars + e

statistic

Main
xb xb fitted values; the default

xbd xb + d_absorbvars

d d_absorbvars

residuals residual

score score; equivalent to residuals

stdp standard error of the prediction (of the xb component)

although predict type newvar is allowed,
the resulting variable will always be of type double.

Examples

. sysuse auto

. reghdfejl price weight length, absorb(rep78)

. reghdfejl price weight length, absorb(rep78) vce(cluster rep78)

. webuse nlswork

. reghdfejl ln_wage age ttl_exp tenure not_smsa south, absorb(idcode year)

. reghdfejl ln_wage age ttl_exp tenure not_smsa south, absorb(year occ_code) cluster(year occ_code)

. reghdfejl ln_wage age ttl_exp tenure not_smsa south, absorb(year occ_code) vce(bs, cluster(occ_code) reps(1000) seed(42) procs(4))

. partialhdfejl ln_wage age ttl_exp tenure not_smsa south, absorb(year occ_code) prefix(_) replace

. reghdfejl _ln_wage _age _ttl_exp _tenure _not_smsa _south, cluster(year occ_code) nocons // same point estimates as in previous regression

. reghdfejl ln_wage (age = ttl_exp tenure) not_smsa south, absorb(year occ_code) cluster(year occ_code) // IV estimation

Stored results

reghdfejl stores the following in e():

Scalars
e(N) number of observations

e(num_singletons) number of singleton observations

e(N_full) number of observations including singletons

e(N_hdfe) number of absorbed fixed-effects

e(tss) total sum of squares

e(tss) total sum of squares after partialling-out

e(rss) residual sum of squares

e(rss) model sum of squares (tss-rss)

e(r2) R-squared

e(r2_a) adjusted R-squared

e(r2_within) Within R-squared

e(r2_a_within) Adjusted Within R-squared

e(df_a) degrees of freedom lost due to the fixed effects

e(rmse) root mean squared error

e(ll) log-likelihood

e(ll_0) log-likelihood of fixed-effect-only regression

e(F) F statistic

e(widstat) For IV regression, Kleibergen-Paap weak identification F

e(rank) rank of e(V)

e(N_clustervars) number of cluster variables

e(N_clust#) number of clusters for the #th cluster variable

e(N_clust) number of clusters; minimum of e(clust#)

e(df_m) model degrees of freedom

e(df_r) residual degrees of freedom

e(sumweights) sum of weights

e(ic) number of iterations

e(converged) 1 if converged, 0 otherwise

e(drop_singletons) 1 if singletons were dropped, 0 otherwise

Macros
e(cmd) reghdfejl

e(cmdline) command as typed

e(depvar) name of dependent variable

e(indepvars) names of independent variables

e(absvars) name of the absorbed variables or interactions

e(clustvar) name of cluster variable

e(clustvar#) name of the #th cluster variable

e(vce) vcetype specified in vce()

e(vcetype) title used to label Std. Err.

e(properties) b V

e(predict) program used to implement predict

e(estat_cmd) program used to implement estat

e(footnote) program used to display footnote

e(marginsnotok) predictions not allowed by margins

e(title) title in estimation output

e(title2) subtitle in estimation output, indicating how many FEs were being absorbed

e(flinejl) Julia command line used to define model formula

e(cmdlinejl) Julia command line used to fit model

Matrices
e(b) coefficient vector

e(V) variance-covariance matrix of the estimators

Functions
e(sample) marks estimation sample

Benchmark

This benchmark creates a data set with 100 million observations and runs regressions with 1 or 2 sets of fixed
effects using areg, reghdfe, and reghdfejl in Stata/MP. The number of processors is set to 1 or 8. THe script is run on a Windows
laptop with an NVIDIA RTX 4070 GPU and an Intel i9-13900H, the latter having 6 performance and 8 (slower) efficiency cores. The log is lightly edited for parsimony.

With 1 set of FE, reghdfe is 2-7 times faster than areg, and reghdfejl is 3-4 faster again, slightly more so with the gpu
option. With 2 sets of FE, reghdfejl is 10 times faster than reghdfe without gpu and 12 times faster with.

On a Mac with an M2 Pro chip--with 8 performance cores and 4 efficiency cores--the absolute times are somewhat higher and the ratios slightly lower (not shown).

Script

scalar N = 100000000

scalar K = 100

set obs `=N'

gen id1 = runiformint(1, N/K)

gen id2 = runiformint(1, K)

drawnorm x1 x2

gen double y = 3 * x1 + 2 * x2 + sin(id1) + cos(id2) + runiform()

set rmsg on

set processors 1

qui areg y x1 x2, a(id1) cluster(id1)

qui reghdfe y x1 x2, a(id1) cluster(id1)

qui reghdfejl y x1 x2, a(id1) cluster(id1)

qui reghdfejl y x1 x2, a(id1) cluster(id1) gpu

qui reghdfe y x1 x2, a(id1 id2) cluster(id1 id2)

qui reghdfejl y x1 x2, a(id1 id2) cluster(id1 id2)

qui reghdfejl y x1 x2, a(id1 id2) cluster(id1 id2) gpu

set processors 8 // requires Stata/MP

qui areg y x1 x2, a(id1) cluster(id1)

qui reghdfe y x1 x2, a(id1) cluster(id1)

qui reghdfejl y x1 x2, a(id1) cluster(id1)

qui reghdfejl y x1 x2, a(id1) cluster(id1) gpu

qui reghdfe y x1 x2, a(id1 id2) cluster(id1 id2)

qui reghdfejl y x1 x2, a(id1 id2) cluster(id1 id2)

qui reghdfejl y x1 x2, a(id1 id2) cluster(id1 id2) gpu

Log

. set processors 1

. qui areg y x1 x2, a(id1) cluster(id1)

t=490.93

. qui reghdfe y x1 x2, a(id1) cluster(id1) dof(none)

t=68.74

. qui reghdfejl y x1 x2, a(id1) cluster(id1)

t=16.86

. qui reghdfejl y x1 x2, a(id1) cluster(id1) gpu

t=15.92

. qui reghdfe y x1 x2, a(id1 id2) cluster(id1 id2) dof(none)

t=315.05

. qui reghdfejl y x1 x2, a(id1 id2) cluster(id1 id2)

t=29.70

. qui reghdfejl y x1 x2, a(id1 id2) cluster(id1 id2) gpu

t=24.48

. set processors 6

. qui areg y x1 x2, a(id1) cluster(id1)

t=99.24

. qui reghdfe y x1 x2, a(id1) cluster(id1) dof(none)

t=44.69

. qui reghdfejl y x1 x2, a(id1) cluster(id1)

t=14.00

. qui reghdfejl y x1 x2, a(id1) cluster(id1) gpu

t=12.23

. qui reghdfe y x1 x2, a(id1 id2) cluster(id1 id2) dof(none)

t=243.90

. qui reghdfejl y x1 x2, a(id1 id2) cluster(id1 id2)

t=27.88

. qui reghdfejl y x1 x2, a(id1 id2) cluster(id1 id2) gpu

t=19.77

Author

David Roodman

Email: david@davidroodman.com

Acknowledgements

More so than for most packages, in writing this one, the author stands on the shoulders of giants. reghdfejl is merely a wrapper
for Matthieu Gomez's FixedEfectModels.jl,
which is itself an implementation---with important innovations---of Sergio Correia's path-breaking reghdfe. reghdfejl's code for
postestimation functionality is copied from reghdfe, as are parts of this help file. The Julia programming language is a free, open-source project.

References

Correia, S. 2015. Singletons, cluster-robust standard rrrors and fixed effects: A bad mixhttps://scorreia.com/research/singletons.pdf.

Correia, S. 2016. A feasible estimator for linear models with multi-way fixed effects. http://scorreia.com/research/hdfe.pdf.