https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#In_least_squares_regression_analysis
I was reading this page and came across the fact that the mean of the predicted target values for an OLS regression is always equal to the mean of the original target values. That is, for the set of predicted values $\{\hat{Y}_1, \hat{Y}_2, ...\}$ and the set of original values $\{Y_1, Y_2, ...\}$, the means of the sets are always equal.
Is there a simple proof of why exactly this holds true?
That is, for the set of predicted values $\{\hat{Y}_1, \hat{Y}_2, ...\}$ and the set of original values $\{Y_1, Y_2, ...\}$, the means of the sets are always equal.
The difference between the predicted values and the original values are the residuals
$$\hat{Y}_i = Y_i + r_i$$
So you can write
$$\begin{array}{} \frac{1}{n} \left(\hat{Y}_1+ \hat{Y}_2+ ...\right) &=& \frac{1}{n} \left(({Y}_1 + r_1)+( {Y}_2+r_2)+ ...\right) \\ &=&\frac{1}{n} \left({Y}_1+ {Y}_2+ ...\right)+\frac{1}{n} \left(r_1+ r_2+ ...\right) &=&\frac{1}{n} \left({Y}_1+ {Y}_2+ ...\right) \end{array}$$
and the last equality is true if the method has by design the following property $\left(r_1+ r_2+ ...\right) =0$ and that is the case for OLS. But note that this is only the case when the regression has an intercept term (as Christoph Hanck's answer explains). The residual term is perpendicular to the regressors. If the intercept is one of the regressors (or more generally as jld mentioned in the comments, if it's in the column space of the regressors) then the perpendicularity has as consequence that $\left(r_1,r_2,...\right) \cdot \left(1,1,...\right) = \left(r_1+ r_2+ ...\right) =0$
In simple words you could say that the $\hat{Y}$ are placed equaly in between the $Y$, as much above as below, and that is why they have the same mean.
In matrix notation, the fitted values can be written as $\hat y=Py$, with the projection matrix $P=X(X'X)^{-1}X'$, wich can be verified by plugging in the definition of the OLS estimator into the formula for the fitted values, $\hat y =X\hat\beta$.
Their mean is, with $\iota$ a vector of ones, $$ \iota'Py/n, $$ as the inner product with $\iota$ just sums up elements, $\iota'a=\sum_ia_i$.
In general, we have $PX=X$, as can be verified by direct multiplication.
Now, if $X$ contains $\iota$, i.e., if you have a constant in your regression, we have $P\iota=\iota$, as one of the columns of the result $PX=X$.
Hence, by symmetry of $P$ (which, again, can be verified directly), $$ \iota'Py/n=\iota'y/n, $$ the mean of $y$. Hence, the statement is true if we have a constant in our regression. It is - see the comment by @jld - however also true if there are columns of $X$ that can be combined into $\iota$. That would for example be the case if we have exhaustive dummy variables but no constant (to avoid the dummy variable trap).
A little numerical illustration:
y <- rnorm(20)
x <- rnorm(20)
lm_with_cst <- lm(y~x)
mean(y)
mean(fitted(lm_with_cst))
lm_without_cst <- lm(y~x-1)
mean(fitted(lm_without_cst))
Output:
> mean(y)
[1] 0.04139399
> mean(fitted(lm_with_cst))
[1] 0.04139399
> mean(fitted(lm_without_cst))
[1] 0.05660456
It's intuitively clear. If you have the correct model as the linear regression, the residuals should be distributed with mean zero. If you take the average on the residuals, you are left only with the predicted values.
For example, if your model is
$y = c + ax + \epsilon$,
where $c$ is constant vector, $a$ is coefficient vector, $x$ is the feature vector, $\epsilon$ is the Gaussian residual vector.
When you take the expectation of $y$ for the mean, you get
$E(y) = E(c + ax + \epsilon) = E(c + ax) = E(\hat{y})$
because $E(\epsilon) = 0$ as the mean of residuals are zero.
