Based on $N$ realizations of two random variables $X \sim N(0,\sigma_X^2)$ and $Y \sim N(0, \sigma_Y^2)$ with correlation $\rho$, I conduct a simple linear regression $Y = \beta_0 + X\beta_1 + \epsilon$.
Using the known distributions of $X$ and $Y$, their correlation $\rho$ and the number of realizations $N$, I am trying to find the sampling distribution of the regression coefficients.
Estimates of the coefficients can be caluclated as:
$\hat \beta_1=r\frac{s_Y}{s_X}$ where $r$ is an estimate of $\rho$ and $s_X, s_Y$ are estimates of $\sigma_X, \sigma_Y$
$\hat \beta_0=\overline{Y} - \hat \beta_1 \overline{X}$ where $\overline{X}, \overline{Y}$ are estimates of the means of $X$ and $Y$
Intuitively, I thought the regression coefficients should follow a normal distribution. However, after some additional thoughts, I am not sure of anymore. Plotting the distribution of $\beta_1$ from 10,000 experimental runs with $N=12, \sigma_X^2=1, \sigma_Y^2=2, \rho =0.6$ (black) vs. fitted normal distribution, a deviation from normal distribution can be seen. This is also the case more iterations (I tested up to 100,000) for which the normal distribution should be quite stable. For more detail on the experiment, see below.
Finding the specific distribution of $\beta_0, \beta_1$ - if there is one - is however only the second step. In the first step, I am trying to find the expected value and variance of the regression coefficients.
It is clear that $E[\hat \beta_0]= 0$ and $E[\hat \beta_1] = \rho \frac{\sigma_Y}{\sigma_X}$.
However, I have trouble deriving the variance of both parameters. I have found some definitions for the coefficient variance in the literature. $Var[\hat \beta_0] = \frac{s_\epsilon^2\sum x_i^2}{N \sum (x_i - \overline{x})^2}$ and $Var[\hat \beta_1] = \frac{s_\epsilon^2}{\sum (x_i - \overline{x})^2} = \frac{s_\epsilon^2}{N \sigma_X^2}$.
However, as stated correctly by whuber, these definitions are conditioned on $X$. A calculation is therefore only possible for specific relizations and the definitions are therefore not applicable to my use case. In the literature, I found that the restriction on fixed $X$ is going back to Fisher (e.g. his work 'Asymptotic distribution of the reduced-rank regression estimator under general conditions' from 1922). However, I did not find a consideration of the case with random $X$ in more recent literature. I have only found refrences to the (non central) Wishart distribution. I'm not sure whether and how the Wishart distribution can be used in my use case, though.
Overall, I am completly stuck on how to derive the variances in the described case, $Var[\hat \beta_0], Var[\hat \beta_0]$ with $X, Y$ being random variables.
For cases with conditioning on $X$, I have found various answers on how to derive the variance of the regression coefficients, for instance https://stats.stackexchange.com/a/89155/48067. However, since not only $\epsilon_i$ is a random variable, the approach described in the answer is not easily transferable to my problem.
Consdering that $\hat \beta = r \frac{cov(X, Y)}{s_X^2}= r \frac{s_Y}{s_X}$, I assume that the distributions of $r, s_Y, s_X$ may for some reason be problematic for the calculation of the variance. I have found that $s_X^2, s_Y^2$ follow a gamma distribution. $s_X, s_Y$ should therefore follow a generalized gamma distribution, which have a well defined variance. However, I am not sure about the quotient of two generalized gamma distributed random variables. $r$ has a skewed distribution on which the Fisher transformation can be used - I am however not sure whether this helps to calculate the variance. Lastly, the sample covariance has a strange distribution, approximately like a shifted gamma distribution.
Overall, I have already derived some information about the distributions of characteristics of the samples. Yet, I have not found a way to use these information to derive the variance of the sample coefficients in my use case.
Can someone point me in the right direction?
--
Some additional detail about my "simulations", as requested:
- Samples of size N (e.g. N=12) are drawn from $X \sim N(0, \sigma_X^2)$ and $Y \sim N(0,\sigma_Y^2)$ (e.g. $\sigma_X^2 = 1$, $\sigma_Y^2 = 2$) so that the samples have a specific correlation $\rho$ (e.g. $\rho = 0.6$). This is achieved by using two uncorrelated random variables $A, B \sim N(0,1)$, and constructing a new (correlated) one as $C = \rho^2 A + \sqrt{1- \rho^2} B$ and correcting the variance of $A$ and $C$ in order to have $X$ and $Y$
- A regression of $Y$ on $X$ is executed leading to coefficients $b_0$ and $b_1$
- Steps one and two are repeated (for instance 10,000 times)
- The distribution resulting of $b_0$ and $b_1$ resulting from the runs are analyzed
The distributions resulting from the last step are the distributions for which I want to derive exact distribution parameters in order to analyze the uncertainty in the regression parameters resulting from the sampling of the random variables.
