Confidence intervals around functions of estimated parameters Suppose I estimate the following model,
$$
 y = \hat{\beta_0} + \hat{\beta_1} X_1 + \hat{\beta_2} X_2\ ,
$$
and am additionally interested in the quantity
$$
 \gamma = \sqrt{\beta_0^2 - 4\beta_1\beta_2}\ .
$$
An estimate of the above, $\hat{\gamma}$, is clearly given by inserting the estimated $\hat\beta_i$'s in the above formula.
But how does one go about computing the confidence interval around such a "compound" quantity that is not directly estimated?
My attempt.
My "empirical" approach is as follows. I generate a multivariate normal distribution of size $n$ using the estimated means $\hat\beta_i$'s and the covariance matrix of the fit. I consequently have a distribution of size $n$ for the $\gamma$. The mean of this distribution $\bar\gamma$ is naturally then close to the $\hat\gamma$ computed above, and I now additionally have a standard deviation of this distribution $\sigma$.
Using $\sigma$ to generate a confidence interval.
If I use the standard definition $\bar\gamma \pm t(\alpha, n-1) \sigma / \sqrt{n}$, where $t$ is the PPF of the $t$-distribution, I run into the following quandary: since I pick $n$ myself, I can in fact make the standard error as small as I like. This is not what I want: I am looking for a "natural" confidence interval which depends solely on the uncertainty in my estimates of the $\beta$'s. Therefore, my idea is the following: $\sigma$ stays positive, and seems to "converge" to its true value for large $n$: hence, 95% of the values in my $\gamma$ distribution lie within the interval $\bar\gamma \pm 2\sigma$.
Is this approach valid?

Usually we take normality assumption for linear regression models. That is, $y_i\sim N(\beta^Tx_i,\sigma^2)$. From this assumption we derive the asymptotic distribution of $\hat{\beta}$, which is also normal: $\hat{\beta}\overset{.}{\sim}N(\beta,\sigma^2(X^TX)^{-1})$. Let us denote the covariance matrix of $\hat{\beta}$ as $Var(\hat{\beta})=V=\sigma^2(X^TX)^{-1}$.
Next, we take $\gamma=g(\beta)=\sqrt{\beta_0^2-4\beta_1\beta_2}$ so $\hat{\gamma}=g(\hat{\beta})$. The function $g$ is differentiable and let us assume it is non-zero. Equipped with these, we can apply the delta method to get the following estimator:
$$\hat{\gamma}\sim N\left(\gamma,\nabla g(\hat{\beta})^T\cdot V\cdot \nabla g(\hat{\beta})\right)$$
Where $\nabla g(\hat{\beta})$ is the gradient vector of $g$ :
$$\nabla g(\beta)=\begin{pmatrix} \frac{\beta_0}{\sqrt{\beta_0^2-4\beta_1\beta_2}} \\ \frac{-2\beta_2}{\sqrt{\beta_0^2-4\beta_1\beta_2}} \\ \frac{-2\beta_1}{\sqrt{\beta_0^2-4\beta_1\beta_2}}
\end{pmatrix}$$
Now, you can easily write the value of $\nabla g$ at $\hat{\beta}$ and then substitute $\hat{\beta}=(X^TX)^{-1}X^Ty$ and obtain a nicer form for the variance of $\hat{\gamma}$ (this does require some linear algebra work). Finally, you can write down your CI.
