0
$\begingroup$

Let us assume this simple linear model: $Y|X=\beta_0+\beta_1X+\epsilon $

where $X \sim N(\mu,\sigma^2)$ and $\epsilon \sim N(0,\sigma_{\epsilon}^2)$

Suppose also that $X$ and $\epsilon$ have all the required properties in order to satisfy Gauss Markov theorem in which OLS estimates are Blue.


My purpose is to setup a 95% confidence intervals on regression coefficients and test them by Montecarlo.

To do so I generate $m=10000$ samples of $Y$ and $X$ with varies numerosities $n=10,100,1000,10000$ using $MATLAB$ software (I did this also with R and I got same result).

I used the following notation: $\boldsymbol{\beta}=[\beta_0,\beta_1] \qquad $ $\chi_{n,m}=[1_{n,1}\quad X_{n,m-1}] $

Estimates are $\hat{\boldsymbol{\beta}}=(\chi^t\chi)^{-1}\chi^ty$

I used the variable: $s^2=\frac{\hat{e}^t \hat{e}}{n-2} \quad$ with $\hat{e}$ being residual vector, as estimate of $\sigma_{\epsilon}^2$ in order to find an estimate of the variance-covariance matrix: $\quad\sigma_{\epsilon}^2 (\chi^t \chi)^{-1}$

Is well known that $\hat{\boldsymbol{\beta}} \sim N(\boldsymbol{\beta},\sigma_{\epsilon}^2 (\chi^t \chi)^{-1})$ providing our "nice" assumptions being true. In this framework we can setup a t-student based confidence interval on $\boldsymbol{\beta}$ by using $s^2$: $$ \boldsymbol{\hat{\beta}} \pm t_{1-\alpha/2,n-2}\sqrt{diag(s^2(\chi^t \chi)^{-1})}$$


Basically (on MATLAB) I define two counting variables, one for $\beta_0$ which is increased by 1 if $\beta_0$ falls into its c.i. and same for $\beta_1$ , repeated for $10000$ interaction.

I succesfully get 95% coverage on $\beta_1$ for every choice of $n$.

Things goes unexpectedly bad with $\beta_0$ : I have overall 95% coverage for $n=10$ which drops at ~80% on $n=100$ , ~10% with $n=1000$ and almost 0% at $n=10000$ which seems like the c.i. for the intercept is not appropriate for large $n$.

Any suggestions?

For the sake of completeness I reported a sample code in R using lm and confint function:

# Variable definition
n=100 ; mue=1 ; sigmae=10; m=2;
        mux=5 ; sigmax=3

b0=1.2; b1=-4;

countb0=0 ; countb1=0; 

for (i in 1:10000){
  # Model generation
  eps=rnorm(n,mue,sigmae)
  x1=rnorm(n,mux,sigmax)

    y=b0+b1*x1+eps

  # Linear fit
  fit=lm(y~x1)
  Ci=confint(fit)

  if (Ci[1,1]<=b0 & b0<=Ci[1,2]){
    countb0=countb0+1
  }
  if (Ci[2,1]<=b1 & b1<=Ci[2,2]){
    countb1=countb1+1
  }
}
$\endgroup$
  • 1
    $\begingroup$ Do the calculation correctly? Seriously, there's probably a bug but you don't provide much information to diagnose it. What happens when you do the same simulation with the simpler model that omits $X$, so that $y=\beta_0+\epsilon$? That should be easier to check and diagnose. $\endgroup$ – whuber Jan 8 at 0:08
  • $\begingroup$ I did it, same issue. I get same result even if I use the lm() function in R, extracting confidence interval with confint function; I added R code in main question $\endgroup$ – omega Jan 8 at 0:42
  • 2
    $\begingroup$ You have the wrong value of $b_0:$ it's essential that the mean of $\varepsilon$ be zero, for otherwise that mean must be added to $b_0.$ This would have been simpler to find by following my suggestion to drop $x$ from the model. $\endgroup$ – whuber Jan 8 at 2:22
  • 1
    $\begingroup$ thank you for the help $\endgroup$ – omega Jan 8 at 6:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.