rm(list = ls())

# generate heteroskedastic data

library(MASS)

n <- 200

# set parameter values

b <- as.matrix(c(1,1))
a <- as.matrix(c(.5,.5))

ols <- function(Y, X){
       n <- length(Y)
       X <- cbind(rep(1,n), X)
       k <- ncol(X)
       XXinv <- solve(t(X)%*%X)
       beta <- XXinv%*%t(X)%*%Y
       e <- Y - X%*%beta
       sig2 <- as.numeric(t(e)%*%e/(n-k))
       se <- sqrt(diag((sig2*XXinv)))

       
       list(coef=beta, sig2=sig2, se=se, e=e)
       }

gls <- function(Y, X){
       n <- length(Y)
       Y <- ????
       X <- ????
       k <- ncol(X)
       XXinv <- solve(t(X)%*%X)
       beta <- XXinv%*%t(X)%*%Y
       e <- Y - X%*%beta
       sig2 <- t(e)%*%e/(n-k)
       se <- sqrt(diag(XXinv)*sig2)
       
       list(coef=beta, sig2=sig2, se=se, e=e)
       }

# number of simulations
 sim <- 1000

 bsim <- rep(0,sim)

 bglssim <- rep(0,sim)

 for(j in 1:sim){

  x  <- rnorm( n )
  sigmai <- a[2]*x^2	
  e <- mvrnorm(1, rep(0,n),diag(sigmai))
  xb <- b[1] + b[2]*x 
  y <- xb + e

  result <- ols(y,x)
  bsim[j]<-result$coef[2]

  resultgls <- gls(y,x)
  bglssim[j]<-resultgls$coef[2]

 }

cat("\n")
cat("Bias in OLS Beta:", mean(bsim)-b[2], "\n")
cat("Bias in GLS Beta:", mean(bglssim)-b[2], "\n")
cat("Standard Deviation of OLS beta:", sd(bsim), "\n")
cat("Standard Deviation of GLS beta:", sd(bglssim), "\n")