Say I have some data points randomly distributed about a linear model. I fit a straight line through it, and obtain best fit parameters $m$ and $b$ with an error $\delta_m$ and $\delta_b$ associated to them.
Say I wish to evaluate the fit at a given $x$, call it $x_0$. This obviously gives, as a point estimate, $y_0 = mx_0+b$, but what about the associated error $\delta_{y_0}$?
One might think that we should obtain $\delta_{y_0}$ as such \begin{align} \delta_{y_0} = \Big(\frac{\delta_m}{m}\Big)mx_0 + \delta_b = \delta_mx_0+\delta_b \end{align} But, as I understand it, this has the unfortunate consequence of making it so that $\delta_y$ is not invariant under a translation of the coordinates of $x$. What I mean is, if we had applied a transformation $x \to x+\alpha$, $\alpha \in \mathrm{R}$ to all data points $\underline{\text{before}}$ doing the fit, then, while $\delta_m$ and $\delta_b$ would not have changed, $\delta_{y_0}$, now given by $\delta_{y_0} = \delta_m(x_0+\alpha)+\delta_b$, would have changed, which should not happen.
My intuition is that computing the error via \begin{align} \delta_{y_0} = \delta_m(x_0-\langle x \rangle)+\delta_b \end{align} where $\langle x \rangle$ is the average value of $x$ taken over all data points, is a potential solution to the problem.
What is the correct approach for finding $\delta_{y_0}$? I'm especially interested in the case where we are performing an extrapolation of the line fit, ie, $x_0$ is not in the interval of $x$ values for which data was collected.