To see equality, let us first derive the FE estimator.
Define the residual-maker matrix
\begin{align*}
\underset{(M\times M)}{\mathbf{Q}}&:=\mathbf{I}_M-\mathbf{1}_M(\mathbf{1}_M'\mathbf{1}_M)^{-1}\mathbf{1}_M'\\
&=\mathbf{I}_M-\left(%
\begin{array}{ccc}
1/M & \cdots & 1/M \\
\vdots & \ddots & \vdots \\
1/M & \cdots & 1/M \\
\end{array}%
\right)\mathbf{1}_M\mathbf{1}_M', \end{align*}
where $M$ denotes the number of observations per individual unit in the panel.
Premultiplication with $\mathbf{Q}$ centers the $\mathbf{y}_i$ and $\mathbf{Z}_i$ around their averages over $m$,
\begin{align*}
\mathbf{Q}\mathbf{y}_i&=\mathbf{y}_i-\mathbf{1}_M\mathbf{1}_M'\mathbf{y}_i/M\\&=\mathbf{y}_i-\mathbf{1}_M\overline{y_{i}}.
\end{align*}
The also implies that every time invariant variable from the set of regressors $\mathbf{Z}_i$ turns into a column of zeros, and hence is eliminated from the data.
This is a serious disadvantage of the FE estimator. Consider the example of wage regressions for a panel of employees. Variables such as gender or schooling are of primary interest, but (typically) do not change over time (anymore).
As $\mathbf{Q}\mathbf{1}_M=\mathbf{0}$, we have that, using the error-component model $\mathbf{y}_i=\mathbf{Z}_i\mathbf{\delta}+\mathbf{1}_M\alpha_i+\mathbf{\eta}_{i}$, where $\eta_i$ denotes the $M$-vector of idiosyncratic time-varying errors,
\begin{align*}
\mathbf{Q}\mathbf{y}_i&=\mathbf{Q}\mathbf{F}_i\mathbf{\beta}+\mathbf{Q}\mathbf{\eta}_{i}\qquad i=1,\ldots,n\\
\tilde{\mathbf{y}}_i&\equiv\tilde{\mathbf{F}}_i\mathbf{\beta}+\tilde{\mathbf{\eta}}_{i},
\end{align*}
where $\mathbf{F}_i$ is the $(M\times L_b)$-matrix of the observations on the time variant regressors. Stacking the observations over the $n$ units gives
$$
\underset{(Mn\times 1)}{\tilde{\mathbf{y}}}:=\left(%
\begin{array}{c}
\tilde{\mathbf{y}}_1 \\
\vdots \\
\tilde{\mathbf{y}}_n \\
\end{array}%
\right)\qquad\underset{(Mn\times L_b)}{\tilde{\mathbf{F}}}:=\left(%
\begin{array}{c}
\tilde{\mathbf{F}}_1 \\
\vdots \\
\tilde{\mathbf{F}}_n \\
\end{array}%
\right)
$$
The FE estimator is simply OLS applied to these $Mn$ observations:
\begin{align*}
\widehat{\mathbf{\beta}}_{\text{FE}}&=(\tilde{\mathbf{F}}'\tilde{\mathbf{F}})^{-1}\tilde{\mathbf{F}}'\tilde{\mathbf{y}}
\end{align*}
To see the equality between FE and least squares dummy variables, stack the observations a bit further:
\begin{equation}
\underset{(Mn\times 1)}{\mathbf{y}}:=\left(%
\begin{array}{c}
\mathbf{y}_1 \\
\vdots \\
\mathbf{y}_n \\
\end{array}%
\right)\;\underset{(Mn\times L_b)}{\mathbf{F}}:=\left(%
\begin{array}{c}
\mathbf{F}_1 \\
\vdots \\
\mathbf{F}_n \\
\end{array}%
\right)
\end{equation}
and
\begin{equation}
\underset{(Mn\times 1)}{\mathbf{\eta}}:=\left(%
\begin{array}{c}
\mathbf{\eta}_1 \\
\vdots \\
\mathbf{\eta}_n \\
\end{array}%
\right)\;
\underset{(n\times 1)}{\mathbf{\alpha}}:=\left(%
\begin{array}{c}
\alpha_1 \\
\vdots \\
\alpha_n \\
\end{array}%
\right).
\end{equation}
Further, let
$$
\underset{(Mn\times n)}{\mathbf{D}}:=\mathbf{I}_n\otimes\mathbf{1}_M=\left(%
\begin{array}{ccc}
\mathbf{1}_M & & \mathbf{O} \\
& \ddots & \\
\mathbf{O}& & \mathbf{1}_M \\
\end{array}
\right)
$$
Then, the linear panel data model from under an error component assumption in matrix notation is obtained as
$$
\mathbf{y}=\mathbf{D}\mathbf{\alpha}+\mathbf{F}\mathbf{\beta}+\mathbf{\eta},
$$
a dummy-variable model.
That is, we can also obtain an estimator of $\mathbf{\beta}$ from an OLS regression on the regressors and $n$ individual specific effects.
Now, note that the Frisch-Waugh-Lovell Theorem says that the OLS estimator of $\mathbf{\beta}$ can be found by regressing $\mathbf{M}_{\mathbf{D}}\mathbf{y}$ on $\mathbf{M}_{\mathbf{D}}\mathbf{F}$, where
$$\underset{(Mn\times Mn)}{\mathbf{M}_{\mathbf{D}}}:=\mathbf{I}-\mathbf{D}(\mathbf{D}'\mathbf{D})^{-1}\mathbf{D}'$$
Using symmetry and idempotency of $\mathbf{M}_{\mathbf{D}}$ gives
\begin{equation}
\widehat{\mathbf{\beta}}_{\text{LSDV}}=(\mathbf{F}'\mathbf{M}_{\mathbf{D}}\mathbf{F})^{-1}\mathbf{F}'\mathbf{M}_{\mathbf{D}}\mathbf{y}
\end{equation}
Now,
\begin{align*}
\mathbf{M}_{\mathbf{D}}&=\mathbf{I}_{Mn}-(\mathbf{I}_n\otimes\mathbf{1}_M)[(\mathbf{I}_n\otimes\mathbf{1}_M)'(\mathbf{I}_n\otimes\mathbf{1}_M)]^{-1}(\mathbf{I}_n\otimes\mathbf{1}_M)'\\
&=\mathbf{I}_{n}\otimes\mathbf{I}_{M}-(\mathbf{I}_n\otimes\mathbf{1}_M)[(\mathbf{I}_n\otimes\mathbf{1}_M')(\mathbf{I}_n\otimes\mathbf{1}_M)]^{-1}(\mathbf{I}_n\otimes\mathbf{1}_M')\\
&=\mathbf{I}_{n}\otimes\mathbf{I}_{M}-(\mathbf{I}_n\otimes\mathbf{1}_M)[\mathbf{I}_n\otimes\mathbf{1}_M'\mathbf{1}_M]^{-1}(\mathbf{I}_n\otimes\mathbf{1}_M')\\
&=\mathbf{I}_{n}\otimes\mathbf{I}_{M}-(\mathbf{I}_n\otimes\mathbf{1}_M)[\mathbf{I}_n\otimes M]^{-1}(\mathbf{I}_n\otimes\mathbf{1}_M')\\
&=\mathbf{I}_{n}\otimes\mathbf{I}_{M}-(\mathbf{I}_n\otimes\mathbf{1}_M)\left[\mathbf{I}_n\otimes \frac{1}{M}\right](\mathbf{I}_n\otimes\mathbf{1}_M')\\
&=\mathbf{I}_{n}\otimes\mathbf{I}_{M}-(\mathbf{I}_n\otimes\mathbf{1}_M)\left[\mathbf{I}_n\otimes \frac{1}{M}\mathbf{1}_M'\right]\\
&=\mathbf{I}_{n}\otimes\mathbf{I}_{M}-\mathbf{I}_n\otimes\mathbf{1}_M\frac{1}{M}\mathbf{1}_M'\\
&=\mathbf{I}_{n}\otimes\left(\mathbf{I}_{M}-\frac{1}{M}\mathbf{1}_M\mathbf{1}_M'\right)\\
&=\mathbf{I}_n\otimes\mathbf{Q}
\end{align*}
Thus,
\begin{align*}
\mathbf{M}_{\mathbf{D}}\mathbf{F}&=(\mathbf{I}_n\otimes\mathbf{Q})\mathbf{F}\\
&=\left(%
\begin{array}{ccc}
\mathbf{Q} & & \\
& \ddots & \\
& & \mathbf{Q} \\
\end{array}
\right)\mathbf{F}\\
&=\tilde{\mathbf{F}},
\end{align*}
so that $$\widehat{\mathbf{\beta}}_{\text{LSDV}}=\widehat{\mathbf{\beta}}_{{FE}}.$$
Incidentally, while the notation works with balanced panel data, the result also goes through in the unbalanced case, as one can either check with more complicated notation or this numerical illustration:
library(plm)
# panel dimensions
n <- 10
m <- sample(2:4, n, replace=T) # unbalanced panel
# some data
alpha <- runif(n)
beta <- -2
y <- X <- y.d <- X.d <- c()
D <- matrix(0, sum(m), n) # for the dummy variable matrix
row.counter <- 0
for (i in 1:n) {
X.n <- runif(m[i],i,i+1)
X.d <- c(X.d, X.n - mean(X.n))
X <- c(X,X.n)
y.n <- alpha[i] + X.n*beta + rnorm(m[i])
y <- c(y, y.n)
y.d <- c(y.d, y.n - mean(y.n))
D[(row.counter+1):(row.counter+m[i]), i] <- rep(1, m[i])
row.counter <- row.counter + m[i]
}
Output:
> # plm
> paneldata <- data.frame(rep(1:n, times=m), unlist(sapply(m, function(i) 1:i)), y, X) # first two columns are for plm to understand the panel .... [TRUNCATED]
> FE <- plm(y~X, data = paneldata, model = "within")
> # results:
> coef(FE) # the slope coefficient
X
-2.331847
> fixef(FE) # the intercepts
1 2 3 4 5 6 7 8 9 10
0.99396 2.30328 1.90957 2.22670 1.09438 3.10411 2.03265 4.39759 4.42384 4.15294
> # FWL
> lm(y.d~X.d-1) # just the slope in this formulation
Call:
lm(formula = y.d ~ X.d - 1)
Coefficients:
X.d
-2.332
> # LSDV
> lm(y~D+X-1) # intercepts and slope
Call:
lm(formula = y ~ D + X - 1)
Coefficients:
D1 D2 D3 D4 D5 D6 D7 D8 D9 D10 X
0.994 2.303 1.910 2.227 1.094 3.104 2.033 4.398 4.424 4.153 -2.332
Yes including both violates certain statistical properties. The pggls documentation indirectly states exactly that:
Conversely, this structure is assumed identical across groups and thus
general FGLS estimation is inefficient under groupwise
heteroskedasticity. Note also that this method requires estimation of
T(T+1)/2 variance parameters, thus efficiency requires N > > T (if
effect="individual", else the opposite).
FGLS requires estimator of covariance matrix of regression disturbances. For individual effects panel data model:
$$y_{it} = x_{it}\beta + c_i + u_{it},$$
where $c_{i}$ is an individual effect, it is assumed that $u_{it}$ are independent from $u_{jt}$ for each $i\neq j$, so you are left with $\frac{T(T+1)}{2}$ covariances $\text{cov}(u_{it},u_{is})$.
For the time effects panel data model:
$$y_{it} = x_{it}\beta + d_t + u_{it},$$
where $d_t$ is a time effect it is assumed that $u_{it}$ are independent from $u_{is}$ for each $t\neq s$, so you are left with $\frac{N(N+1)}{2}$ covariances $\text{cov}(u_{it},u_{is})$.
Now if you have both time and individual effects:
$$y_{it} = x_{it}\beta + c_i + d_t + u_{it},$$
the question arises which covariances $\text{cov}(u_{it},u_{js})$ are zero? If you assume that all of them are not zero, you are left with $\frac{NT(NT+1)}{2}$ unknown parameters with $NT$ data points, which makes the problem non feasible.
Note 1. Independence assumption mention can be relaxed to zero covariances.
Note 2. In both individual and time effect there $NT$ data points. FGLS is an asymptotic procedure and it requires that the number of data points must increase, while the number of parameters remains fixed. For the individual effects hence the $N$ must increase, and for time effects $T$ must increase. More often than not these requirements are satisfied for totally different data sources, hence the answer to your problem depends on your data source. Is is more likely that $N$ is increasing or $T$? Since you mention time dummies, I suspect that $N$ is increasing, hence I suggest using time dummies.
Best Answer
Ok, found this in internet:
"Time dummy is a variable which equals 1 for a given year and 0 for all other years. It allows to control for time-specific fixed effects i.e. shocks which impact is restricted to a given time-period, affects or panel units and are not controlled by other explanatory variables. " source: https://www.researchgate.net/post/Is_anyone_familiar_with_Time_Trends_vs_Time_Dummies
And: " Year effects (more simply known as “year dummies”" Source: https://www.dartmouth.edu/~ethang/Lectures/Class17/Always%20Control%20for%20Year%20Effects%20in%20Panel%20Regressions.pdf
So i think the use of Time dummies, are Time effects! But if time effects are always fixed I dont know.