# >>>> Tests taken from Gossner and Schlag (2011),
# Finite Sample Nonparametric Tests for Linear Regressions
# - implemented for the two unbiased estimates OLS and minimax
# - requires at least two covariates (independent variables)

# notes: 
# (i) remember to introduce a constant covariate unless you have 
# 	good reasons to impose that the regression line should go through the origin
# (ii) fixed effects if necessary need to be introduced manually in the data set
#      (the program cannot accommodate for random effects)
# (iii) needs packages quadprog, Rlab and Rglpk (version 0.3-5)

# code written by Gareth Liu-Evans and Karl Schlag (date Febuary 7, 2012)

# for testing H0: betaj <= betabarj against H1: betaj > betabarj, or 
# for testing H0: betaj >= betabarj against H1: betaj > betabarj 
# where j is index of the coefficient, depending on IE
IE <- ">="  # determines which inequality is used under H0 (insert "<=" or ">=")

alpha <- 0.025 # (enter) significance level of test
j <- 3 # (enter) j is the coefficient to be tested
betabarj <- 0.185 # (enter) betabarj is the critical value for the null hypothesis

betaj1 <- betabarj - 0.268 # (enter)
# the type II errors are evaluated and compared within this softare 
# for betaj>=betaj1 (so need betaj1 > betabarj) when IE = "<="
# for betaj<=betaj1 (so need betaj1 < betabarj) when IE = ">="


## setwd("C:/Users/Schlag/Documents/karl/gareth/dufloetal") 
#         [this is the directory where your data files are stored, such as T4A2Y.dat and T4A2Y.dat]

w1Y <- 0 #(enter) minimal possible value for Y as defined by setting, not as determined by actual data
w2Y <- 1 #(enter) maximal possible value for Y as defined by setting, not as determined by actual data

Y <- as.matrix(read.table("T4A6Y.dat"))
X <- as.matrix(read.table("T4A6X.dat")) 

# using values with headers that are separated by comma
# Y <- as.matrix(read.table("TJendogenous.dat",header=TRUE))
# X <- as.matrix(read.table("TJexogenous.dat",header=TRUE,sep=","))

# using data that comes from stata format
# library(foreign)
# Y <- as.matrix(read.dta("JL Y.dta"))
# X <- as.matrix(read.dta("JL X.dta")) 

# to obtain the lower endpoint of the 95% confidence interval 
# let alpha <- 0.025, set IE <- "<="
# and find smallest value of betabarj such that null hypothesis cannot be rejected,
# this value of betabarj is the lower endpoint

# to obtain the upper endpoint of the 95% confidence interval 
# let alpha <- 0.025 set IE <- ">="
# and find largest value of betabarj such that null hypothesis cannot be rejected,
# this value of betabarj is the upper endpoint

## source("Nonstandardized_and_Bernoullitests_v031.r") 
#         [this is the command to start the program]


# ......... fine tuning:

# weight in Bernoulli test put on d that equals upper limit of interval
# do not change unless error message
lambda <- 1  # OLS (default=1)
lambdamm <- 1  # MM (default=1)
monte <- 1000 # number of loops in the Bernoulli test, increase to get more accuracy (default=1000)

# parameters to ensure positive definiteness
qq <- 0.0001 #OLS (default=0.0001)
qqmm <- 0.0001 #MM (default=0.0001)

# ------------------------------------------------------------------------
# do not change anything below here!
# ------------------------------------------------------------------------

if (min(Y)<w1Y | max(Y)>w2Y) {
print("your values of Y are not within the range [w1Y,w2Y]!")
STOP
}

if (ncol(X)==1) {
print("you need at least two covariates, add a constant")
STOP
}

# check if X has full rank
if (det(t(X)%*%X)==0) {
print("columns are not independent, need to have full rank, eliminate some") 
STOP
}

# adjustments to deal with both inequalities

if (IE == "<=") {
w11Y <- w1Y
w22Y <- w2Y 
IEc <- 1
IEa <- ">="
betabarj0 <- betabarj
betaj11 <- betaj1
if (betaj1 <= betabarj) {
print ("you need to choose betaj1 > betabarj")
STOP }
} else if (IE == ">=") {
w11Y <- -w2Y
w22Y <- -w1Y 
Y <- -1* Y
IEc <- -1
IEa <- "<="
betabarj0 <- -1*betabarj
betaj11 <- betabarj0+betabarj-betaj1
if (betaj1 >= betabarj) {
print ("you need to choose betaj1 < betabarj")
STOP }
} else {
print ("you need to choose IE either as <= or >=")
STOP 
}


XROWS <- nrow(X)
XCOLS <- ncol(X)
YROWS <- XROWS


cj <- j

Summaryinitial<- c("alpha ", "coefficient index (j in the paper)", "n","m", alpha, j, XROWS, XCOLS)
dim(Summaryinitial)<- c(4,2)

mydatainitial<- as.data.frame(Summaryinitial)
names(mydatainitial)<- c(" ", "Data input")
print(mydatainitial)
 
cat("\n NON-STANDARDIZED TEST \n")

Y <- Y /(w22Y-w11Y) # puts Y in [ww,ww+1] where ww = w11Y/(w22Y-w11Y)
X <- X /(w22Y-w11Y) # rescales X so that beta remains unchanged
ww <- w11Y /(w22Y-w11Y)

Betahat <- solve((crossprod(X)))%*%crossprod(X,Y)
Betahatj <- Betahat[j]

Tau <- X%*%solve(crossprod(X))
Tau_j <- Tau[,j]
TAUj_2 <- sum((Tau_j)^2)
TAUj_inf <- max(abs(Tau_j))


# Type 1 OLS

ej <- rep(0,times=XCOLS)
ej[j]<- 1


