## Demonstrate the Central Limit Theorem for beta.hat.OLS:
## For samples "large enough", beta.hat.OLS ~ N(beta, sigma2.beta)
## *regardless* of how u is distributed.

## the DGP is y = alpha + beta*x + u, set e.g. alpha=2, beta=3

## To demonstrate (not "prove"!) this result, we simulate a large
## number of draws from "our" DGP, using non-Normal errors

## set the number of replications
reps <- 100

## set the sample size to a "large" number
n <- 5000

## set "true" parameter values
alpha <- 2
beta <- 3

## initialize betas vector
beta.hats <- rep(NA, reps)

## draw x (only once: we maintain the hyp. "x non stochastic")
x <- rnorm(n)

## now make 100 samples from "our" DGP

for(i in 1:reps) {
  ## draw error for this sample, from the uniform distribution
  u <- runif(n, min=-0.3, max=0.3)
  ## make y according to DGP
  y <- alpha + beta*x + u
  ## estimate linear model over this i-th sample and get beta.hat.OLS
  beta.hats[i] <- coef(lm(y ~ x))[2]
}

## the "typical" error we generated looks like this:
u0 <- runif(n, min=-0.3, max=0.3)
hist(u0, freq=FALSE)

## compare with a Normal:
## graphical comparison of quantiles (theoretical vs. empirical)
qqnorm(u0)
## formal Jarque-Bera test
library(tseries)
jarque.bera.test(u0)

## the vector of 100 estimates:
hist(beta.hats)

## check normality:
qqnorm(beta.hats)
jarque.bera.test(beta.hats)


## what about generating errors from another distribution? Try a Chi2

## make 100 samples as above, errors are sampled from a (recentered)
## Chi2(5) this time:
## (notice: betas etc. are overwritten)

for(i in 1:reps) {
  ## draw error for this sample, from the uniform distribution
  u <- rchisq(n, df=5)
  ## recenter u:
  u <- u - mean(u)
  ## make y according to DGP
  y <- alpha + beta*x + u
  ## estimate linear model over this i-th sample and get beta.hat.OLS
  beta.hats[i] <- coef(lm(y ~ x))[2]
}

## the "typical" error we generated this time is very asymmetric:
u0 <- rchisq(n, df=5)
u0 <- u0 - mean(u0)
hist(u0, freq=FALSE)

## compare with a Normal:
## graphical comparison of quantiles (theoretical vs. empirical)
qqnorm(u0)
## formal Jarque-Bera test
jarque.bera.test(u0)

## the vector of 100 estimates:
hist(beta.hats)

## check normality:
qqnorm(beta.hats)
jarque.bera.test(beta.hats)