1
$\begingroup$

I am trying to understand the idea of the sampling distribution of an OLS estimator. To do so I am trying to visualise this distribution obtained with simulations (I use R for this).

So I generate the universe of observations with 100,000 randomly distributed observations of my independent variable univ.x. Then I also generate the universe of observations for the dependent variable (univ.y) according to the following (exact) linear relationship: univ.y=1.5+3*univ.x.

> univ.x <- rnorm(100000,0,1)
> univ.y <- 1.5 + 3*univ.x
> summary(lm(univ.y~univ.x))

Call:
lm(formula = univ.y ~ univ.x)

Residuals:
       Min         1Q     Median         3Q        Max 
-2.900e-15 -3.000e-16 -2.000e-16  0.000e+00  8.498e-12 

Coefficients:
            Estimate Std. Error   t value Pr(>|t|)    
(Intercept) 1.50e+00   9.75e-17 1.538e+16   <2e-16 ***
univ.x      3.00e+00   9.73e-17 3.083e+16   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.083e-14 on 99998 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 9.507e+32 on 1 and 99998 DF,  p-value: < 2.2e-16

So far so good. At this point, I want to sample 1,000 observations from the universe and I would expect to find the same linear relationship with some random error. But I don't, in fact I find only the intercept but the slope is 0 and not 3. I have also repeated this 100 times with the same result.

What am I missing here?

Here is the code:

> x <- sample(univ.x, 1000)
> y <- sample(univ.y, 1000)
> summary(lm(y ~ x))

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.9094 -2.0161 -0.0021  1.9563  8.9696 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.505428   0.096626  15.580   <2e-16 ***
x           0.009807   0.094319   0.104    0.917    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.055 on 998 degrees of freedom
Multiple R-squared:  1.083e-05, Adjusted R-squared:  -0.0009912 
F-statistic: 0.01081 on 1 and 998 DF,  p-value: 0.9172

Visualising the sampling distribution of OLS estimates [solved]

Thanks to the comments I solved my problem. I leave the code here for interested users and learners.

univ.x <- rnorm(100000,0,1)
univ.y <- 1.5 + 3*univ.x + rnorm(100000,0,3)
summary(lm(univ.y~univ.x))
plot(univ.x,univ.y)

b <- vector(mode = "numeric", 1000)
for (i in 1:1000){
        samp <- sample(1:100000, 100)
        x <- univ.x[samp]
        y <- univ.y[samp]
        model <- lm(y~x)
        b[i] <- coef(model)[2]
}
hist(b)

Sampling distribution of OLS estimator

$\endgroup$
1
$\begingroup$

You are sampling 1000 of the $x_i$ at random and 1000 of the $y_i$ at random, so you are not necessarily sampling 1000 pairs $(x_i,y_i)$ that fulfill the linear relationship, or that are connected in any way.

Moreover, the data you generated has no randomness in the $y_i$ respect to the $x_i$. Once you knew one observation of $x_i$ you would know the corresponding $y_i$ for sure, right? This is not the real-life setting that OLS or regression is designed for.

The true regression equation you would get from real life data is $$y_i=\beta x_i+\varepsilon_i,$$ where $\varepsilon_i\sim N(0,\sigma2)$. So to obtain randomness you also need to generate the random errors $\varepsilon_i$, all of them drawn from the same distribution $N(0,\sigma2)$.

$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.