# adjusting tau_j st Dmat positive definite, 
# rescale and then eliminate 0 entries
# increases bound so test remains exact

Tau_jB <- rep(0,times=XROWS)
i <- 1
while (i <= XROWS) {
if (Tau_j[i]>=0) 
Tau_jB[i] <- max(Tau_j[i]/min(0.01,TAUj_inf),qq*min(1,TAUj_inf)) 
else Tau_jB[i] <- min(Tau_j[i]/min(0.01,TAUj_inf),-qq*min(1,TAUj_inf)) 
i <- i+1
}

DmatOLS <- 2*t(X)%*%((Tau_jB^2)*X)
dvecOLS <- (1+2*ww)*colMeans((Tau_jB^2)*X)*XROWS
Amat <- t(matrix(rbind(-t(ej), X, -X),ncol=XCOLS))
bvec0 <- as.vector(c(-betabarj0, rep(ww, times=XROWS), rep(-(ww+1), times=XROWS)))


if (det(DmatOLS)<=10^-15) {

print(paste("Dmat not positive definite for OLS, Det=",det(DmatOLS)))
sigmasqbar_beta0jOLSvalue <- TAUj_2/4-(1/XROWS)*(min(0,betabarj0-(1/2)*sum(Tau_j)))^2

} else {

sigmasqbar_beta0jOLSQP <- solve.QP(Dmat=DmatOLS, dvec=dvecOLS, Amat, bvec=bvec0, meq=0, factorized=FALSE)

zsol0OLS <- sigmasqbar_beta0jOLSQP$solution

sigmasqbar_beta0jOLSvalue <- sum((Tau_j)^2*((X%*%zsol0OLS-ww)*(1+ww-X%*%zsol0OLS)))
}

#
#
#
#
#
#
# Hoeffding and Cantelli OLS
#
#
#
#
#
#


Hoeftbar <- (-0.5*TAUj_2*log(alpha))^0.5

Cantellitbar <- ((sigmasqbar_beta0jOLSvalue)*(1-alpha)/alpha)^(0.5)

#
#
#
#
# Bhattacharyya OLS
#
#
#
#
# 

if (Cantellitbar*(Cantellitbar - TAUj_inf)/sigmasqbar_beta0jOLSvalue > 1) {

Bhattbar <- sqrt(sigmasqbar_beta0jOLSvalue * (1 + sqrt( 3 *(1-alpha) / alpha ) ) )

if (Bhattbar^2 * TAUj_inf < sigmasqbar_beta0jOLSvalue * ( TAUj_inf + 3 * Bhattbar)) {

r <- uniroot(function(x) ((3*sigmasqbar_beta0jOLSvalue-TAUj_inf^2)*sigmasqbar_beta0jOLSvalue)/((3*sigmasqbar_beta0jOLSvalue-TAUj_inf^2)*(sigmasqbar_beta0jOLSvalue+x^2)+(x^2-x*TAUj_inf-sigmasqbar_beta0jOLSvalue)^2) - alpha, c(0,Bhattbar))

Bhattbar <- r$root
}

} else Bhattbar <- Inf

#
#
#
#
#
# Berry-Esseen OLS
#
#
#
# 


h<- function(wb1){
htemp <- 1000*(1-pnorm((tbartemp-wb1[2])/(sigmasqbar_beta0jOLSvalue+wb1[1]^2)^0.5)+0.56*2*TAUj_inf/((27^0.5)*wb1[1]))/pnorm(wb1[2]/wb1[1])
htemp
}


wb1start <- c(0.1, 0.1)
tbartemp <- Hoeftbar
BEtbar <- Inf
grid <- 0.00001
lowerBE <- c(0.000001,0.000001)

BEtype1 <- optim(wb1start, h, gr = NULL, 
      method = "L-BFGS-B", 
      lower = lowerBE, upper = c(10,10), control = list(fnscale=1))

BE <- (BEtype1$value/1000 < alpha)

while(BE){
BEtype1 <- optim(wb1start, h, gr = NULL, 
      method = "L-BFGS-B", 
      lower = lowerBE, upper = c(10,10), control = list(fnscale=1))

if (BEtype1$value/1000 > alpha) (BE <- FALSE) 
else {
BEtbar <- tbartemp
# print(BEtbar)
tbartemp <- tbartemp - grid
wb1start <- BEtype1$par
}
}

#
# 
#
#
#
#Final results for Non-Standardized test OLS cutoffs
#
#
#
#
#

tbarmin <- min(c(Cantellitbar, Bhattbar, Hoeftbar, BEtbar))

Summary1<- c("betabarj", "Betahatj", "Cantellitbar", "Bhattbar", "Hoeftbar", "BEtbar", round(betabarj, digits=5), round(IEc*Betahatj, digits=5), round(Cantellitbar, digits=5), round(Bhattbar, digits=5), round(Hoeftbar, digits=5), round(BEtbar, digits=5))
dim(Summary1)<- c(6,2)

mydata<- as.data.frame(Summary1)
names(mydata)<- c("statistic name (OLS version)", "value")
print(mydata)

if(Betahatj-betabarj0>=tbarmin) {
cat("Nonstandardized test >>> REJECTS H_0\n\n") 
RejectNonstandardizedOLS <- "YES"} else {
cat("Non-standardized test does not reject H_0\n\n")
RejectNonstandardizedOLS <- "NO" }

# Type 2 OLS

if (det(DmatOLS)<=10^-15) {

sigmasqbar_betajOLSvalue <- TAUj_2/4-(1/XROWS)*(betaj11-(1/2)*sum(Tau_j))^2

} else {

bvec1 <- as.vector(c(-betaj11, rep(ww, times=XROWS), rep(-(ww+1), times=XROWS)))

sigmasqbar_betajOLSQP <- solve.QP(Dmat=DmatOLS, dvec=dvecOLS, Amat, bvec=bvec1, meq=1, factorized=FALSE)

zsolOLS <- sigmasqbar_betajOLSQP$solution

sigmasqbar_betajOLSvalue <- sum((Tau_j)^2*((X%*%zsolOLS-ww)*(1+ww-X%*%zsolOLS)))
}


