#  source("Hsimptests.s")

#  Mar 30, 2001
#  Computes test statistics for testing
#     H0:  b1 = 0
#  in the simple regression model, which has been fit using
#  Huber's method.
#  ######### Assumes var(x1)=1

#  The test statistics are:

#  * two Wald test statistics 
#  * fsh (Shrader and Hettsmansberger)

#  Normally, sig is computed as 1.483*median(abs(resful)).

Hsimptests  <- function(bhat, resful, resred, sig, kbend=1.345){
                       
n        <- length(resful)
kvaluesq <- (kbend*sig)^2
psisq    <- ifelse(resful^2<kvaluesq, resful^2, kvaluesq)
psiprime <- ifelse(resful^2<kvaluesq, 1,0)
meanpsiprime <- mean(psiprime)
sumpsisq <- sum(psisq)
tauhat   <- sqrt(sumpsisq/n)/meanpsiprime
kappa    <- 1 - 2/n + 2/sum(psiprime)

#    Compute the Wald statistics

zstat    <- bhat/tauhat
zstata   <- (sqrt(1-2/n)/kappa) * zstat

#    Compute fsh

lambda   <- sumpsisq/(2*(n-2)*meanpsiprime)

rhof <- ifelse(abs(resful)<kbend*sig,
                 resful^2/2, abs(resful)*kbend*sig- (kbend*sig)^2/2)
strful <- sum(rhof)

# Reduced model

rhor <- ifelse(abs(resred)<kbend*sig,
                 resred^2/2, abs(resred)*kbend*sig- (kbend*sig)^2/2)
strred <- sum(rhor)

fsh <- (strred-strful)/lambda
fsha <- (n/(n-2)) * fsh

#zpvalue  <- 2*pnorm(-abs(zstat))
#tpvalue  <- 2*pt(-abs(zstat),n-2)
#zpvaluea <- 2*pnorm(-abs(zstata))
#tpvaluea <- 2*pt(-abs(zstata),n-2)
#fmpvalue <- 1- pf(fsh, 1, n-2)
#return(zstat,zstata,fsh,zpvalue, tpvalue, zpvaluea, tpvaluea, fmpvalue)

return(zstat,zstata,strred,strful,lambda,fsh,fsha)
}