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?
 A: 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.
A: It seems like you are estimating the discriminant of a quadratic function, ie. your function is
$$y = \hat{\beta_0} + \hat{\beta_1} X_1 + \hat{\beta_2} X_2 =
 \hat{\beta_0} + \hat{\beta_1} X + \hat{\beta_2} X^2$$
We can rewrite this in terms of the roots
$$y = \hat{\beta_2} (X-r_1)(X-r_2)$$
And now we have $$\gamma = \beta_2 (r_2-r_1)$$
If you estimate $\beta_2$, $r_1$, and $r_2$ with a non-linear estimation method and obtain an estimate for the approximate multivariate normal distribution of the estimates (based on an estimate of the Fisher information matrix), then you can describe an estimate distribution for $\gamma$ as a product of two correlated non central normal distributions.
A: There's one big problem with your proposed method, although it can be fairly easily remedied. In particular, you compute the estimate of variance of $\gamma$ as $\frac{V(\gamma_i)}{n}$ (where $V(\gamma_i)$ is variance of your samples generated by first sampling from the multi variance normal and then pushing through your function), when the standard error should really just be $V(\gamma_i)$.
But furthermore, I would suggest you use the simulated percentiles to create the CI rather than using the mean + se formula.
And one final note, we note that the sqrt function is not defined for negative values. This means if a random draw gives you $4\beta_1 \beta_2 > \beta_0^2$, your estimate is not real valued. Assuming that this is not what you want to do, I would say you should consider simply dropping the simulated values when this is true. That may seem adhoc, but it perfectly aligns with defining a restricted parameter space where $4\beta_1 \beta_2 \le \beta_0^2$.