gOLS <- function(Bhattbar) {

if(((Bhattbar^2)/sigmasqbar_betajOLSvalue)-Bhattbar*TAUj_inf/sigmasqbar_betajOLSvalue - 1 <=0) gOLStemp <- 1 else
if(sigmasqbar_betajOLSvalue < (Bhattbar^2)*TAUj_inf/(TAUj_inf+3*Bhattbar))
gOLStemp <- (3*sigmasqbar_betajOLSvalue^2)/(4*(sigmasqbar_betajOLSvalue^2)-2*sigmasqbar_betajOLSvalue*Bhattbar^2+Bhattbar^4)
else gOLStemp <- ((3*sigmasqbar_betajOLSvalue-TAUj_inf^2)*sigmasqbar_betajOLSvalue)/((3*sigmasqbar_betajOLSvalue-TAUj_inf^2)*(sigmasqbar_betajOLSvalue+Bhattbar^2)+(Bhattbar^2-Bhattbar*TAUj_inf-sigmasqbar_beta0jOLSvalue)^2)
gOLStemp
}


hOLS<- function(wb1){
htempOLS <- 1000*(1-pnorm(((betaj11-betabarj0-tbarmin)-wb1[2])/(sigmasqbar_betajOLSvalue+wb1[1]^2)^0.5)+0.56*2*TAUj_inf/((27^0.5)*wb1[1]))/pnorm(wb1[2]/wb1[1])
htempOLS
}


TYPEII_BEfuncOLS <- optim(wb1start, hOLS, gr = NULL, 
      method = "L-BFGS-B", 
      lower = lowerBE, upper = c(10,10), control = list(fnscale=1)) 


TYPEIIOLS_C<- min(1,(sigmasqbar_betajOLSvalue)/(sigmasqbar_betajOLSvalue+(betaj11-betabarj0-tbarmin)^2))

TYPEIIOLS_Y <- min(1,gOLS(betaj11-betabarj0-tbarmin))

TYPEIIOLS_BE <- min(1,TYPEII_BEfuncOLS$value/1000)

TYPEIIOLS_H <- min(1,exp(-2*(betaj11-betabarj0-tbarmin)^2/TAUj_2))

if(betaj11<=betabarj0+tbarmin){
TYPEIIOLS_C <- 1
TYPEIIOLS_Y <- 1
TYPEIIOLS_BE <- 1
TYPEIIOLS_H <- 1 }



if (tbarmin == Cantellitbar) {
CutoffNonstandardizedOLS <- "C "
ModelNonstandardizedOLS <- "Non-Standardized test" }

if (min(TYPEIIOLS_C, TYPEIIOLS_Y, TYPEIIOLS_BE, TYPEIIOLS_H) == TYPEIIOLS_C) {
TYPEIINonstandardizedOLS <- TYPEIIOLS_C
TYPEIINonstandardizedOLSname <- "C " }

if (tbarmin == Bhattbar ) {
CutoffNonstandardizedOLS <- "Bh "
ModelNonstandardizedOLS <- "Non-Standardized test" }

if (min(TYPEIIOLS_C, TYPEIIOLS_Y, TYPEIIOLS_BE, TYPEIIOLS_H) == TYPEIIOLS_Y) {
TYPEIINonstandardizedOLS <- TYPEIIOLS_Y
TYPEIINonstandardizedOLSname <- "Bh "}

if (tbarmin == Hoeftbar ) {
CutoffNonstandardizedOLS <- "H "
ModelNonstandardizedOLS <- "Non-Standardized test" }

if (min(TYPEIIOLS_C, TYPEIIOLS_Y, TYPEIIOLS_BE, TYPEIIOLS_H) == TYPEIIOLS_H) {
TYPEIINonstandardizedOLS <- TYPEIIOLS_H
TYPEIINonstandardizedOLSname <- "H " }

if (tbarmin == BEtbar) {
CutoffNonstandardizedOLS <- "BE "
ModelNonstandardizedOLS <- "Non-Standardized test" }

if (min(TYPEIIOLS_C, TYPEIIOLS_Y, TYPEIIOLS_BE, TYPEIIOLS_H) == TYPEIIOLS_BE) {
TYPEIINonstandardizedOLS <- TYPEIIOLS_BE
TYPEIINonstandardizedOLSname <- "BE " }

Summary2OLS<- c("beta_j under alternative", "TYPE II bound Cantelli", "TYPE II bound Hoeffding", " TYPE II bound Bhattacharyya", "TYPE II bound Berry-Esseen", round(betaj1, digits=5), round(TYPEIIOLS_C, digits=5), round(TYPEIIOLS_H, digits=5), round(TYPEIIOLS_Y, digits=5), round(TYPEIIOLS_BE, digits=5))
dim(Summary2OLS)<- c(5,2)

mydata2OLS<- as.data.frame(Summary2OLS)
names(mydata2OLS)<- c("Type II related variables for OLS", "values")
print(mydata2OLS)


# -------------------- MM

obj <- c(rep(0,XROWS),1)
mat <- matrix(rbind(cbind(diag(XROWS), rep(-1,XROWS)), cbind(diag(XROWS), rep(1,XROWS)), cbind(t(X), rep(0,XCOLS)), cbind(t(rep(0,XROWS)), 1)), nrow=2*XROWS+XCOLS+1)
dir <- c(rep("<=", XROWS), rep(">=", XROWS), rep("==", XCOLS), ">=")
rhs <- c(rep(0, 2*XROWS), ej, 0)
XROWSplusone <- XROWS+1
sol <- Rglpk_solve_LP(obj, mat, dir, rhs, types = NULL, max = FALSE, bounds = list(lower=list(ind=1:XROWSplusone, val=rep(-XROWS,XROWSplusone)),
                          upper=list(ind=1:XROWSplusone, val=rep(XROWS,XROWSplusone))), verbose =FALSE)

