####################### JEL and AJEL ratio tests ######################
rm(list=ls())
library(emplik)
MC<-10000     # no. of repition
type15<-rep()
typea15<-rep()
sam=c(20 ,40, 60, 80, 100) 

for(d in 1:5){
  n<-sam[d]
  Jkn<-rep()
  adJkn<-rep()
  k1= 1/choose(n,3)
  k11=1/choose((n-1),3)
  
  for(r in 1:MC){
    #Alternative distributions
        x<-rt(n,3)               
      # x<- rnorm(n,0,1)
      # x<-rgamma(n,2,1)
      # x<-rlaplace(n,0,1)
      # x<-rbeta(n, 2, 2)
      # x<- runif(n,0,1)
        
    xs= sort(x)
    s1<-0
    for(i in 1:(n-2)){
      for(j in (i+1):(n-1)){
        for(k in (j+1):(n )){
          s1<-s1+1*(0.5*((xs[i]*xs[j]-1)/xs[j]) < xs[k])
        }
      }
    }
    delta1= k1*s1
    s1<-0
    for(i in 1:(n-2)){
      for(j in (i+1):(n-1)){
        for(k in (j+1):(n )){
          s1<-s1+1*(0.5*((xs[i]*xs[k]-1)/ xs[k]) < xs[j])
        }
      }
    }
    delta2=  k1*s1
    s1<-0
    for(i in 1:(n-2)){
      for(j in (i+1):(n-1)){
        for(k in (j+1):(n )){
          s1<-s1+1*(0.5*((xs[j]*xs[k]-1)/ xs[k]) < xs[i])
        }
      }
    }
    delta3= k1*s1
    t<- ((1/3)*(delta1+delta2+delta3))-0.5
    
    ########## 
    
    v<-rep()  
    tm<-rep()  
    for(m in 1:n){
      x1<-x[-m]
      xs1= sort(x1)
      s11<-0
      for(i in 1:(n-3)){
        for(j in (i+1):(n-2)){
          for(k in (j+1):(n-1)){
            
            s11<-s11+1*(0.5*((xs1[i]*xs1[j]-1)/xs1[j]) < xs1[k])
            
          }
        }
      }
      delta11=k11*s11
      s11<-0
      for(i in 1:(n-3)){
        for(j in (i+1):(n-2)){
          for(k in (j+1):(n-1)){
            s11<-s11+1*(0.5*((xs1[i]*xs1[k]-1)/xs1[k]) <  xs1[j])
            
          }
        }
      }
      delta12=k11*s11
      s11<-0
      for(i in 1:(n-3)){
        for(j in (i+1):(n-2)){
          for(k in (j+1):(n-1)){
            
            s11<-s11+1*(0.5*((xs1[j]*xs1[k]-1)/xs1[k]) <  xs1[i])
            
          }
        }
      }
      delta13=k11*s11
      
      tm[m]<- ((1/3)*(delta11+delta12+delta13))-0.5
      v[m]= (n*t) - (n-1)*tm[m]                   ##### Pseduo Value   ######
    }
    
    ########### JEL and AJEL ******##
    
    av1<- -max(1,log(n,exp(1))/2)*mean(v)
    av<-c(v,av1)
    
    Jkn[r]<- el.test(v,mu=0)$'-2LLR'	   # Log likelihood ratio value obtained for JEL
    adJkn[r]<- el.test(av,mu=0)$'-2LLR'	   # Log likelihood ratio value obtained for Ad-JEL

  }
  type15[d]  <-mean(1*(abs(Jkn[Jkn!='NaN'])>3.84))
  typea15[d] <-mean(1*(abs(adJkn[adJkn!='NaN'])>3.84))
}
type15 
typea15