##
##	R CODE FOR MONTE CARLO SIMULATION OF TLS&OLS COEFFICIENT BIAS 
##

# CLEAR MEMORY, FORBID SCI NOTATION
rm(list=ls(all=TRUE))				
options(scipen=999)	

# SET RANDOM NUMBER SEED
set.seed(101)

# REGRESSION SAMPLE SIZE; NUMBER OF REPS 
sample = 200
rep = 500	

# TRUE BETA 
beta = 0

# SIGNAL CORRELATION SIGN AND STRENGTH (c IN PAPER)
cs = seq(from=-0.9, to=0.9, by=0.1)

# X NOISE COEFFICIENT (s IN PAPER)
xn = seq(from=0.0, to=1.8, by=0.2)

# TLS FUNCTION
tls = function(w1,w2,y)		{
	v = svd(cbind(w1,w2,y))$v
	b = -v[1,3] / v[3,3]
	return(b)		}


# MATRICES FOR STORING OUTPUT
beta1 = matrix(0,  nrow=rep, ncol=3)		# to hold simulated slope coeffs and corrs
bias.o = matrix(0, nrow=length(cs), ncol=length(xn))		# to hold ols bias estimates 
bias.t = matrix(0, nrow=length(cs), ncol=length(xn))		# to hold tls bias estimates
corrs  = matrix(0, nrow=length(cs), ncol=length(xn))    # to hold corr coeffs

rownames(corrs)=cs
colnames(corrs)=xn
rownames(bias.o)=cs
colnames(bias.o)=xn
rownames(bias.t)=cs
colnames(bias.t)=xn

##
##	SIMULATIONS
##

# Set var(x) = 1
u = 0.5*sqrt(12)				

# Noise scalar on y
noise.y = 1					

# START LOOP ACROSS SIGNAL CORRELATIONS
for (i.corr in 1:length(cs))		{

# START LOOP ACROSS X SIGNAL-NOISE RATIO
for (j.snr in 1:length(xn))		{

# Noise scalar on x
noise.x = xn[j.snr]*noise.y				

# START LOOP ACROSS MC SIMS
for (i in 1:rep)  	{

# GENERATE NOISE TERMS
	u1 = rnorm(sample,0,1)	* noise.x	
	u2 = rnorm(sample,0,1)	* noise.x
	uy = rnorm(sample,0,1)  * noise.y

# Define X, W and Y terms
	x1 = runif(sample,-u,u) 
	x2 = runif(sample,-u,u)/2 + x1*cs[i.corr]	
	y0  = beta*x1 + x2 + uy
	w1 = x1+u1
	w2 = x2+u2
	w1 = w1-mean(w1)
	w2 = w2-mean(w2)
	y0 = y0-mean(y0)
	
# REGRESSIONS
	beta1[i,1] = lm(y0  ~ 0+w1+w2)$coef[1]
	beta1[i,2] = tls(w1,w2,y0)
	beta1[i,3] = cor(w1,w2)
			} # END MC RUN

# TABLE beta1 CONTAINS #rep ROWS OF OLS & TLS ESTIMATES
bias.o[i.corr,j.snr] = mean(beta1[,1])
bias.t[i.corr,j.snr] = mean(beta1[,2])
corrs[i.corr,j.snr]  = mean(beta1[,3])

} # END SIGNAL-NOISE STEPS
} # END CORRELATION STEPS




colors=c("gray0","gray20","gray40","gray60","gray80","red",
         "darkorange","khaki1","lightgreen","dodgerblue")

## FIG 3: GRAPH BIAS AGAINST CORRELATIONS
# Only if true beta=0
if(beta==0) {
#
# SELECT ELEMENTS FOR GRAPHING
bias.1 = matrix(0, nrow=length(cs), ncol=length(xn))		
for (i.corr in 1:length(cs))		{
  for (j.snr in 1:length(xn))		{
    bias.1[i.corr,j.snr]=bias.t[i.corr,j.snr]*(corrs[i.corr,j.snr]< -0.2)
  }
}
bias.1[bias.1 == 0] = NA
yy=5
pdf('bias.v.corr.pdf', width=6, height=6)
plot(bias.1[,1]~corrs[,1], type="l",lty=1,lwd=2,col=colors[1], ylim=c(-2,yy), xlim=c(-.9,0), 
     xlab="Correlation between w's", ylab="Bias", main="TLS Bias")
for (i in 2:10) {
  par(new=TRUE)
  colr=paste("gray",i*10, sep="")
  plot(bias.1[,i]~corrs[,i], type="l",lty=1,lwd=2,col=colors[i],ylim=c(-2,yy), xlim=c(-.9,0), xaxt="n", yaxt="n", xlab="", ylab="")
}
abline(h=0,lty=2)
abline(h=1,lty=2)
legend("topright",legend=xn,col=colors,lty=1,lwd=2,ncol=2,
       bty="o",title="noise factor s", bg="white")
dev.off()
}




## FIG 1 OR 2: GRAPH MAIN COMPARISON
filename = paste("TLS.v.OLS.",beta,".pdf", sep="")
pdf(filename, width=6, height=6)
# TLS
par(mfrow= c(1,2))
par(mar=c(5,5,5,0.5))
yy=5
plot(bias.t[,1]~cs, type="l",lty=1,lwd=2,col=colors[1], ylim=c(-yy,yy), xlab="X correlation term c", ylab="Beta", 
     main="TLS")
for (i in 2:10) {
  par(new=TRUE)
  colr=paste("gray",i*10, sep="")
  plot(bias.t[,i]~cs, type="l",lty=1,lwd=2,col=colors[i], ylim=c(-yy,yy), xaxt="n", yaxt="n", xlab="", ylab="")
}
abline(h=0,lty=2)
abline(h=1,lty=2)
abline(v=0,lty=3)

# OLS
par(mar=c(5,0.5,5,5), yaxt="n")
plot(bias.o[,1]~cs, type="l",lty=1,lwd=2,col=colors[1], ylim=c(-yy,yy), xlab="X correlation term c", ylab="", 
     main="OLS")
for (i in 2:10) {
  par(new=TRUE)
  colr=paste("gray",i*10, sep="")
  plot(bias.o[,i]~cs, type="l",lty=1,lwd=2,col=colors[i], ylim=c(-yy,yy), xaxt="n", yaxt="n", xlab="", ylab="")
}
abline(h=0,lty=2)
abline(h=1,lty=2)
abline(v=0,lty=3)
legend("bottom",legend=xn,col=colors,lty=1,lwd=2,ncol=2,
       bty="o",title="noise factor s", bg="white")
dev.off()


## SAVE RESULTS TO FILES
write.csv(bias.o, file= paste("bias.b", beta,".ols.csv",sep=""))  
write.csv(bias.t, file= paste("bias.b", beta,".tls.csv",sep=""))  
write.csv(corrs,  file= paste("corrs.b",beta,".csv",sep=""))  
round(bias.o,3)
round(bias.t,3)
round(corrs,3)