#  source("bsimptests.s")

#  Mar 30, 2001
#  Computes test statistics for testing
#     H0:  b1 = 0
#  in the simple regression model, which has been fit using
#  bisquare. 
#  ######### 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)).

bsimptests  <- function(bhat, resful, resred, sig, tunec=4.685){
                       
n        <- length(resful)


cvalue   <- tunec*sig
cvaluesq <- cvalue^2

psisq    <- ifelse(abs(resful)<cvalue, 
                     (resful*(1-(resful/cvalue)^2)^2)^2, 
                                   0)
psiprime <- ifelse(abs(resful)<cvalue, 
                   (1-5*(resful/cvalue)^2)*(1-(resful/cvalue)^2),0)
meanpsiprime <- mean(psiprime)
sumpsisq <- sum(psisq)
tauhat   <- sqrt(sumpsisq/n)/meanpsiprime

kappa    <- 1 + 2*(n-1)*var(psiprime)/(sum(psiprime))^2

#    Compute the Wald statistics

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

#    Compute fsh

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

rho <- ifelse(abs(resful)<cvalue, 
              resful^2/2-resful^4/(2*cvaluesq)+ resful^6/(6*cvalue^4),
              cvaluesq/6)
strful  <- sum(rho)

# Reduced model

rhor <- ifelse(abs(resred)<cvalue, 
              resred^2/2-resred^4/(2*cvaluesq)+ resred^6/(6*cvalue^4),
              cvaluesq/6)

strred  <- sum(rhor)

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

#zpvalue  <- 2*pnorm(-abs(zstat))
#zpvaluea <- 2*pnorm(-abs(zstata))
#tpvalue  <- 2*pt(-abs(zstat),n-2)
#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, fsh,fsha)

}