Data set

n = 30
x.mat = matrix(rnorm(n,mean=0,sd=1),nrow=n,ncol=1)
err.vec = rnorm(n,mean=0,sd=(1+x.mat^2/2))
y.vec = 1 + x.mat + err.vec
data = data.frame(y.vec,x.mat)
o.fit = lm(y.vec~x.mat)
sd=(1+x.mat^2/2)
g.fit = lm(y.vec~x.mat, data = data, weights = 1/sd)

Simulation

ols.vec = NULL
gls.vec = NULL

for(i in 1:500 ) {
 n = 30
 x.mat = matrix(rnorm(n,mean=0,sd=1),nrow=n,ncol=1)
 err.vec = rnorm(n,mean=0,sd=(1+x.mat^2/2))
 y.vec = 1 + x.mat + err.vec
 ols.vec[i] = lm(y.vec~x.mat)$coefficient[2]
 gls.vec[i] = lm(y.vec~x.mat, weights = 1/sd)$coefficient[2]
}

mean

mean(ols.vec)
## [1] 0.9914118
mean(gls.vec)
## [1] 0.9981318

variance

var(ols.vec)
## [1] 0.2225147
var(gls.vec)
## [1] 0.2422072

histogram

hist(ols.vec)

hist(gls.vec)