# Issue about confidence interval on OLS intercept

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
}
}

• 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. – whuber Jan 8 at 0:08
• 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 – omega Jan 8 at 0:42
• 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. – whuber Jan 8 at 2:22
• thank you for the help – omega Jan 8 at 6:09