#  source("fullmonte.s")
#  updated 8/27/01

fullmonte  <- function(m, n, b, dist, design, seed)
{
#  Arguments:
#    m      - the Monte Carlo sample size
#    n      - the sample size
#    b      - the true value of the slope.  
#             The null hypothesis is b=0.
#    dist   - a character variable designating the distribution up to
#             a scale factor (All distributions as scaled so as to have
#             unit variance.)
#             "stdnorm" = Normal(0,1)
#             "t6"      = t with 6 df
#             "cont"    = .9N(0,s) + .1N(0,10s), where s = sqrt(1/10.9)
#    design - the pattern of x's 
#             design = 1 means equally spaced
#    seed   - an integer between 0 and 1023 or a negative integer.
#             if 0 <= seed <= 1023, seed is used in set.seed to change
#             the state of the generator.
#             if seed < 0, the state of the generator is under user
#             control to allow accumulation of results.
#             For user control, the user should set .Random.seed prior 
#             to invoking the full monty:
#                  .Random.seed <- myoldseed
#                  updates <- fullmonte(m, n, b, dist, design,-1)
#             and then after the run
#                  myoldseed  <- .Random.seed 
#             so as to be able to pick up from where he left off.

#  Tests of the hypothesis in simple linear regression
#
#    model:  y = b0 + b1*x + E
#
#  Hypothesis: b1=0

#  We will consider the hypothesis test using test statistics
#  derived from
#  4 methods of fitting:
#    * "ls": ordinary least squares  
#      -- lsfit
#    * "h1": M regression using Huber and MAD 
#      -- rlmx(.., scale.est="MAD")
#    * "h2": M regression using Huber's proposal 2 
#      -- rlmx(.., scale.est="proposal 2")
#    * "b1": M regression using a bisquare weight
#      -- rlmx(.., psi = psi.bisquare)

#  rlmx is a modification of Venables and Ripley's rlm.

#  For each of the
#  3 robust methods, we will consider 
#  4 test statistics:
#    * a Wald's t-type statistic, "zstat"
#    * a modified Wald's t-type statistic, "zstata"
#    * the F-type statistic of Shrader and Hettsmansberger, "fsh"
#    * an adjusted F-type statistic "fsha"
#    * for each of the Wald's statistics, we will consider 
#       two distributions:
#       * t
#       * normal

#   Thus, there are 19 tests.

if (seed >0) set.seed(seed)
kbend <- 1.345  ###### this is important, because it defines the fits
tunec <- 4.685

#   The critical values for the 3 distributions used for 
#   approximating those of the test statistics under the null.

zcrit05  <- qnorm(1-.05/2)
tcrit05  <- qt(1-.05/2,n-2)
fcrit05  <- qf(1-.05, 1, n-2)
zcrit01  <- qnorm(1-.01/2)
tcrit01  <- qt(1-.01/2,n-2)
fcrit01  <- qf(1-.01, 1, n-2)

#  x is a matrix with the first column of 1's and a scaled second column
if (design == 1) {
   x      <- c(1:n)
   x      <- x/sqrt((n-1)*var(x))
   xconst <- rep(1,n)
   x      <- cbind(xconst,x)
}

#  The 32 counters for the 32 tests:

nrejlst05   <- 0   # ls
nrejlst01   <- 0   

nrejh1z05   <- 0   # Huber MAD
nrejh1t05   <-0
nrejh1za05  <- 0
nrejh1ta05  <-0
nrejh1fsh05 <- 0
nrejh1fsha05 <- 0
nrejh1z01   <- 0 
nrejh1t01   <-0
nrejh1za01  <- 0
nrejh1ta01  <-0
nrejh1fsh01 <- 0
nrejh1fsha01 <- 0

nrejh2z05   <- 0   # Huber Proposal 2
nrejh2t05   <-0
nrejh2za05  <- 0
nrejh2ta05  <-0
nrejh2fsh05 <- 0
nrejh2fsha05 <- 0
nrejh2z01   <- 0 
nrejh2t01   <-0
nrejh2za01  <- 0
nrejh2ta01  <-0
nrejh2fsh01 <- 0
nrejh2fsha01 <- 0

nrejb1z05   <- 0   # M bisquare
nrejb1t05   <-0
nrejb1za05  <- 0
nrejb1ta05  <-0
nrejb1fsh05 <- 0
nrejb1fsha05 <- 0
nrejb1z01   <- 0 
nrejb1t01   <-0
nrejb1za01  <- 0
nrejb1ta01  <-0
nrejb1fsh01 <- 0
nrejb1fsha01 <- 0

###### Begin

for (i in 1:m) {
bad.Random.seed<-.Random.seed ### looking for problems
if (dist=='stdnorm') errors<-rnorm(n)
  else if (dist=='t6') errors<-rt(n,6)*sqrt(2/3)
  else if (dist=='cont') {
     errors<-rnorm(n)
     errors<-ifelse(runif(n)<.9, errors, 10*errors)/sqrt(10.9)
    }
    else return(message='Incorrect distribution specified')

y<-1+b*x[,2]+errors

lsfitf  <- lsfit(x[,2],y)

rlmxh1f <- rlmx(x,y,k=kbend,scale.est="MAD",k2=kbend)
rlmxh2f <- rlmx(x,y,k=kbend,scale.est="proposal 2",k2=kbend)
rlmxb1f <- rlmx(x,y,psi = psi.bisquare,c=tunec)

bhath1  <- rlmxh1f$coef[2]
bhath2  <- rlmxh2f$coef[2]
bhatb1  <- rlmxb1f$coef[2]

#   t statistic from ls fit
tlsfitf    <- lsfitf$coef[2]/sqrt((n-1)*var(lsfitf$residuals)/(n-2))

#   t statistics computed in rlmx  *** we're not using these currently
trlmxh1f   <- summary(rlmxh1f)$coef[2,3]
trlmxh2f   <- summary(rlmxh2f)$coef[2,3]
trlmxb1f   <- summary(rlmxb1f)$coef[2,3]

sigrlmxh1f  <- summary(rlmxh1f)$sigma
sigrlmxh2f  <- summary(rlmxh2f)$sigma
sigrlmxb1f  <- summary(rlmxb1f)$sigma

rrlmxh1f <- summary(rlmxh1f)$residuals
rrlmxh2f <- summary(rlmxh2f)$residuals
rrlmxb1f <- summary(rlmxb1f)$residuals

rlmxh1r <- rlmx(x[,1],y,k=kbend,const.scale = sigrlmxh1f,k2=kbend)
rlmxh2r <- rlmx(x[,1],y,k=kbend,const.scale = sigrlmxh2f,k2=kbend)
rlmxb1r <- rlmx(x[,1],y,psi = psi.bisquare,const.scale = sigrlmxb1f,c=tunec)

rrlmxh1r <- summary(rlmxh1r)$residuals
rrlmxh2r <- summary(rlmxh2r)$residuals
rrlmxb1r <- summary(rlmxb1r)$residuals

hhh1 <- Hsimptests(bhath1, rrlmxh1f, rrlmxh1r, sigrlmxh1f, kbend=kbend)
hhh2 <- Hsimptests(bhath2, rrlmxh2f, rrlmxh2r, sigrlmxh2f, kbend=kbend)
bbb1 <- bsimptests(bhatb1, rrlmxb1f, rrlmxb1r, sigrlmxb1f, tunec=tunec)
#############################################
if (hhh1$fsh<0) cat("Huber 1",bad.Random.seed,"\n")
if (hhh2$fsh<0)cat("Huber 2",bad.Random.seed,hhh2,"\n")
if (bbb1$fsh<0) cat("bisquare",bad.Random.seed,"\n")
#  Update the 32 counters

if (abs(tlsfitf)>tcrit05)     nrejlst05  <- nrejlst05 + 1
if (abs(hhh1$zstat)>zcrit05)  nrejh1z05   <- nrejh1z05 + 1
if (abs(hhh1$zstat)>tcrit05)  nrejh1t05   <- nrejh1t05 + 1
if (abs(hhh1$zstata)>zcrit05) nrejh1za05  <- nrejh1za05 + 1
if (abs(hhh1$zstata)>tcrit05) nrejh1ta05  <- nrejh1ta05 + 1
if (hhh1$fsh>fcrit05)         nrejh1fsh05 <- nrejh1fsh05 + 1
if (hhh1$fsha>fcrit05)        nrejh1fsha05<- nrejh1fsha05 + 1
if (abs(hhh2$zstat)>zcrit05)  nrejh2z05   <- nrejh2z05 + 1
if (abs(hhh2$zstat)>tcrit05)  nrejh2t05   <- nrejh2t05 + 1
if (abs(hhh2$zstata)>zcrit05) nrejh2za05  <- nrejh2za05 + 1
if (abs(hhh2$zstata)>tcrit05) nrejh2ta05  <- nrejh2ta05 + 1
if (hhh2$fsh>fcrit05)         nrejh2fsh05 <- nrejh2fsh05 + 1
if (hhh2$fsha>fcrit05)        nrejh2fsha05<- nrejh2fsha05 + 1
if (abs(bbb1$zstat)>zcrit05)  nrejb1z05   <- nrejb1z05 + 1
if (abs(bbb1$zstat)>tcrit05)  nrejb1t05   <- nrejb1t05 + 1
if (abs(bbb1$zstata)>zcrit05) nrejb1za05  <- nrejb1za05 + 1
if (abs(bbb1$zstata)>tcrit05) nrejb1ta05  <- nrejb1ta05 + 1
if (bbb1$fsh>fcrit05)         nrejb1fsh05 <- nrejb1fsh05 + 1
if (bbb1$fsha>fcrit05)        nrejb1fsha05<- nrejb1fsha05 + 1

if (abs(tlsfitf)>tcrit01)     nrejlst01  <- nrejlst01 + 1
if (abs(hhh1$zstat)>zcrit01)  nrejh1z01   <- nrejh1z01 + 1
if (abs(hhh1$zstat)>tcrit01)  nrejh1t01   <- nrejh1t01 + 1
if (abs(hhh1$zstata)>zcrit01) nrejh1za01  <- nrejh1za01 + 1
if (abs(hhh1$zstata)>tcrit01) nrejh1ta01  <- nrejh1ta01 + 1
if (hhh1$fsh>fcrit01)         nrejh1fsh01 <- nrejh1fsh01 + 1
if (hhh1$fsha>fcrit01)        nrejh1fsha01<- nrejh1fsha01 + 1
if (abs(hhh2$zstat)>zcrit01)  nrejh2z01   <- nrejh2z01 + 1
if (abs(hhh2$zstat)>tcrit01)  nrejh2t01   <- nrejh2t01 + 1
if (abs(hhh2$zstata)>zcrit01) nrejh2za01  <- nrejh2za01 + 1
if (abs(hhh2$zstata)>tcrit01) nrejh2ta01  <- nrejh2ta01 + 1
if (hhh2$fsh>fcrit01)         nrejh2fsh01 <- nrejh2fsh01 + 1
if (hhh2$fsha>fcrit01)        nrejh2fsha01<- nrejh2fsha01 + 1
if (abs(bbb1$zstat)>zcrit01)  nrejb1z01   <- nrejb1z01 + 1
if (abs(bbb1$zstat)>tcrit01)  nrejb1t01   <- nrejb1t01 + 1
if (abs(bbb1$zstata)>zcrit01) nrejb1za01  <- nrejb1za01 + 1
if (abs(bbb1$zstata)>tcrit01) nrejb1ta01  <- nrejb1ta01 + 1
if (bbb1$fsh>fcrit01)         nrejb1fsh01 <- nrejb1fsh01 + 1
if (bbb1$fsha>fcrit01)        nrejb1fsha01<- nrejb1fsha01 + 1

}
output <- matrix(rep(0,38),ncol=2)
output[1,1] <- nrejlst01/m
output[2,1] <- nrejh1z01/m
output[3,1] <- nrejh1t01/m
output[4,1] <- nrejh1za01/m
output[5,1] <- nrejh1ta01/m
output[6,1] <- nrejh1fsh01/m
output[7,1] <- nrejh1fsha01/m
output[8,1] <- nrejh2z01/m
output[9,1] <- nrejh2t01/m
output[10,1] <- nrejh2za01/m
output[11,1] <- nrejh2ta01/m
output[12,1] <- nrejh2fsh01/m
output[13,1] <- nrejh2fsha01/m
output[14,1] <- nrejb1z01/m
output[15,1] <- nrejb1t01/m
output[16,1] <- nrejb1za01/m
output[17,1] <- nrejb1ta01/m
output[18,1] <- nrejb1fsh01/m
output[19,1] <- nrejb1fsha01/m

output[1,2] <- nrejlst05/m
output[2,2] <- nrejh1z05/m
output[3,2] <- nrejh1t05/m
output[4,2] <- nrejh1za05/m
output[5,2] <- nrejh1ta05/m
output[6,2] <- nrejh1fsh05/m
output[7,2] <- nrejh1fsha05/m
output[8,2] <- nrejh2z05/m
output[9,2] <- nrejh2t05/m
output[10,2] <- nrejh2za05/m
output[11,2] <- nrejh2ta05/m
output[12,2] <- nrejh2fsh05/m
output[13,2] <- nrejh2fsha05/m
output[14,2] <- nrejb1z05/m
output[15,2] <- nrejb1t05/m
output[16,2] <- nrejb1za05/m
output[17,2] <- nrejb1ta05/m
output[18,2] <- nrejb1fsh05/m
output[19,2] <- nrejb1fsha05/m

dimnames(output)<-list(c(    "           t",
                             "      MAD Hz",
                             "      MAD Ht",
                             "     MAD Haz",
                             "     MAD Hat",
                             "      MAD HF",
                             "     MAD HFa",
                             "  Huber 2 Hz",
                             "  Huber 2 Ht",
                             " Huber 2 Haz",
                             " Huber 2 Hat",
                             "  Huber 2 HF",
                             " Huber 2 HFa",
                             " bisquare Hz",
                             " bisquare Ht",
                             "bisquare Haz",
                             "bisquare Hat",
                             " bisquare HF",
                             "bisquare HFa"),
                        c("  0.01  ","0.05  "))
return(output)
}