The reason is very similar to the reason that the residual terms will equal zero.

https://stats.stackexchange.com/questions/189584

For a simple linear regression (not for OLS in general which we will show later) you have the formula:

$$y_i = \hat\alpha +  \hat\beta x_{i} + \epsilon_i$$

and for the average this becomes:

$$\bar{y_i} = \alpha + \beta \bar{x_{i}} + \bar{\epsilon_i}$$

And the OLS regression with an intercept included will be such that the residual terms $\bar{\epsilon}$ are equal to zero. 

The reason is that we must have the derivatives of the residual terms equal to zero

$$ \frac{\partial}{\partial \hat\alpha} \sum \epsilon_i^2  =0$$

which will lead to 

$$ \begin{array}{}
\frac{\partial}{\partial \hat\alpha} \sum \epsilon_i^2 &=& \frac{\partial}{\partial \alpha} \sum (y_i - \hat\alpha -  \hat\beta x_{i})^2\\ &=&   \sum -2 \hat\alpha (y_i - \hat\alpha +  \hat\beta x_{i})\\ &=&  -2 \hat\alpha \sum \epsilon_i  &=& 0 \end{array}$$

----------

### Intuitively 

If the line does not go through the point $\bar{x},\bar{y}$ then this will be equivalent to the sum of the errors being unequal to zero, and we will be able to improve the sum of squared residuals $\sum \epsilon_i^2$ by shifting the line up or down.

When the line goes through the point $\bar{x},\bar{y}$, then a change of the parameter $\hat{\alpha}$ does not improve the fit (the derivative is zero).

[![comparison][1]][1]

--------

### It depends on the point of view

The OLS regression curve only passes necessarily through the point $\hat{x},\hat{y}$ when you have *linear* regression. For OLS that is used to express a non-linear function with an intercept term you will still have the situation that the sum of residuals will be equal to zero, but you do not have anymore the relationship $$\bar{y_i} = \alpha + \beta \bar{x_{i}} + \bar{\epsilon_i}$$

An example is a quadratic cuvre

[![quadratic fit][2]][2]

However, you could express this regression of $y$ with a quadratic curve (which is seemingly a non-linear function) as a linear function of $x$ and $x^2$ by means of viewing it as multivariate linear regression. Then it will *also* pass through the means.

 $$\bar{y_i} = \alpha + \beta_1 \overline{x_{i}} + \beta_2 \overline{x_{i}^2} + \bar{\epsilon_i}$$

because now you express it as a multidimensional case.

[![enter image description here][3]][3]



  [1]: https://i.sstatic.net/4VACq.png
  [2]: https://i.sstatic.net/3ppGZ.png
  [3]: https://i.sstatic.net/4S46D.png