If you use the log returns, you're essentially making the assumption that there is no conditional variation in the mean. In some circumstances you may want to explicitly model both, but other times it may be sufficient to assume a constant mean and focus on the conditional variance. Depends on what you're trying to do.
In addition, if you fit a GARCH model with raw log returns, then you're also implicitly assuming the mean is zero. Centering the data may be important if the mean is large (i.e. especially in lower frequency data).
But is the one step ahead predictor not already defined as the value $\hat \sigma$ of the volatility that minimizes the MSE?
If you estimate the GARCH model using maximum likelihood then the fitted values $\hat\sigma_t$ are the likelihood-maximizing values (subject to the GARCH(1,1) functional form) which need not coincide with MSE-minimizing values. That depends on the distribution assumed for the likelihood calculation.
Also, when fitting the model on a data sample indexed from $1$ to $T$, the fitted value $\hat\sigma^2_t$ for $t<T$ utilizes the information not only from $1\leqslant\tau<t$ but also from $t\leqslant\tau\leqslant T$. Meanwhile, to fairly evaluate the forecast performance, you should not allow the estimate to be calculated using future data; instead, you need to check the forecast accuracy out of sample.
So why compute this measure if it going to be the minimum across models anyway?
Because you may want to know what the actual values of MSE is. Knowing that the MSE is minimal does not tell you what its value is. (See also the answer to the previous question.)
Moreover do I understand correctly that the practical way to compute the $MSE$...
That will be the mean squared forecast error. If you want in-sample MSE, just use the fitted values from the model estimated on the whole sample. The former should give an unbiased estimate of model performance, while the latter will be too optimistic in that respect.
Also I understand from this paper (Bollerslev 1998) that utilizing the squared daily return to approximate the realized volatility leads to noise.
Very important point indeed! The response variable of the GARCH model is measured with noise when squared errors are used as proxies; this noise may be quite substantial. When trying to assess model fit, the measurement error associated with the dependent variable may cause quite some trouble.
For example a better estimate of realized daily volatility would be the sum of 30 minutes squared returns of that day.
On the first thought, that could be a valid option. If that was proposed in the Andersen and Bollerslev (1998) paper, then it must be fine.
So to get a "better" MSE I could substitute every σi with the sum of 30 minutes squared returns of that day instead of simply the daily squared return? Is this how you use realized volatility to evaluate the goodness of your forecasts?
That makes sense to me.
I will end this rambling by asking for a good reference in evaluating the accuracy of the forecasts using realized volatility...
Regarding references, I think the Andersen and Bollerslev (1998) paper is quite relevant and complete. Also see Patton & Sheppard "Evaluating volatility and correlation forecasts" (2009) (free version here) and other works by Patton.
Best Answer
Lets start with the set up of a simple GARCH(1,1) model $$ y_t=\sigma_t\varepsilon_t,\;\;\varepsilon_t \stackrel{iid}{\sim}N(0,1) $$ $$ Var(y_t | y_{t-1},\sigma_{t-1}) = \sigma^2_t = \alpha_0 + \alpha_1 y^2_{t-1}+ \beta_1\sigma^2_{t-1} $$ So we know $$ y_t|y_{t-1},\sigma_{t-1} \sim N(0,Var(y_t | y_{t-1})) \quad(1) $$
$(1)$ gives you the 1-step ahead predictive/forecast density and the best predictive density in the GARCH(1,1). Hopefully that makes sense.