Taumm_j<- sol$solution
Taumm_j<- Taumm_j[1:XROWS]
Betahatmmj=crossprod(Taumm_j, Y)
TAUjmm_2=sum((Taumm_j)^2)
TAUjmm_inf=max(abs(Taumm_j))


#Type 1 MM


# adjusting tau_j st Dmat positive definite, 
# rescale and then eliminate 0 entries
# increases bound so test remains exact

Taumm_jB <- rep(0,times=XROWS)
i <- 1
while (i <= XROWS) {
if (Taumm_j[i]>=0) 
Taumm_jB[i] <- max(Taumm_j[i]/min(1,TAUjmm_2),qqmm*min(1,TAUjmm_2)) 
else Taumm_jB[i] <- min(Taumm_j[i]/min(1,TAUjmm_2),-qqmm*min(1,TAUjmm_2)) 
i <- i+1
}

Dmatmm <- 2*t(X)%*%((Taumm_jB^2)*X)
dvecmm <- (1+2*ww)*colMeans((Taumm_jB^2)*X)*XROWS

if (det(Dmatmm)<=10^-15) {

print(paste("Dmat not positive definite for MM, Det=",det(Dmatmm)))
sigmasqbar_beta0jmmvalue <- TAUjmm_2/4-(1/XROWS)*(min(0,betabarj0-(1/2)*sum(Taumm_j)))^2

} else {

sigmasqbar_beta0jmmQP <- solve.QP(Dmat=Dmatmm, dvec=dvecmm, Amat, bvec=bvec0, meq=0, factorized=FALSE)

zsol0mm <- sigmasqbar_beta0jmmQP$solution 

sigmasqbar_beta0jmmvalue <- sum((Taumm_j)^2*((X%*%zsol0mm-ww)*(1+ww-X%*%zsol0mm)))
}

#
#
#
#
#
#
# Hoeffding and Cantelli MM
#
#
#
#
#
#


Hoeftbarmm <- (-0.5*TAUjmm_2*log(alpha))^0.5

Cantellitbarmm <- ((sigmasqbar_beta0jmmvalue)*(1-alpha)/alpha)^(0.5)


#
#
#
#
# Bhattacharyya MM
#
#
#
#
# 

if (Cantellitbarmm*(Cantellitbarmm - TAUjmm_inf)/sigmasqbar_beta0jmmvalue > 1) {

Bhattbarmm <- sqrt(sigmasqbar_beta0jmmvalue * (1 + sqrt( 3 *(1-alpha) / alpha ) ) )

if (Bhattbarmm^2 * TAUjmm_inf < sigmasqbar_beta0jmmvalue * ( TAUjmm_inf + 3 * Bhattbarmm)) {

r <- uniroot(function(x) ((3*sigmasqbar_beta0jmmvalue-TAUj_inf^2)*sigmasqbar_beta0jmmvalue)/((3*sigmasqbar_beta0jmmvalue-TAUj_inf^2)*(sigmasqbar_beta0jmmvalue+x^2)+(x^2-x*TAUj_inf-sigmasqbar_beta0jmmvalue)^2) - alpha, c(0,Bhattbarmm))

Bhattbarmm <- r$root
}

} else Bhattbarmm <- Inf

#
#
#
#
#
# Berry-Esseen MM
#
#
#
# 

hmm<- function(wb1){
hmmtemp=1000*(1-pnorm((tbartempmm-wb1[2])/(sigmasqbar_beta0jmmvalue+wb1[1]^2)^0.5)+0.56*2*TAUjmm_inf/((27^0.5)*wb1[1]))/pnorm(wb1[2]/wb1[1])
hmmtemp
}

BEtype1 <- optim(wb1start, h, gr = NULL, 
      method = "L-BFGS-B", 
      lower = lowerBE, upper = c(10,10), control = list(fnscale=1))

BE <- (BEtype1$value/1000 < alpha)

while(BE){
BEtype1 <- optim(wb1start, h, gr = NULL, 
      method = "L-BFGS-B", 
      lower = lowerBE, upper = c(10,10), control = list(fnscale=1))

if (BEtype1$value/1000 > alpha) (BE <- FALSE) 
else {
BEtbar <- tbartemp
print(BEtbar)
tbartemp <- tbartemp - grid
wb1start <- BEtype1$par
}
}



wb1mmstart <- c(0.1, 0.1)
tbartempmm <- Hoeftbarmm
BEtbarmm <- Inf

BEtype1mm <- optim(wb1mmstart, hmm, gr = NULL, 
      method = "L-BFGS-B", 
      lower = lowerBE, upper = c(10,10), control = list(fnscale=1))

BE <- (BEtype1mm$value/1000 < alpha)

while(BE){
BEtype1mm <- optim(wb1mmstart, hmm, gr = NULL, 
      method = "L-BFGS-B", 
      lower = lowerBE, upper = c(10,10), control = list(fnscale=1))

if (BEtype1mm$value/1000 > alpha) (BE <- FALSE) 
else {
BEtbarmm <- tbartempmm
# print(BEtbarmm)
tbartempmm <- tbartempmm - grid
wb1mmstart <- BEtype1$par
}
}

#
# 
#
#
#
# Final results for non-standardized test cutoff MM
#
#
#
#
#

tbarminmm <- min(c(Cantellitbarmm, Bhattbarmm, Hoeftbarmm, BEtbarmm))


# Type 2 MM


