##### R-code for Example 8.7. It writes the related data to the file
##### "example87.dat" and makes "fig89.ps"

set.seed(100)

#### Generate 50 observations from the t(3)/sqrt(3) distribution, 
#### and another 50 observations from t(3)/sqrt(3)+1

x1=rt(50,3)/sqrt(3)
x2=rt(50,3)/sqrt(3)+1

x = c(x1,x2)
n = length(x)

#### Compute T_{max,n}, for 15 <= n <= 100

Tmaxn = rep(0,100) 
Ujn = matrix(0,n,n) 
Tjn = matrix(0,n,n) 

for(j in 1:14){       # Compute Ujn[j,15], for 1 <= j <= 14.

   Ujn[j,15]=0

   for(i in 1:j){
      for(i1 in (j+1):15){
         Ujn[j,15]=Ujn[j,15]+as.numeric(x[i] > x[i1])
         Ujn[j,15]=Ujn[j,15]-as.numeric(x[i] < x[i1])
      }
   }

   Tjn[j,15]=Ujn[j,15]/sqrt(j*(15-j)*16/3)

}

Tmaxn[15]=max(abs(Tjn[1:14,15]))

for(t in 16:n){

   for(j in 1:(t-1)){  # Compute Ujn[j,t], for 1 <= j <= t-1.

      Ujn[j,t]=Ujn[j,t-1]
      for(i in 1:j){
         Ujn[j,t]=Ujn[j,t]+as.numeric(x[i] > x[t])
         Ujn[j,t]=Ujn[j,t]-as.numeric(x[i] < x[t])
      }

      Tjn[j,t]=Ujn[j,t]/sqrt(j*(t-j)*(t+1)/3)

   }

   Tmaxn[t]=max(abs(Tjn[1:(t-1),t]))
}

write.table(round(cbind(x,Tmaxn),digits=3), 
            file="example87.dat",col.names=T, row.names=F)

postscript("fig89.ps",width=7,height=3.5,horizontal=F)

par(mfrow=c(1,2),mar=c(5,4,1,1))

ii = seq(1,n)

plot(ii,x,type="p",pch=16,xlab="n",ylab=expression(X[n]),mgp=c(2,1,0),
     xlim=c(0,101), ylim=c(-5,4.1), cex=0.8)
title(xlab="(a)",cex=0.9)

hn=c(2.947,2.910,2.862,2.860,2.869,rep(2.851,2),rep(2.862,2),rep(2.870,2),rep(2.875,2),
     rep(2.883,2),rep(2.879,5),rep(2.894,5),rep(2.900,5),rep(2.906,5),rep(2.908,10),
     rep(2.914,10),rep(2.917,10),rep(2.918,10),rep(2.920,10),2.922)

plot(ii[15:n],Tmaxn[15:n],type="o",lty=1,pch=16,xlab="n",ylab=expression(T[list(max,n)]),
     mgp=c(2,1,0),xlim=c(0,101), ylim=c(0,7),cex=0.8)
lines(ii[15:n],hn,lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

graphics.off()