Linear probability model, dummy variables and the same standard errors on all estimates I am fitting the linear probability model,
$$
Y_i=\beta_0 + \sum_{j=1}^J \beta_j G_{ji} +\varepsilon_i
$$
where $Y_i \in \{0,1\}$ and $G_{ji} \in \{0,1\}$, for $j=1,\ldots,J$ and $\sum_{j=1}^J G_{ji}=1$, i.e., the $G_{ij}$ are groups and each individual is placed into one of these groups. 
I use OLS to get estimates for $\beta_j$, $j=1,\ldots,J$ (leaving out one group to avoid multi-linearity) and I get my estimates. The weird thing is that the standard errors are the same for all $\beta_j$, $j=1,\ldots,J$, i.e. $se(\widehat \beta_1)=se(\widehat \beta_{2})=\ldots= se(\widehat \beta_J)$. 
Is this normal that the standard errors are the same? It's not the first time that this has happened to me and I'm not sure if I'm doing something wrong. 
 A: Recall that the standard errors are the diagonal elements of the matrix
$$
\hat\sigma^2(X'X)^{-1}
$$
As pointed out by @repmat, this result requires that each group is of equal size, i.e., that 
$$\sum_iG_{ji}=c$$
for $j=1,\ldots,J$.
In that case, you can easily check that 
$$
X'X=n
\begin{pmatrix}
1&1/J&\cdots&\cdots&\cdots&1/J\\
1/J&1/J&0&\cdots&\cdots&0\\
1/J&0&1/J&\ddots&\cdots&0\\
\vdots&0&0&\ddots&\ddots&0\\
\vdots&\vdots&\cdots&\ddots&\ddots&0\\
1/J&0&\cdots&\cdots&0&1/J\\
\end{pmatrix},
$$
assuming that the first column contains the constant term. Assuming that one of the redundant dummies has been dropped, so that $X'X$ is of dimension $J\times J$, the inverse is
$$
(X'X)^{-1}=\frac{1}{n}\begin{pmatrix}
J&-J&\cdots&&\cdots&-J\\
-J&2J&J&\cdots&\cdots&J\\
\vdots&J&2J&J&&J\\
&\vdots&J&2J&\ddots&\vdots\\
\vdots&\vdots&&\ddots&2J&J\\
-J&J&\cdots&\cdots&J&2J
 \end{pmatrix}
$$
We see that the diagonal elements are identical (except for the one on the constant), yielding identical standard errors.
The nature of the $y_i$ is irrelevant, as these enter the standard errors only through $\hat\sigma^2$, which however obviously just multiplies each element of $(X'X)^{-1}$ and hence does not affect the result that the diagonal elements are identical.
This result also shows that the squared standard error on the constant is one half that of the dummies.
P.S.: As suggested by Alecos in the comments, we may define the block matrix
$$
X'X=
n\begin{pmatrix}
A&B\\C&D
\end{pmatrix},
$$
with $A=1$, $B=(1/J,\ldots, 1/J)$, $C=B'$ and $D=I/J$ and use a result on partitioned inverses that the lower right block has the inverse
$$
D^{-1}+D^{-1}C(A-BD^{-1}C)^{-1}BD^{-1}
$$
to see that the result is as presented as above.
UPDATE: Regarding the discussion in the comments of @repmat's answer, the numerical equivalence does not exactly hold for robust standard errors. 
This is because the "meat" matrix $X'\Sigma_uX$ of the robust variance estimator
$$
(X'X)^{-1}X'\Sigma_uX(X'X)^{-1}
$$
has a diagonal 
$$
\begin{pmatrix}
\sum_{i=1}^n\hat{u}_i^2\\
\sum_{i=1,\,i\in j=1}^n\hat{u}_i^2\\
\vdots\\
\sum_{i=1,\,i\in j=J-1}^n\hat{u}_i^2\\
\end{pmatrix}
$$
(assuming that group $J$ has been dropped to avoid multicollinearity), and there is in general no reason to believe that the sums of squared residuals belonging to the different groups are identical.
The differences are minor, though (at least under homoskedasticity). Here is a modification of his illustration:
set.seed(42)
n <- 999
library(sandwich)
library(lmtest)

year1 <- data.frame(rep(1, n/3))
year2 <- data.frame(rep(2, n/3))
year3 <- data.frame(rep(3, n/3))

require(plyr)
years <- rbind.fill(year1, year2, year3)
years[is.na(years)] <- 0
years <- as.factor(rowSums(years))

y <- round(runif(n),0)
reg <- lm(y ~ years)

coeftest(reg, vcov = vcovHC(reg, "HC0"))

> 

t test of coefficients:

             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.486486   0.027390 17.7616   <2e-16 ***
years2      -0.027027   0.038678 -0.6988   0.4849    
years3      -0.012012   0.038717 -0.3103   0.7564    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

A: This will happen if, and only if, the two (or more variables) have the same variance, or in other words that all groups are equally large (in terms of 1's). The nature of $y$ does not matter. 
Here is an example from R:
set.seed(42)
year1 <- data.frame(rep(1, 333))
year2 <- data.frame(rep(2, 333))
year3 <- data.frame(rep(3, 333))

require(plyr)
years <- rbind.fill(year1, year2, year3)
years[is.na(years)] <- 0
years <- as.factor(rowSums(years))

y <- round(runif(999),0)
coef(summary(lm(y ~ years)))


               Estimate Std. Error    t value     Pr(>|t|)
(Intercept)  0.48648649 0.02739570 17.7577679 1.667287e-61
years2      -0.02702703 0.03874337 -0.6975910 4.855958e-01
years3      -0.01201201 0.03874337 -0.3100404 7.565951e-01

If the groups are not equally large (has the same numbers of 1's), the errors will be similar but not identical (as above).
The intuition is that you are just estimating means in different groups. Two dummy variables will always have the same variance if the amount of 1's are the same. Because the formula for the variance is:
$$
p (1-p)
$$
Where $p$ is the (sample) mean of the variable. 
Edit: To be clear, you can add other non dummy variables to the regression above. The two dummy variables will still have the same variance.
A: You could ask yourself the question 'Why would I expect one of them to be much larger than the others?'. I do not think you are doing anything wrong as far as specifying your model matrix is concerned. I do wonder though whether using linear regression for an outcome which is binary is a good thing. Most people would use logistic regression here.
