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)