# source("shrinkiv.r")
# This is an R file. It creates the figures reported in "A Stein-Like 2SLS Estimator" by Bruce E. Hansen

library("MASS")

rep  = 100000
grid = 40

ns  = c(100,800)
r2s = c(.01,.1,.4)
ms  = c(1,2,3,4)

nn  = length(ns)
r2n = length(r2s)
mn  = length(ms)

mse=array(0,dim=c(nn,r2n,mn,grid,3))

for (ni in 1:nn){
for (ri in 1:r2n){
for (mi in 1:mn){

n  = ns[ni]
m  = ms[mi]
r2 = r2s[ri]
rhos = (0:(grid-1))/grid
c = sqrt(r2/(1-r2))

for (rhoi in 1:grid) {
rho = rhos[rhoi]
tau = (m-2)*(m>2) + (m==2) + (m==1)*(1/4)

store = matrix(0,rep,3)

for (i in 1:rep) {

z = matrix(rnorm(n*m),n,m)
v = matrix(rnorm(n*m),n,m)
y = matrix(rnorm(n),n,1)*sqrt(1-rho^2) + v%*%matrix(1,m,1)*(rho/sqrt(m))
x = z*c + v
z1= matrix(1,n,1)
x = cbind(x,z1)
z = cbind(z,z1)

xxi = solve(crossprod(x))
bols = xxi%*%crossprod(x,y)
eols = y - x%*%bols
sigols = crossprod(eols)/(n-m-1)
bols = bols[1:m,1]
v1 = xxi[1:m,1:m]*matrix(sigols,m,m)

zz = solve(crossprod(z))
zx = crossprod(z,x)
zzi = ginv(t(zx)%*%zz%*%zx)
b2sls = zzi%*%(t(zx)%*%zz%*%crossprod(z,y))
e2sls = y - x%*%b2sls
sig2sls = crossprod(e2sls)/(n-m-1)
b2sls = b2sls[1:m,1]
v2 = zzi[1:m,1:m]*matrix(sig2sls,m,m)

h = t(b2sls-bols)%*%ginv(v2-v1)%*%(b2sls-bols)
w = (tau/h)*(h>tau) + (h<=tau)
bstein = bols%*%w + b2sls%*%(1-w)

store[i,1]=crossprod(bols)
store[i,2]=crossprod(b2sls)
store[i,3]=crossprod(bstein)

}
mse[ni,ri,mi,rhoi,1] = median(store[,1])
mse[ni,ri,mi,rhoi,2] = median(store[,2])
mse[ni,ri,mi,rhoi,3] = median(store[,3])

}

windows()
ltype=c(1,5,2)
lthick=c(1,1,2)
lcolor=c(1,2,3)
mser = mse[ni,ri,mi,,]
mmax =max(mse[ni,ri,mi,,2])*1.5
xlabel = "Degree of Endogeneity (rho)"
ylabel = "Median Squared Error"
mlegend = c("OLS","2SLS","Stein")
figname <- paste("figure.",n,".",ri,".",mi,".eps",sep="")

matplot(rhos,mser,type="l",xlab=xlabel,ylab=ylabel,lty=ltype,lwd=lthick,col=lcolor,ylim=c(0,mmax))
legend("bottomright",legend=mlegend,lwd=lthick,lty=ltype,col=lcolor,bty="n")
savePlot(file=figname,type="eps",dev.cur())

}
}
}
save(mse,file="mse")