Demeaning with two (n) fixed effects in panel regressions I would like to ask for your help concerning the following issue.
In a nutshell, I would like to know, how demeaning works in a panel regression with two (separate) fixed effects.
Let us assume I have the textbook balanced panel with $N=100$ individuals $i$ and $T=20$ time periods $t$ (say, 100 firms observed over 20 years, neither does a firm enter the sample, nor does a firm leave it). I observe regressand $Y$ and the two regressors $X$ and $Z$.
To estimate a two-way-fixed effects model with both individual- and time-fixed effects, I could estimate (being only interested in coefficients $\mathbf{\beta}$)
$Y_{i,t} = \lambda^i + \lambda^t + \beta_{1} \cdot X_{i,t}+ \beta_{2}\cdot Z_{i,t}+\varepsilon_{i,t},$
where $\lambda^{m}$ is meant as a short-hand notation for $\sum_{m=1}^M \lambda_m \cdot \mathbf{I}(m \in group~ m)$, with $\mathbf{I}(.)$ being the indicator function, which takes upon the value of 1, if the expression in brackets is true and 0 else. I.e., I include a dummy-variable for every individual $i$ and time period $t$.
But for computational convenience, I would like to get rid off the dummies and estimate the model by demeaning. 
However, I am stuck which means to subtract and which terms to add.
Thus, let us visit the easier example of exclusively individual fixed effects.
Instead of estimating
$Y_{i,t} = \lambda^i + \alpha_{1} \cdot X_{i,t}+ \alpha_{2}\cdot Z_{i,t}+\eta_{i,t}$
I estimate 
$Y_{i,t}-\bar{Y}_i = \alpha_{1} \cdot (X_{i,t}-\bar{X}_i)+ \alpha_{2} \cdot (Z_{i,t}-\bar{Z}_i) + (\eta_{i,t}-\bar{\eta}_{i}),$
where "upper-bar with index i" denotes means taken over i. E.g., for $i=1$, for each of the $T=20$ observations I have for $i=1$, I subtract the mean of $i=1$ computed over that $T=20$ observations.
My question is:
With two separate fixed effects: Which means do I subtract? What do I add back? E.g., simply additionally subtracting an $\bar{Y}_t$ for $Y$ would presumably not do the trick.
Is there a formula for demeaning wich generalized to the case of integer $f>0$ fixed effects?
Yours sincerely,
Sinistrum
 A: If you use demean approach (which is theoretically right), then you have to do demean your data both cross sectionally and time series (irrespective of the order). See how it is works. 
Assume following regression model: 
$$y_{it} = u_i + \nu_t + \beta X_{it} + e_{it} \,\,\,\,\, i = 1, 2, \dots, n \,\,\,\,\, T = 1, 2, \dots, t \tag{1}$$ 
First demean cross sectionally. Mean equation is
$$\bar{y_i} = u_i + \bar{v} + \beta \bar{X_i} + \bar{e_i}\,\,\,\,\, i = 1, \dots, n \tag{2}$$
(notice, the $\bar{v}$ will be same for each cross sectional) 
Subtracting (1) - (2)
$$y_{it} - \bar{y_i} = \nu_t - \bar{v} + \beta (X_{it} - \bar{X_i}) + (e_{it} - \bar{e_i}) \tag{3}$$
(see there is no fixed effect in equation 3)
Now, take mean of equation 1 for each $t$, and the mean equation is, 
$$\bar{y}_t = \bar{u} + \nu_t + \beta \bar{X_t} + \bar{e_t} \,\,\,\,\, T = 1, 2, \dots, t\tag{4}$$
Now, subtract equation 4 from equation 3, we get: 
$$y_{it} - \bar{y_i}  - \bar{y_t} = \beta (X_{it} - \bar{X}_i - \bar{X}_t)  -(\bar{v} + \bar{u}) + (e_{it} - \bar{e_i} - \bar{e_t})  \tag{5}$$
In this way, there is no fixed effects and time effects in equation 5. 
A: Been trying to figure this out myself and saw the only other answer to this question is wrong here. You need to subtract time and group means but then add the overall mean back in. See Greene (2012) on Fixed Time and Group Effects (section 11.4.4). You can try it out yourself and see that just subtracting the time and group means does not give the correct result.
A: The point is to rid the process of the individual-specific and time-specific nuisance parameters (in this case the $\alpha_i$ and $\theta_t$ in notation below). This idea dates back all the way to Neyman and Scott's (1948) Econometrica paper. To do so, you should subtract the individual-specific and time-specific means and add back the grand mean, as follows (suppose $x_{it}$ is scalar without loss of generality):
$y_{it} = \alpha_i + \theta_t + \beta x_{it} + u_{it}$
$\bar{y}_{i}:=\frac{1}{T}\sum_{t=1}^Ty_{it} = \alpha_i + \bar{\theta} + \beta \bar{x}_{i} + \bar{u}_{i}$
$\bar{y}_{t}:=\frac{1}{N}\sum_{i=1}^Ty_{it} = \bar{\alpha} + \theta_t + \beta \bar{x}_{t} + \bar{u}_{t}$
$\bar{y}:=\frac{1}{NT}\sum_{t=1}^T\sum_{i=1}^N y_{it} = \bar{\alpha} + \bar{\theta} + \beta \bar{x} + \bar{u}$,
now, define $\dot{y}_{it}:= y_{it}-\bar{y}_{i}-\bar{y}_{t} + \bar{y}$ to obtain an equation which can be readily used to estimate $\beta$ under suitable restrictions which I won't state here as it is beyond OP's scope:
$\dot{y}_{it} = \beta \dot{x}_{it} + \dot{u}_{it}$
The current top answer isn't really appropiate because the nuisance parameters are still in the model, that is why adding back the grand mean is important. For more modern references, see standard graduate texts on panel data, e.g., Analysis of Panel Data by Hsiao pg. 62, Econometric Analysis of Panel Data by Baltagi pg. 27-28.
A: There also exists an iteration procedure using the equation shown by Neeraj (Source: https://doi.org/10.1177/1536867X1501500318). This procedure seems also to lead to the same results if there is an unbalanced panel. I assume that the formula shown by Greene is also subject to some restriction as adressed by Helix123.
