saws<-function(beta,u,omega,test=1:p,alpha=0.05,output=1){
### calculate Small Sample Adjustments for Wald-type tests using Sandwich estimators (saws)
### input:
###   u = k by p matrix
###  omega = k by p by p array
###  beta = p by 1 matrix or vector
###  test = q by 1 matrix of indices, where q<=p, 
###      indices denote parameters from 1 to p
###       AS OF NOW THIS ONLY TESTS beta[j]=0 or not for j in test
###  alpha = nominal significance level (default=.05) Note:p-values are 2-sided
###  output = 
###     0 only output recommended test (Sandwich using F distn
###            with degrees of freedom=1, d.tilde.no.H)
###     1 (default) output Model based, Standard Sandwich and Sandwich using F distn
###            with degrees of freedom=1, d.tilde.no.H
###     2 output all tests
beta.names<-names(beta)
beta<-as.vector(beta)
p<-length(beta)

beta<-matrix(beta,p,1)
K<-dim(u)[[1]]

	vminv<- apply(omega,c(2,3),sum)
	vm<- solve(  vminv )
  vs<- vm %*% (t(u) %*% u) %*% vm

### calculate pre- the correction matrix that premultiplies U*_i
### and the adjusted Var - va
H<-array(NA,c(K,p,p))
Psi.tilde.with.H<-matrix(0,K*p,K*p)
Psi.hat.with.H<-matrix(0,K*p,K*p)

psi.hat.with.H.sum<-matrix(0,p,p)

Psi.tilde.no.H<-matrix(0,K*p,K*p)
Psi.hat.no.H<-matrix(0,K*p,K*p)

psi.hat.no.H.sum<-matrix(0,p,p)

G<-matrix(0,p*K,p*K)

for (i in 1:K){
	omegaivm<-omega[i,,] %*% vm
   H[i,,]<- diag( 1/sqrt( 1- pmin(.75,diag( omegaivm ))  ) )

	index<- ((i-1)*p + 1): (i*p)
   Psi.hat.with.H[index,index]<-H[i,,] %*%  matrix(u[i,],p,1) %*% matrix(u[i,],1,p) %*% H[i,,]
   psi.hat.with.H.sum<-psi.hat.with.H.sum+Psi.hat.with.H[index,index]
   Psi.hat.no.H[index,index]<- matrix(u[i,],p,1) %*% matrix(u[i,],1,p)
   psi.hat.no.H.sum<-psi.hat.no.H.sum+Psi.hat.no.H[index,index]

    G[index,index]<- diag(p) 
    G[index,]<- G[index,] - matrix(rep(omegaivm,K),p,p*K)
}




## calculate d.tilde.with.H and d.hat.with.H
d.tilde.with.H<-rep(0,p)
d.hat.with.H<-rep(0,p)

d.tilde.no.H<-rep(0,p)
d.hat.no.H<-rep(0,p)

for (j in test){
 D.with.H<-matrix(0,p*K,p*K)
 D.no.H<-matrix(0,p*K,p*K)

 Psi.tilde.with.H<-matrix(0,K*p,K*p)

 wi<-rep(0,K)
 for (i in 1:K){

     index<-	((i-1)*p + 1): (i*p)
     D.with.H[index,index]<- t(H[i,,]) %*% matrix(vm[,j],p,1) %*% matrix(vm[j,],1,p) %*% (H[i,,])
     D.no.H[index,index]<-  matrix(vm[,j],p,1) %*% matrix(vm[j,],1,p) 

	wi[i]<- ( solve( vminv - omega[i,,]) - vm )[j,j]
	Psi.tilde.with.H[index,index]<- wi[i]* (psi.hat.with.H.sum)
	Psi.tilde.no.H[index,index]<- wi[i]* (psi.hat.no.H.sum)

}
Psi.tilde.with.H<-Psi.tilde.with.H/sum(wi)
Psi.tilde.no.H<-Psi.tilde.no.H/sum(wi)

### to save memory we use Psi# to store a big matrix
Psi.tilde.with.H<- Psi.tilde.with.H %*% t(G) %*% D.with.H %*% G
Psi.hat.with.H<- Psi.hat.with.H %*% t(G) %*% D.with.H %*% G
Psi.tilde.no.H<- Psi.tilde.no.H %*% t(G) %*% D.no.H %*% G
Psi.hat.no.H<- Psi.hat.no.H %*% t(G) %*% D.no.H %*% G

### use properties of eigenvalues to  get    d[j]<- trace(Psi)^2 /  trace(Psi %*% Psi )
d.tilde.with.H[j]<- ( sum( diag(Psi.tilde.with.H) )**2 )/ sum(as.vector(Psi.tilde.with.H)*as.vector(t(Psi.tilde.with.H)))
d.hat.with.H[j]<- ( sum( diag(Psi.hat.with.H) )**2 )/ sum(as.vector(Psi.hat.with.H)*as.vector(t(Psi.hat.with.H)))

d.tilde.no.H[j]<- ( sum( diag(Psi.tilde.no.H) )**2 )/ sum(as.vector(Psi.tilde.no.H)*as.vector(t(Psi.tilde.no.H)))
d.hat.no.H[j]<- ( sum( diag(Psi.hat.no.H) )**2 )/ sum(as.vector(Psi.hat.no.H)*as.vector(t(Psi.hat.no.H)))

}

## now do sandwiching 
psi.hat.with.H.sum<- vm %*% psi.hat.with.H.sum %*% vm


### Test names=type of variance, distribution of (Wald)^2 stat, degrees of freedom for F
###  Chi2=chi squared with 1 degree of freedom
test.names<-c("Vm,Chi2","Vs,Chi2","Vs,F,K-p","Vs,F,d.tilde.no.H","Vs,F,d.hat.no.H","Va,Chi2","Va,F,d.tilde.with.H",
"Va,F,d.hat.with.H")

numM<-length(test.names)

   twosided.pvalues<-matrix(NA,length(test),numM)
   ci.length<-matrix(NA,length(test),numM)
	rse <- sqrt(diag(vs))[test]
	arse1 <- sqrt(diag(psi.hat.with.H.sum))[test]
    nse <- sqrt(diag(vm))[test]
    coef<-beta
    beta<-beta[test]

   ## for each element of the sd vectors, take the maximum 
   	maxvec<-function(x,y){ x[y>x]<-y[y>x] ; x }
	maxan1<-maxvec(arse1,nse)
	
	maxrn<-maxvec(rse,nse)
	
	
	df.tilde.with.H<-d.tilde.with.H[test]
	df.hat.with.H<-d.hat.with.H[test]
	
	df.tilde.no.H<-d.tilde.no.H[test]
	df.hat.no.H<-d.hat.no.H[test]
	
	
	
		### Calculate P-values 
	# Vm,Chi2 = model-based  
	twosided.pvalues[,1]<-  1 - pchisq( (beta/nse)**2,1)
	# Vs,Chi2 = naive sandwich
	twosided.pvalues[,2]<-  1 - pchisq( (beta/rse)**2,1)
	# Vs,F,K-p = naive sandwich with F distrn
	twosided.pvalues[,3]<-  1 - pf( (beta/rse)**2,1,K-p)
	# Vs,F,d.tilde.no.H = naive sandwich with F distrn
	twosided.pvalues[,4]<-  1 - pf( (beta/rse)**2,1,df.tilde.no.H)
	# Vs,F,d.hat.no.H = naive sandwich with F distrn
	twosided.pvalues[,5]<-  1 - pf( (beta/rse)**2,1,df.hat.no.H)
	# Va,Chi2 =  adjusted sandwich var using Chi square
	twosided.pvalues[,6]<- 1- pchisq( (beta/arse1)**2,1)
	# Va,F,d.tilde.with.H = adjusted sandwich var using F distribution
	twosided.pvalues[,7]<- 1- pf( (beta/arse1)**2,1,df.tilde.with.H)
	# Va,F,d.hat.with.H = adjusted sandwich var using F distribution
	twosided.pvalues[,8]<- 1- pf( (beta/arse1)**2,1,df.hat.with.H)





 dimnames(twosided.pvalues)<-list(NULL,test.names)

	### Calculate Confidence Interval Lengths
	
	
	# Vm,Chi2 = model-based  
	ci.length[,1]<-  sqrt(qchisq(1-alpha,1))*nse
	# Vs,Chi2 = naive sandwich
	ci.length[,2]<-  sqrt(qchisq(1-alpha,1))*rse
	# Vs,F,K-p = naive sandwich with F distrn
	ci.length[,3]<-  sqrt(qf(1-alpha,1,K-p))*rse
	# Vs,F,d.tilde.no.H = naive sandwich with F distrn
	ci.length[,4]<-  sqrt(qf(1-alpha,1,df.tilde.no.H))*rse
	# Vs,F,d.hat.no.H = naive sandwich with F distrn
	ci.length[,5]<-  sqrt(qf(1-alpha,1,df.hat.no.H))*rse
	# Va,Chi2 =  adjusted sandwich var using Chi square
	ci.length[,6]<-  sqrt(qchisq(1-alpha,1))*arse1
	# Va,Chi2 =  adjusted sandwich var using F,d.tilde.with.H
	ci.length[,7]<-sqrt(qf(1-alpha,1,df.tilde.with.H))*arse1
	# Va,Chi2 =  adjusted sandwich var using F,d.hat.with.H
	ci.length[,8]<-sqrt(qf(1-alpha,1,df.hat.with.H))*arse1

 dimnames(ci.length)<-list(NULL,test.names)


### bundle up results into matrices for output
	Vm.Chi2<-matrix(c(beta,beta-ci.length[,1],beta+ci.length[,1],
	twosided.pvalues[,1]),length(beta[test]),4,
	dimnames=list(beta.names,c("Coefficients",paste("Lower ",100*(1-alpha),"% C.L."),
	paste("Upper ",100*(1-alpha),"% C.L."),"2-Sided P-value")))


	Vs.Chi2<-matrix(c(beta,beta-ci.length[,2],beta+ci.length[,2],
	twosided.pvalues[,2]),length(beta[test]),4,
	dimnames=list(beta.names,c("Coefficients",paste("Lower ",100*(1-alpha),"% C.L."),
	paste("Upper ",100*(1-alpha),"% C.L."),"2-Sided P-value")))

	Vs.F.K.minus.p<-matrix(c(beta,beta-ci.length[,3],beta+ci.length[,3],rep(K-p,length(beta[test])),
	twosided.pvalues[,3]),length(beta[test]),5,
	dimnames=list(beta.names,c("Coefficients",paste("Lower ",100*(1-alpha),"% C.L."),
	paste("Upper ",100*(1-alpha),"% C.L."),"DF estimates","2-Sided P-value")))

	Vs.F.d.tilde.no.H<-matrix(c(beta,beta-ci.length[,4],beta+ci.length[,4],df.tilde.no.H,
	twosided.pvalues[,4]),length(beta[test]),5,
	dimnames=list(beta.names,c("Coefficients",paste("Lower ",100*(1-alpha),"% C.L."),
	paste("Upper ",100*(1-alpha),"% C.L."),"DF estimates","2-Sided P-value")))

	Vs.F.d.hat.no.H<-matrix(c(beta,beta-ci.length[,5],beta+ci.length[,5],df.hat.no.H,
	twosided.pvalues[,5]),length(beta[test]),5,
	dimnames=list(beta.names,c("Coefficients",paste("Lower ",100*(1-alpha),"% C.L."),
	paste("Upper ",100*(1-alpha),"% C.L."),"DF estimates","2-Sided P-value")))

	Va.Chi2<-matrix(c(beta,beta-ci.length[,6],beta+ci.length[,6],
	twosided.pvalues[,6]),length(beta[test]),4,
	dimnames=list(beta.names,c("Coefficients",paste("Lower ",100*(1-alpha),"% C.L."),
	paste("Upper ",100*(1-alpha),"% C.L."),"2-Sided P-value")))

	Va.F.d.tilde.with.H<-matrix(c(beta,beta-ci.length[,7],beta+ci.length[,7],df.tilde.with.H,
	twosided.pvalues[,7]),length(beta[test]),5,
	dimnames=list(beta.names,c("Coefficients",paste("Lower ",100*(1-alpha),"% C.L."),
	paste("Upper ",100*(1-alpha),"% C.L."),"DF estimates","2-Sided P-value")))

	Va.F.d.hat.with.H<-matrix(c(beta,beta-ci.length[,8],beta+ci.length[,8],df.hat.with.H,
	twosided.pvalues[,8]),length(beta[test]),5,
	dimnames=list(beta.names,c("Coefficients",paste("Lower ",100*(1-alpha),"% C.L."),
	paste("Upper ",100*(1-alpha),"% C.L."),"DF estimates","2-Sided P-value")))


if (output==0){
	out<-list(Vs.F.d.tilde.no.H=Vs.F.d.tilde.no.H)
}
else if (output==1){
		out<-list(Vm.Chi2=Vm.Chi2,Vs.Chi2=Vs.Chi2,Vs.F.d.tilde.no.H=Vs.F.d.tilde.no.H)
}
else {
	out<-list(Vm.Chi2=Vm.Chi2,
	Vs.Chi2=Vs.Chi2,
	Vs.F.K.minus.p=Vs.F.K.minus.p,
	Vs.F.d.tilde.no.H=Vs.F.d.tilde.no.H,
	Vs.F.d.hat.no.H=Vs.F.d.hat.no.H,
	Va.Chi2=Va.Chi2,
	Va.F.d.tilde.with.H=Va.F.d.tilde.with.H,
	Va.F.d.hat.with.H=Va.F.d.hat.with.H)
	}
out
}