####################################################################
##Execute following on dos command prompt (omit the "##"):
##"S:\\R\\rw2001\\bin\\Rterm.exe" -q --no-restore --no-save < bias.r
####################################################################

options(echo = FALSE)

png(file="bias.png", width=680, height=450, pointsize=12, bg="white", res=200)

#----------------------------------------------------------------
# Inputs:
n <- 5;
m <- 3;
w <- c(-0.07343, 0.29371, 0.55944, 0.29371, -0.07343);
a <- 0;
var <- 0.1*seq(1,100,10);
beta <- 0.01*seq(1:300);
#----------------------------------------------------------------

delta <- (1/(m*a+0.5*beta*m*(m+1)))*(a*sum(w[(m+1):n]) +
          beta*(sum(seq(1:n)*w[1:n]) - sum(seq(1:m)*w[1:m])));

#For a given beta[[i]]; asymmetric weights: u1[[i]], u2[[i]], u3[[i]]:
u1 <- w[[1]] + delta;
u2 <- w[[2]] + delta;
u3 <- w[[3]] + delta;

for(j in 1:length(var)) {

  R <- 0;
  VarR <- 0;
  StdErr <- 0;
  for(i in 1:length(beta)) {

    u <- c(u1[[i]], u2[[i]], u3[[i]]);

    R[[i]] = sum(u*(a + beta[[i]]*seq(1:m))) - sum(w*(a + beta[[i]]*seq(1:n)));

    VarR[[i]] = (sum(delta[[i]]*(a + beta[[i]]*seq(1:m))))**2 + 
                var[[j]]*(sum((delta[[i]])**2) + sum((w[(m+1):n])**2)) -
                2*(sum(delta[[i]]*(a + beta[[i]]*seq(1:m)))*
                sum(w[(m+1):n]*(a + beta[[i]]*seq((m+1):n)))) +
                (sum(w[(m+1):n]*(a + beta[[i]]*seq((m+1):n))))**2;

    StdErr[[i]] = sqrt(VarR[[i]]); 

  }

  if(j==1) {
    plot(beta,(R+StdErr),type="l",ylim=c(-3,3),
         xlab="Beta",ylab="E(R) +/- StdErr(R)",cex.lab=1.4);
  }

  lines(beta,(R+StdErr));
  lines(beta,(R-StdErr));
}

lines(beta,rep(0,length(beta)),lty=3);

dev.off();
q();