# jh 2013.03.30

require(scatterplot3d)

delay=1000000; # adjust this to suit the speed of YOUR computer..

samples=2; # per scenario (column)

par(mfrow=c(samples,2),mar = c(0.01,0.01,0.1,0.1))

B.age=0.4; B.work=0.3;SD.e=2;

age=c(45,45,45,55,55,55,65,65,65);

#order=seq(1,0); # balanced ones on left 
order=seq(0,1); # balanced ones on right

output="";
for (sample in (1:samples) ) 
{s = toString(sample);
for (imbalanced in order)
{

  for (i in (1:delay)) {x=sqrt(2)}
    work=(1-imbalanced)*c(10,20,30,10,20,30,10,20,30) + (imbalanced)*c(8,10,12, 18,20,22, 28,30,32);
    loss=B.age*(age-25) + round(B.work*work+rnorm(9,0,SD.e),1);
    d=data.frame(loss,age,work);
    my.lm <- lm(loss ~ work + age);
    summary(my.lm)
    b.work=toString(round(round(100 * my.lm$coefficients[2],2)/100,2));
    b.age =toString(round(round(100 * my.lm$coefficients[3],2)/100,2));
    a.lm <- lm(loss ~ age); w.lm <- lm(loss ~ work);
    b.age1=toString(round(round(100 * a.lm$coefficients[2],2)/100,2));
    b.work1=toString(round(round(100 * w.lm$coefficients[2],2)/100,2));
    title=paste("work & age:", b.work,"&",b.age,"[work:",b.work1," age:",b.age1,"]");
    if (imbalanced==0)  output=paste(output,s,"&", b.work,",",b.age," & ", b.work1, " & ", b.age1, "& &" )
    else                output=paste(output,       b.work,",",b.age," & ", b.work1, " & ", b.age1, "\\");
    myplot3d=scatterplot3d(work,age,loss, col.axis = "blue", zlim=c(5,40), xlim=c(5,35), 
    col.grid = "lightblue",  angle=55, type="h", highlight.3d=TRUE, 
    scale.y=0.7, pch=16,  main=title,cex.main=0.75,cex.axis=0.6)                 
    myplot3d$plane3d(my.lm)
  }
}