if (det(Dmatmm)<=10^-15) {

sigmasqbar_betajmmvalue <- TAUjmm_2/4-(1/XROWS)*(betaj11-(1/2)*sum(Taumm_j))^2

} else {

bvec1 <- as.vector(c(-betaj11, rep(ww, times=XROWS), rep(-(ww+1), times=XROWS)))

sigmasqbar_betajmmQP <- solve.QP(Dmat=Dmatmm, dvec=dvecmm, Amat, bvec=bvec1, meq=1, factorized=FALSE)

zsolmm <- sigmasqbar_betajmmQP$solution 

sigmasqbar_betajmmvalue <- sum((Taumm_j)^2*((X%*%zsolmm-ww)*(1+ww-X%*%zsolmm)))
}


gmm <- function(Bhattbar) {

if(((Bhattbar^2)/sigmasqbar_betajmmvalue)-Bhattbar*TAUjmm_inf/sigmasqbar_betajmmvalue - 1 <=0) gmmtemp <- 1 else
if(sigmasqbar_betajmmvalue < (Bhattbar^2)*TAUjmm_inf/(TAUjmm_inf+2*Bhattbar))
gmmtemp <- (3*sigmasqbar_betajmmvalue^2)/(4*(sigmasqbar_betajmmvalue^2)-2*sigmasqbar_betajmmvalue*Bhattbar^2+Bhattbar^4)
else gmmtemp <- ((3*sigmasqbar_betajmmvalue-TAUjmm_inf^2)*sigmasqbar_betajmmvalue)/((3*sigmasqbar_betajmmvalue-TAUjmm_inf^2)*(sigmasqbar_betajmmvalue+Bhattbar^2)+(Bhattbar^2-Bhattbar*TAUjmm_inf-sigmasqbar_beta0jmmvalue)^2)
gmmtemp
}

gOLS <- function(Bhattbar) {

if(((Bhattbar^2)/sigmasqbar_betajOLSvalue)-Bhattbar*TAUj_inf/sigmasqbar_betajOLSvalue - 1 <=0) gOLStemp <- 1 else
if(sigmasqbar_betajOLSvalue < (Bhattbar^2)*TAUj_inf/(TAUj_inf+3*Bhattbar))
gOLStemp <- (3*sigmasqbar_betajOLSvalue^2)/(4*(sigmasqbar_betajOLSvalue^2)-2*sigmasqbar_betajOLSvalue*Bhattbar^2+Bhattbar^4)
else gOLStemp <- ((3*sigmasqbar_betajOLSvalue-TAUj_inf^2)*sigmasqbar_betajOLSvalue)/((3*sigmasqbar_betajOLSvalue-TAUj_inf^2)*(sigmasqbar_betajOLSvalue+Bhattbar^2)+(Bhattbar^2-Bhattbar*TAUj_inf-sigmasqbar_beta0jOLSvalue)^2)
gOLStemp
}



hmm<- function(wb1){
htempmm<-1000*(1-pnorm(((betaj11-betabarj0-tbarminmm)-wb1[2])/(sigmasqbar_betajOLSvalue+wb1[1]^2)^0.5)+0.56*2*TAUj_inf/((27^0.5)*wb1[1]))/pnorm(wb1[2]/wb1[1])
htempmm
}


TYPEII_BEfuncmm <- optim(wb1mmstart, hmm, gr = NULL, 
      method = "L-BFGS-B", 
      lower = lowerBE, upper = c(10,10), control = list(fnscale=1)) 


TYPEIImm_C<- min(1,(sigmasqbar_betajmmvalue)/(sigmasqbar_betajmmvalue+(betaj11-betabarj0-tbarminmm)^2))

TYPEIImm_Y <- min(1,gmm(betaj11-betabarj0-tbarminmm))

TYPEIImm_BE <- min(1,TYPEII_BEfuncmm$value/1000)

TYPEIImm_H<- min(1,exp(-2*(betaj11-betabarj0-tbarminmm)^2/TAUjmm_2))

if(betaj11<=betabarj0+tbarminmm){
TYPEIImm_C <- 1
TYPEIImm_Y <- 1
TYPEIImm_BE <- 1
TYPEIImm_H <- 1}

if (tbarminmm == Cantellitbarmm) {
CutoffNonstandardizedmm <- "C"
ModelNonstandardizedmm <- "Non-Standardized test" }

if (min(TYPEIImm_C, TYPEIImm_Y, TYPEIImm_BE, TYPEIImm_H) == TYPEIImm_C) {
TYPEIINonstandardizedmm <- TYPEIImm_C
TYPEIINonstandardizedmmname <- "C"}

if (tbarminmm == Bhattbarmm ) {
CutoffNonstandardizedmm <- "Bh"
ModelNonstandardizedmm <- "Non-Standardized test" }

if (min(TYPEIImm_C, TYPEIImm_Y, TYPEIImm_BE, TYPEIImm_H) == TYPEIImm_Y) {
TYPEIINonstandardizedmm <- TYPEIImm_Y
TYPEIINonstandardizedmmname <- "Bh" }

if (tbarminmm == Hoeftbarmm ) {
CutoffNonstandardizedmm <- "H"
ModelNonstandardizedmm <- "Non-Standardized test" }

if (min(TYPEIImm_C, TYPEIImm_Y, TYPEIImm_BE, TYPEIImm_H) == TYPEIImm_H) {
TYPEIINonstandardizedmm <- TYPEIImm_H
TYPEIINonstandardizedmmname <- "H" }

if (tbarminmm == BEtbarmm) {
CutoffNonstandardizedmm <- "BE"
ModelNonstandardizedmm <- "Non-Standardized test" }

if (min(TYPEIImm_C, TYPEIImm_Y, TYPEIImm_BE, TYPEIImm_H) == TYPEIImm_BE) {
TYPEIINonstandardizedmm <- TYPEIImm_BE
TYPEIINonstandardizedmmname <- "BE" }

