RE=function(data,x,y,est1,est2){
  data=as.name(data); x=x; y=y
  est1=as.name(est1); est2=as.name(est2)
  s2x=var(x); s2y=var(y)
  sx=sd(x); sy=sd(y)
  r=sum(y)/sum(x)
  rho=cor(x,y)
  plot(x,y,main=data)
  if(est1=='srs'){
    var1=s2y
   } else if(est1=='ratio'){
    var1=s2y+r^2*s2x-2*r*rho*sx*sy
   } else if(est1=='regr'){
    var1=s2y*(1-rho^2)
   } else if(est1=='diff'){
    var1=s2y+s2x-2*rho*sx*sy
  }
  if(est2=='srs'){
      var2=s2y
    } else if(est2=='ratio'){
    var2=s2y+r^2*s2x-2*r*rho*sx*sy
    } else if(est2=='regr'){
      var2=s2y*(1-rho^2)
    } else if(est2=='diff'){
    var2=s2y+s2x-2*rho*sx*sy
    }
  vars=c(var1,var2)
  re=round(var2/var1,4)
  cat("","\n","Estimated Relative Efficiency",'\n',"Data: ",data,"\n",'Estimation methods: ',est1,
      'and',est2,'\n','Var1 =',var1,',','Var2 =',var2,'\n','RE(',est1,',',est2,') = V(',est2,')/V(',est1,') =',re,"\n","")
  results=list(data=data,est1=est1,est2=est2,vars=vars,re=re,rho=rho,r=r)  
}
# to use the function with its call:
# RE(data,x,y,est1,est2)
# data: name of dataset, in quotes
# x: vector (data)
# y: vector (data)
# est1, est2: c('srs','ratio','regr','diff')