Autocorrelation in residuals of a regression model with ARIMA errors (example in Rob Hyndman's book) - Part 2 Part 1 is here
What does the forecaster do when there's correlation in the residuals of an ARIMA model that's used to model the errors from a regression model? Does this mean the forecasting approach - regression model with ARIMA errors - is not suitable and cannot/shouldn't be used or what steps can be taken to produce good points forecasts if the forecaster must necessarily use the regression model with ARIMA errors approach?
 A: There is a difference between forecasting into the future (predicting $y_{t+1}$ based on $y_t$) and contemporaneous prediction (predicting $y_t$ based on $x_t$).
As discussed in the linked question, forecasting into the future necessarily involved lagged dependent variables in the regression. In this case, serial correlation in the residuals indicates serial correlation in error term. This would be problematic.
For contemporaneous prediction in a time series regression with no lagged dependent variables, valid predictions and prediction intervals can be computed under very general serial correlation and heteroskedascity conditions for the error term, under the key exogeneity assumption.
Empirically, as long as the regressors are exogenous, estimates are consistent, and gives consistent predicted value. Prediction errors can be computed by applying HAC procedure to the residuals.
Take the simplest example,
$$
y_t = \beta x_t + \epsilon_t.
$$
As long as exogeneity holds, i.e. $E[x_t \epsilon_t] = 0$, or even under the weaker condition that it holds "in the long run"
$$
\lim_{T \rightarrow \infty}\frac{1}{T} \sum_{t=1}^T E[x_t \epsilon_t] = 0
$$
the regression estimate $\hat{\beta}$ is consistent, and $\hat{\beta} x_t$ is a consistent predictor of $y_t$. In the context of prediction, exogeneity is customarily strengthened to $E[\epsilon_t|x_t] = 0$. So the best prediction is   $E[y_t|x_t] = \beta x_t$.
The population prediction error would just be the long-run variance of $\epsilon_t$.
The corresponding sample quantity can be computed by applying a HAC computation to the residuals.
(One can plug in/assume/forecast future values of $x_{T+2}$ and predict $y_{T+2}$, but this is empirical practice.)
Further Comments
Both no lagged dependent variable and exogeneity are assumptions.
They cannot be verified statistically and their validity rests on empirical justifications.

*

*Exogeneity $E[x_t \epsilon_t]$ is by definition a statement about what is not observed ${\epsilon_t}$, therefore cannot be tested statistically. You have to justify empirically that everything you do not observe is uncorrelated with the regressor $x_t$.
Serial correlation and heteroskedasticity in the residuals are not a problem only if exogeneity holds.


*For example, if the $y_t$ depends on its lagged value $y_{t-1}$ but $y_{t-1}$ is omitted from the regression, then exogeneity would not hold.
In this case, there would be serial correlation and heteroskedasticity in the residuals
Therefore, just like exogeneity, having no lagged dependent variables in the model is a choice. It implies that you have made assumption that $y_t$ does not depend on its lagged value, which then allows you to conclude non-whiteness of residuals is OK.
For example, suppose the true model
is
$$
y_t = \phi y_{t-1} + \beta x_t + \epsilon_t,
$$
and you fit the model
$$
y_t = \beta x_t + \epsilon_t.
$$
If you mistakenly assumed exogeneity, you will conclude that the serial correlation you observe in the residuals is not due to omitted lagged dependent variable (LDV), and mistakenly conclude $\hat{\beta}$, and corresponding predicted value, is consistent.
Data series from these models are observationally indistinguishable. Is the serial correlation in the residuals due to autoregression of the dependent variable or due to serial correlation in an exogenous error term? There's no statistical test that distinguishes the two cases.
Imposing a parametric ARMA structure on $(\epsilon_t)$ would not fix this problem.
(In the quoted example involving electricity demand and temperature, the model could well be correctly specified with no lagged dependent variables.  I don't know nearly enough about the electricity market to say either way.)
Caveat
All this is relevant only if you care about the best prediction $E[y_t|x_t] = \beta x_t$. If you're only interested in the best linear prediction, go ahead, run the regression and use $\hat{\beta} x_t$.
In this case, the bias in $\hat{\beta}$ is not a concern, since you don't really care about estimating the "true model".
The OLS estimate, by construction, consistently estimates the linear correlation between $x$ and $y$.
In situations where you believe lagged variables play a role, they should certainly be included. Serial correlation in the residuals may suggest relevant lagged variables are being omitted, which leads to loss of  predictive power.

Fine print in response to comments:

*

*Best prediction of $y_t$ based on $x_t$ is $E[y_t | x_t]$.
It is "the function $f(x_t)$ of $x_t$ that minimizes $E[( f(x_t)- x_t )^2]$, informally.


*Best linear prediction of $y_t$ based on $x_t$ is $\frac{Cov(x_t, y_t)}{Var(x_t)} x_t$. It is "the linear function $f(x_t)$ of $x_t$ that minimizes $E[( f(x_t)- x_t )^2]$. By construction, regression estimate $\hat{\beta}$ will "always"  consistently estimate $\frac{Cov(x_t, y_t)}{Var(x_t)}$.

A: In practice of forecasting there's very little that is absolute. This is one such case where there is not prescribed course of actions. Presumably you started with a time series regression model $y_t=X_t\beta+\varepsilon_t$ where $\varepsilon_t\sim\mathcal N(0,\sigma^2)$.
Once you looked at residuals $\hat\varepsilon_t$ and noticed that they're autocorelated, you decided to improve the model and apply regARIMA model:
$$y_t=X_t\beta+\varepsilon_t$$
where $\varepsilon_t=\phi_1\varepsilon_{t-1}+u_t$ with $u_t\sim\mathcal N(0,\sigma^2_u)$
Then you find that residuals $\hat u_t$ are autocorrelated. Now what? You could try to fit higher order ARIMA(p,d,q) instead if the first attempt with AR(1). In fact if you pick high enough orders of P,D,Q, I bet that at some point residuals $\hat u_t$ will start looking like a white noise. Should you do this? Maybe, maybe not. It's up to you.
I prefer parsimonious models, and dislike high order models, especially when it comes to differencing D. You also need to be careful with autocorrelation measures since they're sensitive to outliers. For instance, you may have two big events 6 months apart and if the dataset is not large, they'll appear like 6 month frequency seasonality.
