set.seed(234235)

const  = FALSE

n        = 1000
sig2    = 0.5
betas = c(0.5,1.0,2.0,0.5)
gammas = c(1,2,0.0,0.0)
p = length(betas)
k = p - 1
x = matrix(runif(n*k,0,2),n,k)
X = cbind(1,x)

if (const){
  y = X%*%betas+rnorm(n,0,sqrt(sig2))
}else{
  h = exp(X%*%gammas)  
  y = X%*%betas+rnorm(n,0,sqrt(sig2*h))
}

# OLS regression
reg = lm(y~x)
summary(reg)
se.ols = sqrt(diag(summary(reg)$cov.unscaled))*summary(reg)$sigma
coef.ols = reg$coef

# Computing hetereskedasticity-robust standard errors
se.hr = rep(0,ncol(X))
for (i in 1:p){
  reg1 = lm(X[,i]~X[,-i]-1)
  se.hr[i] = sqrt(sum((reg1$res^2)*(reg$res^2))/(sum(reg1$res^2))^2)
}
round(cbind(betas,coef.ols,se.ols,se.hr),3)

# Hetereskedasticity-robust F test
esq     = reg$res^2
R2.e   = summary(lm(esq~x))$r.sq
Ftest   = R2.e/(1-R2.e)*(n-p)/k
pvalue = 1-pf(Ftest,k,n-p)
pvalue

# Generalized least squares
lesq = log(esq)
g    = lm(lesq~x)$fit
hhat = exp(g)
reg.gls = lm(y~x,weights=1/hhat)
se.gls = sqrt(diag(summary(reg.gls)$cov.unscaled))*summary(reg.gls)$sigma
coef.gls = reg.gls$coef


round(cbind(betas,coef.ols,se.ols,se.hr,coef.gls,se.gls),3)