One of the reasons why this feature of
glm is useful is the possibility to perform quasi-maximum likelihood estimations. I cannot be be sure if this was originally the main purpose for not restricting the admissible domain of the dependent variable, but it gives you a good example of a setting where this is very useful.
See, as a similar case, the fractional regression, where the dependent variable is continuous in [0,1] instead of binary but the link function is often chosen to be a normal cdf or a logistic function. As long as the likelihood belongs to the linear exponential family and the range of variation of the dependent variable is the same, the parameters of the conditional mean are consistently estimated (even if the distribution of the dependent variable is misspecified; see GMT(1984)).
fracreg [logit|probit] has been around for a while now and performs this estimation, but an alternative is just
glm y X, link(logit) family(binomial) vce(robust) (and when Papke & Wooldridge (1996) was published the only way in Stata; see Baum(2008)). In both cases the estimation is by quasi-maximum likelihood.
See also at the link that Nick Cox posted in a comment above:
"It turns out that the estimated coefficients of the maximum-likelihood Poisson estimator in no way depend on the assumption that E(yj) = Var(yj), so even if the assumption is violated, the estimates of the coefficients b0, b1, …, bk are unaffected." That's exactly what I am talking about.
Baum (2008) Stata Tip 63: Modeling Proportions, The Stata Journal
Goriereux, Monfort, Trognon (1984), Pseudo Maximum Likelihood Methods: Theory
Papke, Wooldridge (1996), Econometric Methods for Fractional Response Variables With an Application to 401 (K) Plan Participation Rates