Solution (source code)

= Solution

The future <mean> response is $E_*=\alpha+\beta x_*+\gamma z_*$. By invariance of <maximum likelihood estimation>, its estimator is the fitted <mean>
$$
\boxed{\widehat E_*=\widehat\alpha+x_*\widehat\beta+z_*\widehat\gamma.}
$$
Its sampling <variance> is the quadratic form of the coefficient <covariance matrix> with $c=(1,x_*,z_*)^T$:
$$
\boxed{\operatorname{Var}(\widehat E_*)=c^T\operatorname{Var}(\widehat\theta)c
=\sigma^2\left(\frac1n+\frac{x_*^2}{S_{xx}}+\frac{z_*^2}{S_{zz}}\right).}
$$
This estimator is unbiased because each fitted coefficient is unbiased. The future error $\epsilon_*$ does not enter the estimator of the <mean>, so its independent <variance> is not part of this expression.