cat("\n")
cat("NON-STANDARDIZED TEST \n")

Summarymm1<- c("betabarj", "Betahatj", "Cantellitbar", "Bhattbar", "Hoeftbar", "BEtbar", round(betabarj, digits=5), round(IEc*Betahatmmj, digits=5), round(Cantellitbarmm, digits=5), round(Bhattbarmm, digits=5), round(Hoeftbarmm, digits=5), round(BEtbarmm, digits=5))
dim(Summarymm1)<- c(6,2)

mydata<- as.data.frame(Summarymm1)
names(mydata)<- c("statistic name (Minimax version)", "value")
print(mydata)
if(Betahatmmj-betabarj0>=tbarminmm) {
cat("Nonstandardized test >>> REJECTS H_0\n\n") 
RejectNonstandardizedmm <- "YES"} else {
cat("Non-standardized test does not reject H_0\n\n")
RejectNonstandardizedmm <- "NO" }

Summary2mm<- c("beta_j under alternative", "TYPE II bound Cantelli", "TYPE II bound Hoeffding", "TYPE II bound Bhattacharyya", "TYPE II bound Berry-Esseen", round(betaj1, digits=5), round(TYPEIImm_C, digits=5), round(TYPEIImm_H, digits=5), round(TYPEIImm_Y, digits=5), round(TYPEIImm_BE, digits=5))
dim(Summary2mm)<- c(5,2) 

mydata2mm<- as.data.frame(Summary2mm)
names(mydata2mm)<- c("Type II related variables for Minimax", "values")
print(mydata2mm)

#if(betaj11<=betabarj0+tbarminmm){ cat("Doesn't apply (MM).\n\n") }
cat("\n")


# ------------
#
#
#
#
#   Non-Randomised Bernoulli Test (OLS)
#
#
#
#
#
#
a <- 0
b <- XROWS*TAUj_inf
i <- 1
d <- rep(0,times=XROWS)
while(i <= XROWS){
d[i] <- (1-lambda) * max(-Tau_j[i]*ww,-Tau_j[i]*(ww+1)) + lambda * (TAUj_inf - max(Tau_j[i]*ww,Tau_j[i]*(ww+1))) 
i <- i+1
}
ds <- sum(d)

if (((betaj11+ds-a)/(b-a)>1)& (IE == "<=")) {
print(paste("   betaj1 is too large, has to be <= ", round(b-ds,digits=5)))
STOP
}
if (((betaj11+ds-a)/(b-a)>1)& (IE == ">=")) {
print(paste("   betaj1 is too small, has to be >= ", round(ds-b,digits=5)))
STOP
}


Z <- XROWS*(Tau_j*Y+d)
p1 <- (Z-a)/(b-a)

i <- 1
while(i <= XROWS){
if(p1[i]>1) {p1[i] <- 1} 
i <- i+1}

p0 <- 1-p1

W <- rep(0,times=monte*XROWS)
i <- 1
j <- 1
while(j<=monte){
for(i in 1:XROWS){
W[i+(j-1)*XROWS] <- rbern(1, p1[i])
} 
j <- j+1
}

Wbars<- rep(0,times=monte)
i<- 1
while(i<=monte){
temp1<- (i-1)*XROWS+1
temp2<- i*XROWS
Wbars[i]<- mean(W[temp1:temp2])
i<- i+1
}

pbar=(betabarj0+ds-a)/(b-a)

# Bkp is the probability of getting at least k successes
Bkp <- function(k,p){
Bkptemp <- pbinom(k-1,XROWS,p,lower.tail=FALSE)
Bkptemp
}

# we first look for smallest k st Bkp(k,pbar)<=alpha
k1 <- floor(pbar*XROWS)
while ((Bkp(k1,pbar)>alpha) & (k1 <= XROWS)) k1 <- k1+1
# so k1 is this smallest k

if (k1 <= XROWS) {

TYPEIIbernoulli <- function(betaj){
if (kbar/XROWS <= pbar) TYPEIIbernoullitemp <- 1 else 
TYPEIIbernoullitemp <- (1-Bkp(kbar,(betaj+ds-a)/(b-a)))/(1-theta)
TYPEIIbernoullitemp
}

TYPEII_B <- Inf
kb <- XROWS
theta <- Inf
kk <- k1           

while( Bkp(kk,pbar)/alpha > 0.05 ){

alphabar <- Bkp(kk,pbar)

theta <- alphabar/alpha

kbar <- kk
# print(theta)
# print(TYPEIIbernoulli(betaj11))

if (TYPEIIbernoulli(betaj11) < TYPEII_B) {

TYPEII_B <- TYPEIIbernoulli(betaj11)

kb <- kbar #value of kbar to remember for below

}

kk <- kk+1
}

kbar <- kb

alphabar <- Bkp(kbar,pbar)

theta <- alphabar/alpha


r_alphaprimeWbar <- function(Wbar){
if(XROWS*Wbar >= kbar){
r_alphaprimeWbartemp <- 1 
} else if(XROWS*Wbar == kbar-1){
r_alphaprimeWbartemp <- (alphabar-Bkp(kbar,pbar))/(Bkp(kbar-1,pbar)-Bkp(kbar,pbar))
} else {
r_alphaprimeWbartemp <- 0
}
r_alphaprimeWbartemp
}

r_alphaprimeWbars <- rep(0,times=monte)
i<-1
while(i<=monte){
r_alphaprimeWbars[i]<- r_alphaprimeWbar(Wbars[i])
i<- i+1
}


cat("NON-RANDOMISED BERNOULLI TEST (using OLS estimator) \n")
if (lambda < 1) cat(paste("lambda =", lambda," \n"))

if (mean(r_alphaprimeWbars)>=theta){
cat(">>> REJECTS H0  (rejection probability of randomized test =", mean(r_alphaprimeWbars),")")
RejectBernoulliOLS <- "YES"
} else {
cat("does not reject H0, rejection probability of randomized test =",mean(r_alphaprimeWbars))
RejectBernoulliOLS <- "NO"
}
cat("\n")


SummaryBernoulli<- c("theta", "TYPE II bound Bernoulli", round(theta, digits=5), round(min(1,TYPEII_B), digits=5))
dim(SummaryBernoulli)<- c(2,2)
 

mydataBernoulli<- as.data.frame(SummaryBernoulli)
names(mydataBernoulli)<- c("Type II related variables", "values")
print(mydataBernoulli)

} else {

TYPEII_B <- Inf

cat("NON-RANDOMISED BERNOULLI TEST (using OLS estimator) \n")
if (lambda < 1) cat(paste("lambda =", lambda," \n"))
cat(" .... need to reduce lambda to make Bernoulli applicable under OLS")
cat("\n")
}

