Is a moving average model fitted to white noise? I don't understand the following definition of a moving average model from Hyndman 2021, Forecasting: principles and practice

A moving average model uses past forecast errors in a regression-like
model, $$ y_t = c + \varepsilon_t + \theta_1\varepsilon_{t-1} + \theta_2\varepsilon_{t-2} + \dots + \theta_q\varepsilon_{t-q}   $$
where $\varepsilon_t$ is white noise.

Wikipedia makes a similar statement

... and the $\varepsilon _{t},\varepsilon _{t-1},...,\varepsilon _{t-q}$ are white noise error terms.

But if $\varepsilon_t$ is white noise, then how can anything be predicted from it? What should come out from a linear regression model, that one tries to fit to white noise?
Brownlee 2017 on the other hand writes

The residual errors from forecasts on a time series provide another
source of information that we can model. Residual errors themselves
form a time series that can have temporal structure. A simple
autoregression model of this structure can be used to predict the
forecast error, which in turn can be used to correct forecasts.

That I can understand. But it implies, that the errors are not white noise. Otherwise there would be no structure to model. So is the definition of Hyndman and Wikipedia wrong or misleading? Is a moving average model fitted to white noise or must some structure exist in the errors so fitting the moving average model makes sense?
 A: Short answer: Yes, the model is in your described first described setup to white noise and the confusion you are facing probably originates from mixed notation across your sources.
An easy way to gain some intuition is to "split" the observed error from your average component $c$ into two different types: A (theoretically, as the occurrence lays in the past) explainable error component $e_t$ and non-explainable, white noise component $\varepsilon_t$. Note that we can rewrite your first equation as:
\begin{equation}
y_t=c+e_t+\varepsilon_t,
\end{equation}
with $e_t = \theta_1\varepsilon_{t-1}+\theta_2\varepsilon_{t-2}+...+\theta_k\varepsilon_{t-k}$. In order to estimate these kind of models in an appropriate way, one either relies on yule-walker equations or (under some circumstances feasible) reformulations into an $AR(\infty)$-representation. You may find details here. The idea of simply fitting an $AR(k)$-model to the residuals looks conceptually similar, the key difference is however that the difference equation for $e_t$ is formulated based on true, unobservable errors and not based on estimates for this error.
