A new feature of Stata is the factor variable list. See -help fvvarlist- for more information, but briefly, it allows Stata to create dummy variables and interactions for each observation just as the estimation command calls for that observation, and without saving the dummy value. This makes possible such constructs as interacting a state dummy with a time trend without using any memory to store the 50 possible interactions themselves. (You would still need memory for the cross-product matrix).

There are a large number of regression procedures in Stata that avoid calculating fixed effect parameters entirely, a potentially large saving in both space and time. Where analysis bumps against the 9,000 variable limit in stata-se, they are essential. These are documented in the panel data volume of the Stata manual set, or you can use the -help- command for xtreg, xtgee, xtgls, xtivreg, xtivreg2, xtmixed, xtregar or areg. There are additional panel analysis commands in the SSC mentioned here

However, by and large these routines are not coded with efficiency in mind and will be intolerably slow for very large datasets. Worse still, the -xtivreg2- requires additional memory for the de-meaned data turning 20GB of floats into 40GB of doubles, for a total requirement of 60GB.

-xtreg- is the basic panel estimation command in Stata, but it is very slow compared to taking out means. For example:

What if you have endogenous variables, or need to cluster standard errors? Jacob Robbins has written a fast tsls.ado program that handles those complications:

The dof() option on the -reg- command is used to correct the standard errors for degrees of freedom after taking out means. -distinct- is a very fast way of calculating the number of panel units. I warn you against either of

For taking out means, you may use

For IV regressions this is not sufficient to correct the standard errors. Use the -reg- command for the 1st stage regression. Then run the 2nd stage regression using the predicted (-predict- with the xb option) values for the endogenous variables. In econometrics class you will have learned that the coefficients from this sequence will be unbiased, but the standard errors will be inconsistent. The formulas for the correction of the standard errors are known, and not computationally expensive. An easy way to obtain corrected standard errors is to regress the 2nd stage residuals (calculated with the real, not predicted data) on the independent variables. Those standard errors are unbiased for the coefficients of the 2nd stage regression.

xtreg, tsls and their ilk are good for one fixed effect, but what if you have
more than one? Possibly you can take out means for the largest dimensionality effect
and use factor variables for the others. That works untill you reach the 11,000
variable limit for a Stata regression.
Otherwise, there is **-reghdfe-** on SSC which is an interative process
that can deal with multiple high dimensional fixed effects. It used to be
slow but I recently tested a regression with a million observations and
three fixed effects, each with 100 categories. That took 8 seconds
(limited to 2 cores). Increasing the number of categories to 10,000
only tripled the execution time.