Why do residuals in linear regression always sum to zero when an intercept is included? I'm taking a course on regression models and one of the properties provided for linear regression is that the residuals always sum to zero when an intercept is included. 
Can someone provide a good explanation for why this is the case? 
 A: When an intercept is included in multiple linear regression,
$$
\hat{y}_i = \beta_0 + \beta_1x_{i,1} + \beta_2x_{i,2} +…+ \beta_px_{i,p}
$$
In Least squares regression, the sum of the squares of the errors is minimized.
$$
SSE=\displaystyle\sum\limits_{i=1}^n \left(e_i \right)^2= \sum_{i=1}^n\left(y_i - \hat{y_i} \right)^2= \sum_{i=1}^n\left(y_i -\beta_0- \beta_1x_{i,1}-\beta_2x_{i,2}-…- \beta_px_{i,p} \right)^2
$$
Take the partial derivative of SSE with respect to $\beta_0$ and setting it to zero.
$$
\frac{\partial{SSE}}{\partial{\beta_0}} = \sum_{i=1}^n 2\left(y_i -\beta_0- \beta_1x_{i,1}-\beta_2x_{i,2}-…- \beta_px_{i,p} \right)^1 (-1) =-2\displaystyle\sum\limits_{i=1}^ne_i=0
$$
Hence, the residuals always sum to zero when an intercept is included in linear regression.
A: Write the linear model in matrical form as
\begin{align}
y = X\beta + \varepsilon,  
\end{align}
where $y \in \mathbb{R}^n$ is the response vector, $X \in \mathbb{R}^{n \times p}$ is the design matrix, $\varepsilon \in \mathbb{R}^n$ is the error vector.
A key observation is that because the model has intercept, $1$, which is the first column of design matrix $X$, can be written as
$$1 = Xe,$$
where $e$ is a column vector with all zeros but the first component one. Also note, in matrix notation, the sum of residuals is just $1^T(y - \hat{y})$, where $\hat{y} = Hy$ and $H$ is the projection matrix (i.e., "hat matrix") equals to $X(X^TX)^{-1}X^T$.
Therefore,
\begin{align}
& 1^T(y - \hat{y}) = 1^T(I - H)y \\
= & e^TX^T(I - X(X^TX)^{-1}X^T)y \\
= & e^T(X^T - X^TX(X^TX)^{-1}X^T)y \\
= & e^T(X^T - X^T)y \\
= & 0.
\end{align}
A: This follows directly from the normal equations, i.e. the equations that the OLS estimator solves,
$$\mathbf{X}^{\prime} \underbrace{\left( \mathbf{y} - \mathbf{X} \mathbf{b} \right)}_{\mathbf{e}} = 0 $$
The vector inside the parentheses is of course the residual vector or the projection of $\mathbf{y}$ onto the orthogonal complement of the column space of $X$, if you like linear algebra. Now including a vector of ones in the $\mathbf{X}$ matrix, which by the way doesn't have to be in the first column as is conventionally done, leads to
$$\mathbf{1}^{\prime} \mathbf{e} = 0 \implies \sum_{i=1}^n e_i = 0$$ 
In the two-variable problem this is even simpler to see, as minimizing the sum of squared residuals brings us to
$$\sum_{i=1}^n \left(y_i - a - b x_i \right) = 0$$
when we take the derivative with respect to the intercept. From this then we proceed to obtain the familiar estimator
$$a = \bar{y} - b \bar{x}$$
where again we see that the construction of our estimators imposes this condition.
A: In case you are looking for a rather intuitive explanation.
In some sense, the linear regression model is nothing but a fancy mean. To find the arithmetic mean $\bar{x}$ over some values $x_1, x_2, \dots, x_n$, we find a value that is a measure of centrality in a sense that the sum of all deviations (where each deviation is defined as $u_i = x_i - \bar{x}$) to the right of the mean value are equal to the sum of all the deviations to the left of that mean. There is no inherent reason why this measure is good, let alone the best way to describe the mean of a sample, but it is certainly intuitive and practical. The important point is, that by defining the arithmetic mean in this way, it necessarily follows that once we constructed the arithmetic mean, all deviations from that mean must sum to zero by definition!
In linear regression, this is no different. We fit the line such that the sum of all differences between our fitted values (which are on the regression line) and the actual values that are above the line is exactly equal to the sum of all differences between the regression line and all values below the line. Again, there is no inherent reason, why this is the best way to construct a fit, but it is straightforward and intuitively appealing. Just as with the arithmetic mean: by constructing our fitted values in this way, it necessarily follows, by construction, that all deviations from that line must sum to zero for otherwise this just wouldn't be an OLS regession. 
A: A simple derivation using matrix algebra:
$\sum e $ can be written as $1^Te$
Then 
$1^Te = 1^T(M_x y)$ where $M_x$ is the orthogonal matrix. Since $M_x$ is symmetric we can rearrange so that
$(M_x1)^Ty$
which equals zero if $M_x$ and $1$ are orthogonal, which is the case if the matrix of the regressors $x$ contains the intercept (a vector of $1$, indeed).
A: *

*$e_i = y_i - [1, X] [a, b] = y_i - Xb - a = v_i - a$

*$\frac{d}{da} \sum e_i^2 \propto \sum e_i\cdot 1 = \sum v_i - a = 0$ so $\hat{a} = \frac{1}{n}\sum v_i$

*$\sum e_i = \sum_i v_i - a = \sum_i v_i - \frac{n}{n}\sum_i v_i = 0$
..