#
#
#
#
#
#
#   Non-Randomised Bernoulli Test (MM)
#
#
#
#
#
#

a <- 0
bmm <- XROWS*TAUjmm_inf
i <- 1
dmm <- rep(0,times=XROWS)
while(i <= XROWS){
dmm[i] <- (1-lambdamm) * max(-Taumm_j[i]*ww,-Taumm_j[i]*(ww+1)) + lambdamm * (TAUjmm_inf - max(Taumm_j[i]*ww,Taumm_j[i]*(ww+1))) 
i <- i+1
}
dsmm <- sum(dmm)

if (((betaj11+dsmm-a)/(b-a)>1)&(IE == "<=")) {
print(paste("   betaj1 is too large, has to be <= ", round(b-dsmm,digits=5)))
STOP
}
if (((betaj11+dsmm-a)/(b-a)>1)&(IE == ">=")) {
print(paste("   betaj1 is too small, has to be >= ", round(dsmm-b,digits=5)))
STOP
}



Zmm <- XROWS*(Taumm_j*Y+dmm)
p1mm <- (Zmm-a)/(bmm-a)

i <- 1
while(i <= XROWS){
if(p1mm[i]>1) {p1mm[i] <- 1} 
i <- i+1}

p0mm <- 1-p1mm

Wmm<- rep(0,times=monte*XROWS)
i<- 1
j<- 1
while(j<=monte){
for(i in 1:XROWS){
Wmm[i+(j-1)*XROWS] <- rbern(1, p1mm[i])
} 
j<- j+1
}

Wbarsmm<-rep(0,times=monte)
i<- 1
while(i<=monte){
temp1mm<- (i-1)*XROWS+1
temp2mm<- i*XROWS
Wbarsmm[i]<- mean(Wmm[temp1mm:temp2mm])
i<- i+1
}

pbarmm <- (betabarj0+dsmm-a)/(bmm-a)

k1mm <- floor(pbarmm*XROWS)
while ((Bkp(k1mm,pbarmm)>alpha) & (k1mm <= XROWS)) k1mm <- k1mm+1

if (k1mm <= XROWS) {

TYPEIIbernoullimm<-  function(betaj){
if (kbarmm/XROWS <= pbarmm) TYPEIIbernoullitempmm <- 1 else 
TYPEIIbernoullitempmm <- (1-Bkp(kbarmm,(betaj+dsmm-a)/(bmm-a)))/(1-thetamm)
TYPEIIbernoullitempmm
}

TYPEII_Bmm <- Inf
kbmm <- XROWS
thetamm <- Inf
kkmm <- k1mm           

while( Bkp(kkmm,pbarmm)/alpha > 0.05 ){

alphabarmm <- Bkp(kkmm,pbarmm)

thetamm <- alphabarmm/alpha

kbarmm <- kkmm

# print(thetamm)
# print(TYPEIIbernoullimm(betaj11))

if (TYPEIIbernoullimm(betaj11) < TYPEII_Bmm) {

TYPEII_Bmm <- TYPEIIbernoullimm(betaj11)

kbmm <- kbarmm

}

kkmm <- kkmm+1
}

kbarmm <- kbmm

alphabarmm <- Bkp(kbarmm,pbarmm)

thetamm <- alphabarmm/alpha


r_alphaprimeWbarmm <- function(Wbar){
if(XROWS*Wbar >= kbarmm){
r_alphaprimeWbartempmm <- 1 
} else if(XROWS*Wbar == kbarmm-1){
r_alphaprimeWbartempmm <- (alphabarmm-Bkp(kbarmm,pbarmm))/(Bkp(kbarmm-1,pbarmm)-Bkp(kbarmm,pbarmm))
} else {
r_alphaprimeWbartempmm <- 0
}
r_alphaprimeWbartempmm
}

r_alphaprimeWbarsmm<-rep(0,times=monte)
i<- 1
while(i<=monte){
r_alphaprimeWbarsmm[i]<- r_alphaprimeWbarmm(Wbarsmm[i])
i<- i+1
}

cat("\n")
cat("NON-RANDOMISED BERNOULLI TEST (using minimax estimator)\n")
if (lambdamm < 1) cat(paste("lambdamm =", lambdamm," \n"))

if (mean(r_alphaprimeWbarsmm)>= thetamm){
cat(">>> REJECTS H0  (rejection probability of randomized test =",mean(r_alphaprimeWbarsmm),")")
RejectBernoullimm <- "YES"
} else {
cat("does not reject H0, rejection probability of randomized test =",mean(r_alphaprimeWbarsmm))
RejectBernoullimm <- "NO"}
cat("\n")

SummaryBernoullimm<- c("theta", "TYPE II bound Bernoulli", round(thetamm, digits=5), round(min(1,TYPEII_Bmm), digits=5))
dim(SummaryBernoullimm)<- c(2,2)
 

mydataBernoullimm<- as.data.frame(SummaryBernoullimm)
names(mydataBernoullimm)<- c("Type II related variables", "values")
print(mydataBernoullimm)

} else {

TYPEII_Bmm <- Inf

cat("\n")
cat("NON-RANDOMISED BERNOULLI TEST (using minimax estimator)\n")
if (lambdamm < 1) cat(paste("lambdamm =", lambdamm," \n"))
cat(" .... need to reduce lambdamm to make Bernoulli applicable under MM")
cat("\n")
}

