Standard error of the intercept in Frisch-Waugh theorem (de-meaned regression) I am applying in Frisch-Waugh Theorem to partial out a set of fixed effects D and get the OLS estimates and standard errors of the remaining regressors X.
The theorem is more general, but one leading application is to re-center the outcome and the right-hand side variables about group means ("de-meaning"), hence absorbing the group fixed-effects without having to estimate them
(see for instance the discussion here).
Below a MWE in Stata/Mata to clarify what I do.
If you have a solution with another software, it is well accepted, I am mainly interested in the theory behind it. 
Here, y is price, D is turn, and X corresponds to the remining right-hand side variables (including other fixed effects, trunk) and the constant. 
For reference I report also a benchmark when using -areg-.  
cls
clear all
sysuse auto, clear


// Benchmark

areg price gear length i.trunk, absorb(turn) 


// Absorb "manually" in MATA

xi i.turn i.trunk
gen uno = 1

mata

 // import
    y = st_data(., "price")
    X = st_data(., ("gear", "length", "_Itrunk_*", "uno"))  
    D = st_data(., "_Iturn*")

 // demeaned X and y
    M_D = I(rows(y)) - D * qrinv(cross(D,D)) * D'  // "residual maker" M_D = I - D(D'D)^(-1)D' 
    y_dem = M_D * y
    X_dem = M_D * X

 // OLS using de-meaned variables and corresponding standard errors
    b1 = qrsolve(X_dem, y_dem)
    res2 = cross(y_dem - X_dem*b1, y_dem - X_dem*b1)  
    MSE = res2/(rows(X_dem) - (cols(X_dem)+cols(D)-1))
    XX_dem = qrinv(cross(X_dem, X_dem))
    SE = sqrt(diagonal(XX_dem) * MSE)
    SE

 // Compute constant as c = y_bar - X_bar*b1
    e = rows(b1) - 1
    c = mean(y) - (mean(X[., 1..e]) * b1[1..e])'
    c

end

All the coefficients (including the separately computed constant) correspond to the ones reported by -areg-. Also the standard errors are correct, except for the one relative to the constant, which I am not sure how to get.
Any help is highly appreciated.
 A: When de-meaning the model, the intercept is obviously lost. 
If we are interested in estimating the intercept, we need to proceed in two steps.$^1$
First, in order to get the intercept back, one can add $\bar{y}$, the grand total mean of $y$, on the left-hand side of the demeaned model equation, and the intercept $\alpha$ and $\bar{X}'\beta + \bar{u}$ on the right-hand side of the demeaned model. 
It's trivial to see that, by construction, this new equation is equivalent to the demeaned one (since we're adding the same quantity on both sides of the equation).
However, with this trick we now have the intercept in the model, and hence by doing OLS we can estimate it (together with $b$). We also directly get the full variance-covariance matrix and the correct sum of squared residuals.
This procedure is explained more in detail here (see in particular the Derivation section).
As a "sanity check" one can always subtract $\bar{X}'b$ (without the constant) from $\bar{y}$, as I do in the MWE. One should get exactly the same intercept, which is the case.   
Second, in order to get the correct standard errors, we must take into account that we demeaned the variables at the first step, hence we must adjust the degrees of freedom (see e.g. the discussion in this post, and here). This is all we need to do to get all correct standard errors, including that of the constant. 
Since here only one fixed effect is partialled-out, the intercept that is added back to the model is interpreted as the average fixed effect across all groups of the partialled-out variable.
Finally, note that the overall approach is consistent with the rationale behind the FW theorem: we partial out a set of covariates (used to de-mean) and get the OLS estimates and standard errors of the remaining ones without ever having to estimate the former ones. This is particularly useful (sometimes needed) when we want to adjust/control for a set of covariates but we are not really interested in their OLS estimates, which is often the case when we have hundreds or perhaps thousands of fixed effects.
For a more code-based answer, see the cross-posted question here. 

$^1$ This motivated my question. However, strictly speaking, we do not need to estimate the intercept to get the correct OLS estimates and standard errorrs of the non-absorbed covariates: we can just estimate the demeaned model as I do in the MWE.
