I begin with introductin some notation. The matrix $X\in\mathbb R^{n\times p}$ contains $n$ observations of $p$ features and the vector $y\in\mathbb R^n$ contains $n$ observations of the response variable. $\bar x\in\mathbb R^{1\times p}$ contains the means of the columns of $X$ and $\bar y$ is the mean of $y$. Denote the centred versions of $y$ and $X$ by $$ y^c =y-\mathbf 1\bar y \quad\text{and}\quad X^c=X-\mathbf 1\bar x, $$ where $\mathbf 1\in\mathbb R^{n\times 1}$ is a vector of ones.
Consider a linear model $$ y^c=X^c\beta+\varepsilon, $$ where $\beta\in\mathbb R^p$ is a vector of parameters and $\varepsilon\in\mathbb R^n$ is vector of errors such that $\operatorname E\varepsilon=0$ and $\operatorname E[\varepsilon\varepsilon']=\sigma^2I$ with $\sigma^2>0$. $\beta$ is estimated using least squares, i.e. $\hat\beta=((X^c)'X^c)^{-1}(X^c)'y^c$.
Suppose that $z\in\mathbb R^{1\times p}$ contains new observations of the $p$ features. Then the predicted value of $y$ is given by $$ \hat y =\bar y+(z-\bar x)\hat\beta. $$ What is the expression of the variance of the predicted value $\hat y$?
I saw the following expression with no explanation on how it is derived $$ \operatorname{Var}(\hat y) =\sigma^2\Bigl\{\frac1n+(z-\bar x)((X^c)'X^c)^{-1}(z-\bar x)'\Bigr\}. $$ This is the right expression if $\bar y$ and $(z-\bar x)\hat\beta$ are uncorrelated. Indeed, $\operatorname{Var}\bar y=\sigma^2/n$ and $$ \operatorname{Var}((z-\bar x)\hat\beta) =\sigma^2(z-\bar x)((X^c)'X^c)^{-1}(z-\bar x)' $$ using the fact that $\operatorname{Var}(\hat\beta)=\sigma^2((X^c)'X^c)^{-1}$. Are $\bar y$ and $(z-\bar x)\hat\beta$ uncorrelated?
Any help is much appreciated!