##################################################USER OUTPUT

# changed: Error output

cat("\n")

if (det(DmatOLS)<=10^-15) {
print("Upper bound on sigma^2bar used as Dmat not positive definite (OLS)")
}
if (det(Dmatmm)<=10^-15) {
print("Upper bound on sigma^2bar used as Dmat not positive definite (MM)")
}
if (k1 > XROWS){
print(" ... need to decrease lambda to apply Bernoulli under OLS")
}
if (k1mm > XROWS){
print(" ... need to decrease lambdamm to apply Bernoulli under MM")
}



#Selecting the test
cat("\n")
if (min(TYPEIINonstandardizedOLS, TYPEIINonstandardizedmm, TYPEII_B, TYPEII_Bmm) == 1) {
if (IE == "<=") cat("betaj1 needs to be made larger") else cat("betaj1 needs to be made smaller")
cat("\n") } else {
if (min(TYPEIINonstandardizedOLS, TYPEIINonstandardizedmm, TYPEII_B, TYPEII_Bmm) == TYPEIINonstandardizedOLS) {

UserSummary<- c(paste("Estimated coefficient of beta", cj, ":"), paste("Reject H0: beta", cj,IE, betabarj,"?"), "based on:", "Estimator (tau):", paste("using",CutoffNonstandardizedOLS, "cutoff (t-bar):"), paste("Type II bound (", TYPEIINonstandardizedOLSname,") for beta",cj, IEa, round(betaj1,digits=5),":"), round(IEc*Betahatj, digits=5), paste(RejectNonstandardizedOLS, " at alpha = ", alpha), ModelNonstandardizedOLS, "OLS", round(tbarmin, digits=5), round(TYPEIINonstandardizedOLS, digits=5))
dim(UserSummary)<- c(6,2)
 
mydataUserSummary<- as.data.frame(UserSummary)
names(mydataUserSummary)<- c(paste("Summary Results, n=",XROWS,"m=",XCOLS), paste(" Y range [",w1Y,",",w2Y,"]"))
print(mydataUserSummary) }


if (min(TYPEIINonstandardizedOLS, TYPEIINonstandardizedmm, TYPEII_B, TYPEII_Bmm) == TYPEIINonstandardizedmm) {

UserSummary<- c(paste("Estimated coefficient of beta", cj, ":"), paste("Reject H0: beta", cj,IE, betabarj,"?"), "based on:", "Estimator (tau):", paste("using",CutoffNonstandardizedmm, "cutoff (t-bar):"), paste("Type II bound (", TYPEIINonstandardizedmmname,") for beta",cj, IEa, round(betaj1,digits=5),":"), round(IEc*Betahatj, digits=5), paste(RejectNonstandardizedmm, " at alpha = ", alpha),  ModelNonstandardizedmm, "Minimax", round(tbarmin, digits=5), round(TYPEIINonstandardizedmm, digits=5))

dim(UserSummary)<- c(6,2)
 

mydataUserSummary<- as.data.frame(UserSummary)
names(mydataUserSummary)<- c(paste("Summary Results, n=",XROWS,"m=",XCOLS), paste(" Y range [",w1Y,",",w2Y,"]"))
print(mydataUserSummary) }

if (min(TYPEIINonstandardizedOLS, TYPEIINonstandardizedmm, TYPEII_B, TYPEII_Bmm) == TYPEII_B) {

UserSummary<- c(paste("Estimated coefficient of beta", cj, ":"), paste("Reject H0: beta", cj, IE, betabarj,"?"), "based on:", "Estimator:", "theta:", paste("Type II bound for beta", cj, IEa, round(betaj1,digits=5),":"), round(IEc*Betahatj, digits=5), paste(RejectBernoulliOLS, " at alpha = ", alpha), "Bernoulli" , "OLS", round(theta, digits=5), round(TYPEII_B, digits=5))
dim(UserSummary)<- c(6,2)
 

mydataUserSummary<- as.data.frame(UserSummary)
names(mydataUserSummary)<- c(paste("Summary Results, n=",XROWS,"m=",XCOLS), paste(" Y range [",w1Y,",",w2Y,"]"))
print(mydataUserSummary) }

if (min(TYPEIINonstandardizedOLS, TYPEIINonstandardizedmm, TYPEII_B, TYPEII_Bmm) == TYPEII_Bmm) {

UserSummary<- c(paste("Estimated coefficient of beta", cj, ":"), paste("Reject H0: beta", cj, IE, betabarj,"?"), "based on:", "Estimator (tau):", "theta:", paste("Type II bound for beta", cj, IEa, round(betaj1,digits=5),":"), round(IEc*Betahatj, digits=5), paste(RejectBernoullimm, " at alpha = ", alpha), "Bernoulli" , "Minimax", round(thetamm,digits=5), round(TYPEII_Bmm, digits=5))
dim(UserSummary)<- c(6,2)
 

mydataUserSummary<- as.data.frame(UserSummary)
names(mydataUserSummary)<- c(paste("Summary Results, n=",XROWS,"m=",XCOLS), paste(" Y range [",w1Y,",",w2Y,"]"))
print(mydataUserSummary) }}