In order to prove the Gauss-Markov Theorem, we first have to show that the OLS estimate $\hat{\theta}$ is an unbiased estimator. From what Im reading on Internet and some textbooks, these are the main steps (skipping some stages):
\begin{equation} \begin{aligned} \mathbb{E}[\hat{\theta}] = \quad & \mathbb{E}[(X^TX)X^TY] \\ \stackrel{\ast}{=} \quad & \mathbb{E}[(X^TX)X^T(X\theta + \epsilon)] = \theta \end{aligned} \end{equation}
But I can't figure out why on the $\stackrel{\ast}{=}$ step we are allowed to replace $Y$ which from my understanding consists of a known and deterministic vector (actual known labels) with $X\theta + \epsilon = \hat{Y}$ (label prediction of known data points). $Y=\hat{Y}$ iff the linear model is assumed to be exact all the time. But the residual sum of squares goes very rarely (or never) to 0, consequently it should never be the case.
Could someone enlighten the situation? Because Im a bit confused right now. Thanks!