program npuconf, rclass
    version 9.1
    args varname
    syntax varname [if] [in] [,low (real 0) up(real 1) times(integer 10000) theta(real 0.2) alpha(real 0.025) w(real 0.5)]
    marksample touse

    local pseudoalpha = `alpha'*`theta' 
    local v end
    quietly{
        summarize `varname'
    }
    local m=r(N)
    local mean1=r(mean)



        tempvar y
        gen `y'= (`varname' - `low')/(`up' - `low')
        local p = (`w' - `low')/(`up' - `low')
        

    mata:
        P=J(`m',1,`p')
        st_view(Y=., ., "`y'")
        A=Y-P
        rj=0 /* rejection probability under size pseudoalpha */

        
        q=uniform(`m',1)

        for (t=1; t<=`times'; t++)
        {
            q=uniform(`m',1)
            s1=0
            s2=0
            for (i=1; i<=`m' ;i++) 
            {
                if (A[i,1]>q[i,1]*(1-`p'))
                {
                    s2++
                }
                else if (Y[i,1]<q[i,1]*`p')
                {
                    s1++
                }
            }
    
            
                    
    
            h1=0
    
            for (k=0; k<=s2; k++)
            {
                h1 = h1 + (`p'^k)*((1-`p')^(s1+s2-k))*comb(s1+s2,k)
            }    
            if (h1 <= `pseudoalpha')
            {
                rj=rj+(1/`times')
            }
            else
            {
                h = (`p'^s2)*((1-`p')^s1)*comb(s1+s2,s2)
                if (h1 <= `pseudoalpha' + h)
                {
                    rj= rj +((`pseudoalpha' - h1 + h)/h)/`times'
                }
            }
        }
    st_numscalar("r(rj)", rj)
    `v'
    local rj = r(rj)
	local rj1=round(`rj',.001)
	local rj2=round(`rj',.000001)



	local mean2=round(`mean1',.001)
    di as text _col(15) "Nonparametric test for an upper confidence bound (Schlag, 2008)"
    di as text _dup(83) "-"
di in smcl in gr " Variable  "_col(13) "{c |}    Obs" _col(26) "low" /*
        */ _col(35) "up" _col(42) "theta" _col(50) " alpha " _col(61) "w" _col(68) "avg     rj " _n /*
        */ "{hline 12}{c +}{hline 70}"
    di in smcl   " `varname'  " _col(13) "{c |}"in ye"    `m'" _col(27)  "`low'" /*
        */ _col(36) "`up'" _col(43) "`theta'" _col(51) " `alpha' " _col(61) "`w'" _col(70) "`mean2'     `rj1' "
        */"{hline 12}{c +}{hline 70}"

return scalar rej=`rj'
return scalar n=`m'
return scalar lb=`low'
return scalar ub=`up'
return scalar alpha=`alpha'
return scalar theta=`theta'
return scalar psalpha=`pseudoalpha'

    
    if (`rj' >= `theta'){
        di _newline /*
        */ in gr "HO: E(`varname')>=`w' rejected "
        }

    else{
        di _newline /*
        */ in gr "HO: E(`varname')>=`w' not rejected "
    }

 
end