#------------------------------------------------------------------------#
# modified on   :   Tuesday, November 3, 2009 at 14:29
# created on    :   Thursday, July 24, 2008 at 16:25
# author        :   Huixia Wang
# AFFILIATION   :   DEPARTMENT OF STATISTICS, NORTH CAROLINA STATE UNIVERSITY
# EMAIL         :   WANG@STAT.NCSU.EDU
#------------------------------------------------------------------------#

library(quantreg)
library(MASS)




###########################################################################################################
#  ELLDrq:        computes the empirical likelihood ratio (-2loglr) for longitudinal data in quantreg
#                  the lagrange multipler is solved using Wu's modified Newton Raphson's method
#                assuming working independence, block smoothed EL
#       For faster computation, assume balanced design for the simulated data set
#  Input:                                  
#  x          N x p design matrix, N=n*m in the simulation
#########  subject      index of subjects
#  n          number of subjects
#  m           number of repeats for each subject
#  y          N x 1 responses            
#  beta       hypothesized regressor value p*1
#  tau        quantile of interest
#  h          smoothing parameter used to approximate the indicator function I(u<0) \approx 1-Phi(u/h)
#  maxit      maximum iterations                        
#  itertrace  1 to trace iterations (default = 0)       
#  Output:                            
#  logelr     log empirical likelihood ratio
###########################################################################################################

Ku = function(u)
{
    # kernel
    1*(abs(u)<=sqrt(5))*3/(4*sqrt(5))*(1-u^2/5)
}
Gu = function(u)
{
    # CDF (based on kernel Ku)
    u = pmin(u, sqrt(5))
    3/(4*sqrt(5))*(u - u^3/15 + 2*sqrt(5)/3)*(abs(u)<=sqrt(5))
}
psi.h = function(u, hn, tau)
{   
    # the smoothed quantile score function
    # unsmoothed when hn=0, out = I(u<0) - tau
    if(hn==0) out = 1*(u<0) - tau
    if(hn>0) out = 1-Gu(u/hn) - tau
    return(out)
}
rho.h = function(u, hn, tau)
{
    # smoothed check function
    -u*(psi.h(u,hn,tau))
}



############# bartlet correction factor calculation ############
bartlett = function(x, y, n, m, tau=0.5, hn)
{
    # for now, designed for simulated data with balanced design
    # bhat is any consistent estimator of beta such as rq or EL estimate
    # x includes the intercept
    beta = rq(y~x-1, tau=tau)$coef

    x = as.matrix(x)
    N = n*m
    p = ncol(x)

    ### z is n by p matrix!
    if(p==1)
    {
        X = matrix(x, ncol=n)   # m * n matrix
        E = matrix(psi.h(y-x*beta, hn, tau), ncol=n)  # m*n
        z = apply(X*E, 2, sum)
    }
    if(p>1)
    {
        E = psi.h(y-x%*%beta, hn, tau)
        XE = array(cbind(x,E), dim=c(m, n, p+1))
        z = t(apply(XE, 2, function(x) apply(x[,1:p] * x[,p+1], 2, sum)))        
    }
    Vn = t(z)%*%z/n
    alpha.iikk = 1/n * sum((diag(z %*% solve(Vn) %*% t(z)))^2)
    alpha.ikm2 = 1/(n^2) * sum(( z %*% solve(Vn) %*% t(z))^3)
    bart.b = p^(-1) *( alpha.iikk/2 - alpha.ikm2/3)
    return(bart.b)
}


ELLDrq = function(beta, x, y, n, m, tau=0.5, hn=0.1)
{
########## for balanced design in the simulation study only!!!
    # beta=c(1,1); tau=0.5; hn=0.1
    maxit=50; itertrace=0
    x = as.matrix(x)
    N = n*m
    p = ncol(x)
    if(is.na(beta)) beta = as.vector(rep(0,p))
    if(length(beta) != p)  warning("Need beta to be the same dimension with the number of covariates")

    ### z is n by p matrix!
    if(p==1)
    {
        X = matrix(x, ncol=n)   # m * n matrix
        E = matrix(psi.h(y-x*beta, hn, tau), ncol=n)  # m*n
        z = apply(X*E, 2, sum)
    }
    if(p>1)
    {
        E = psi.h(y-x%*%beta, hn, tau)
        XE = array(cbind(x,E), dim=c(m, n, p+1))
        z = t(apply(XE, 2, function(x) apply(x[,1:p] * x[,p+1], 2, sum)))        
    }
    # solve for lambda and pi
    loglr = Lag.Wu(z, maxit, itertrace)
    return(loglr)
}

Lag.Wu = function(z, maxit=50, itertrace=0)
{
# using the modified Newton-Raphson's method of Wu to
# solve for Lagrange multiplier lam such that g(lam) = sum_{i=1}^n z_i(beta) /[1 + lambda^T z_i(beta) ] = 0
# INPUTS
# z         : a n*p matrix, score function, depends on beta
# itertrace : =1/0 track or don't track the iteration process
# maxit     : maximum iteration
#
# OUTPUTS
# loglr     : =-2logEL ratio = -2 \sum_{i=1}^n log(n \hat p_i) (So want to minimize loglr)

# maxit=50; itertrace=1
    z = as.matrix(z)
    n = nrow(z)
    p = ncol(z)
    if(p==1){
        L<--1/max(z)
        R<--1/min(z)
        gsize <-1
        tol <-1e-8
        k = 0
        if(itertrace) cat("lam, logl, loglr, gsize, k", "\n")
        while(gsize>tol  & k<=maxit){
            lam = (L+R)/2
            glam = sum(z/(1+lam*z))   # decreases with lam
            if(glam>0) L = lam
            if(glam<0) R = lam  # glam<0 so need increase glam, decrease lam, so reduce the right end of lam: R
            gsize = abs(glam)
            k = k + 1
            phat <- 1/as.vector((1+z*lam)*n)
            logl = sum(log(phat)) 
            loglr = -2 * sum(log(n*phat)) 
            glam = apply(z/as.vector((1+z*lam)),2,sum)
            gsize = mean(abs(glam))
            if(itertrace) 
                print(c(lam, logl, loglr, gsize, k))
        }
    }
    # check if 1 + lam^Tz_i > 1/n for all i
    #arg = 1 + z%*%lam
    #sum(arg < 0)
    
    # for higher dimensional z
    if(p>1){
        lam = rep(0, p)
        gsize = 1
        tol = 1e-08
        D1 = rep(0, p)
        D2 = matrix(0, p, p)
        eval = 1
        k = 0  # number of iteration
        if(itertrace) cat("lam, logl, loglr, eval=|hessian|, gsize=|g(lam)|, iter", "\n")
        while(eval > tol & k<=maxit){
            arg <- as.vector(1 + z%*%lam)  #n*1 vector
            D1 = apply(z / arg, 2, mean)
            tmp = t(apply(z, 1, function(x) x%*%t(x)))  # dim=n * (p^2)
            tmp = tmp/(arg^2)
            DD <- - matrix(apply(tmp, 2, mean), ncol=p)
           # Delta <- solve(D2, D1, tol=1e-15)
            eig <- eigen(DD)$values
            if(max(abs(eig))/min(abs(eig))>=1e+07){
            lam <- rep(0,p)
            k <- maxit + 1
                }
            if(max(abs(eig))/min(abs(eig))<1e+07){
                #D2<-solve(DD, D1)
                D2 = solve(DD)%*% D1
                eval <- max(abs(D2))
                delta <- D2
                dd <- min(1+(z%*%(lam-delta)))  #check if all (1+z lambda) > 1/n for all i  
                while(dd<1/n){
                    delta <- delta/2
                    dd <- min(1+(z%*%(lam-delta)))
                }
                lam <-lam-delta
            }
            k <- k+1
            phat <- 1/as.vector((1+z%*%lam)*n)
            logl = sum(log(phat)) 
            loglr = -2 * sum(log(n*phat)) 
            glam = apply(z/as.vector((1+z%*%lam)),2,sum)
            gsize = mean(abs(glam))
            if(itertrace) 
                print(c(lam, logl, loglr, eval, gsize, k))
        }
    }
    if(k>=maxit) loglr<-1000  # diverge, or for such beta, g(lam) has no solution
    return(loglr)  
}

            ###########################################################################
            ### EL confidence interval (longitudinal)                                 #
            ###  obtain the CI for each compoenent of the pars separately             #
            ###########################################################################
CI.ELLDrq = function(x, y, n, m, tau=0.5, hn=0.1, alpha=0.05)
{
# tau=0.5; hn=0.1; maxit=50; itertrace=0
    x = as.matrix(x)
    N = length(y)
    p = ncol(x)
    if(p==1) warning("see function-ELLDRQ-v081208.R for the corresponding function")
    
    # calculate the bartlett correction factor
    bart.b = bartlett(x, y, n, m, tau, hn)
    
    if(p>1)
    {
        Cut = qchisq(1-alpha, 1)
        Cut2 = Cut*(1+ bart.b/n)   # bartlet corrected cutoff
        
        # LOG LIKELIHOOD UNDER THE FULL MODEL
        rq1 = summary(rq(y~x-1, tau=tau), se="iid")$coef
        EL.F = optim(par=rq1[,1], fn=ELLDrq, x=x, y=y, n=n, m=m, tau=tau, hn=hn, method="Nelder-Mead")
        beta.full = EL.F$par
        logel.full = EL.F$value
        coeff = coeff2 = matrix(0, nrow=p, ncol=3)
        for(j in 1:p)
        {
            b0 = beta.full[j]
            se = rq1[j, 2]
            # search for the right end of the interval
            t1 = b0
            t2 = b0 + 5*se
            dif1 = t2 - t1
            tol = 1e-08
            while(dif1>tol)
            {
                tt = (t1+t2)/2
                ystar = y - x[,j] * tt
                xstar = x[,-j]
                # obtain the log likelihood under the restricted model under H0: beta(j) = tt
                beta.R = rep(0, p)
                beta.R[-j] = optim(par=beta.full[-j], fn=ELLDrq, x=xstar, y=ystar, n=n, m=m, tau=tau, hn=hn, method="Nelder-Mead")$par
                beta.R[j] = tt
                logel.R = ELLDrq(beta.R, x, y, n=n, m=m, tau=tau, hn=hn)
                loglr = -logel.full + logel.R
                if(loglr > Cut) t2 = tt
                if(loglr < Cut) t1 = tt
                dif1 = t2-t1
               # print(c(t1, t2))
            }
            coeff[j, 1] = b0
            coeff[j, 3] = (t1+t2)/2
            
            # search for the left end of the interval
            t1 = b0
            t2 = b0 - 5*se
            dif1 = t1 - t2
            tol = 1e-08
            while(dif1>tol)
            {
                tt = (t1+t2)/2
                ystar = y - x[,j] * tt
                xstar = x[,-j]
                # obtain the log likelihood under the restricted model under H0: beta(j) = tt
                beta.R = rep(0, p)
                beta.R[-j] = optim(par=beta.full[-j], fn=ELLDrq, x=xstar, y=ystar, n=n, m=m, tau=tau, hn=hn, method="Nelder-Mead")$par
                beta.R[j] = tt
                logel.R = ELLDrq(beta.R, x, y, n=n, m=m, tau=tau, hn=hn)
                loglr = -logel.full + logel.R
                if(loglr > Cut) t2 = tt
                if(loglr < Cut) t1 = tt
                dif1 = t1-t2
                # print(c(t1, t2))
            }
            coeff[j, 2] = (t1+t2)/2
            
            #### BARTLET CORRECTED CI, start with the CI of SEL
            # search for the right end of the interval
            se = (coeff[j,3]- coeff[j,2])/4
            t1 = coeff[j, 3]
            t2 = t1 + se
            dif1 = t2 - t1
            while(dif1>tol)
            {
                tt = (t1+t2)/2
                ystar = y - x[,j] * tt
                xstar = x[,-j]
                # obtain the log likelihood under the restricted model under H0: beta(j) = tt
                beta.R = rep(0, p)
                beta.R[-j] = optim(par=beta.full[-j], fn=ELLDrq, x=xstar, y=ystar, n=n, m=m, tau=tau, hn=hn, method="Nelder-Mead")$par
                beta.R[j] = tt
                logel.R = ELLDrq(beta.R, x, y, n=n, m=m, tau=tau, hn=hn)
                loglr = -logel.full + logel.R
                if(loglr > Cut2) t2 = tt
                if(loglr < Cut2) t1 = tt
                dif1 = t2-t1
               # print(c(t1, t2))
            }
            coeff2[j, 1] = b0
            coeff2[j, 3] = (t1+t2)/2
            # search for the left end of the interval
            t1 = coeff[j, 2]
            t2 = t1 - se
            dif1 = t1 - t2
            tol = 1e-08
            while(dif1>tol)
            {
                tt = (t1+t2)/2
                ystar = y - x[,j] * tt
                xstar = x[,-j]
                # obtain the log likelihood under the restricted model under H0: beta(j) = tt
                beta.R = rep(0, p)
                beta.R[-j] = optim(par=beta.full[-j], fn=ELLDrq, x=xstar, y=ystar, n=n, m=m, tau=tau, hn=hn, method="Nelder-Mead")$par
                beta.R[j] = tt
                logel.R = ELLDrq(beta.R, x, y, n=n, m=m, tau=tau, hn=hn)
                loglr = -logel.full + logel.R
                if(loglr > Cut2) t2 = tt
                if(loglr < Cut2) t1 = tt
                dif1 = t1-t2
                # print(c(t1, t2))
            }
            coeff2[j, 2] = (t1+t2)/2
        }
    }
    colnames(coeff) = colnames(coeff2) = c("coef EL","lower bd","upper bd")
    return(list(CI.SEL = coeff, CI.SELc = coeff2))
}

                    #==========================================================#
                    # EL hypothesis testing, H0: beta(j) = theta (logitudinal) #
                    #       test the j-th component                            #
                    #==========================================================#
test.ELLDrq = function(x, y, n, m, comp, theta, tau=0.5, hn=0.1)
{
    # comp: the index of the coefficient compoenent to be tested, comp=1 means test on intercept
    # for balanced design only!
    x = as.matrix(x)
    p = ncol(x)
    if(p<=1) warning("For univariate case, this can be done by simply calculating loglr with function ELLDrq")
    j = comp
    q = length(comp)
    
    # LOG LIKELIHOOD UNDER THE FULL MODEL
        beta0 = rq(y~x-1, tau=tau)$coef
        EL.F = optim(par=beta0, fn=ELLDrq, x=x, y=y, n=n, m=m, tau=tau, hn=hn, method="Nelder-Mead")
        beta.full = EL.F$par
        logel.full = EL.F$value
        
    # LOG LIKELIHOOD UNDER THE RESTRICTED MODEL
        xstar = x[,-j]
        if(q>1) ystar = y - x[,j]%*%theta
        if(q==1) ystar = y - x[,j]*theta
        beta.R = rep(0, p)
        beta.R[-j] = optim(par=beta.full[-j], fn=ELLDrq, x=xstar, y=ystar, n=n, m=m, tau=tau, hn=hn, method="Nelder-Mead")$par
        beta.R[j] = theta
        logel.R = ELLDrq(beta.R, x, y, n=n, m=m, tau=tau, hn=hn)
        loglr = -logel.full + logel.R
    pval = 1-pchisq(loglr,q)
    
    out = c(beta0[j], loglr, pval)
    names(out) = c(rep("coefEL",q), "logelr", "pval")
    return(out)       
}



################################### Wald type test
Wald = function(y, x, subject, n, m, h, tau, alpha=0.05)
{
# works also for unbalanced design
    zstar = as.matrix(x )  
    # estimation under full model
    rq1 = rq(y~x-1, tau=tau)
    beta = rq1$coeff
    ehat = rq1$resid
    ehat = ehat - quantile(ehat, tau)      ########## is this necessary?????
    psi = tau - 1*(ehat<0)
    B = diag(Ku(ehat/h))
nt=m
p = ncol(x)
 #Lambda: meat part of the sandwich formula
 # S: bread part of the sandwich formula
    if(p==1)
    {
        psi2 = zstar2 = rep(0, n*nt)  # fill 0 in for those missed measurements
        for(i in 1:n)
        {
            ind = which(subject[,i]==1)
            ind2 = ((i-1)*nt+1): (i*nt)
            psi2[ind2[1:length(ind)]] = psi[ind] 
            zstar2[ind2[1:length(ind)]] = zstar[ind]
        }
        psi.m = matrix(psi2, ncol=n)
        S = 1/(N*h)*sum(zstar^2*B)
        zstar.m = matrix(zstar2, ncol=n)
        Obj = rbind(psi.m, zstar.m)
        Lambda = 1/N*sum(apply(Obj, 2, function(x) crossprod(x[-(1:nt)],x[(1:nt)])^2))
    }
    if(p>1)
    {
        psi2 = rep(0, n*nt)
        zstar2 = matrix(0, nrow=n*nt, ncol=p)
        for(i in 1:n)
        {
            ind = which(subject[,i]==1)
            ind2 = ((i-1)*nt+1): (i*nt)
            psi2[ind2[1:length(ind)]] = psi[ind] 
            zstar2[ind2[1:length(ind)],] = zstar[ind,]
        }
        S = 1/(N*h)*t(zstar2)%*%B%*%zstar2
        Obj = cbind(zstar2, psi2)
        Obj = array(Obj, dim=c(nt, n, p+1)) # the first [,,1:p] are zstar.m, and the last one is psi.m
        Lambda = 1/N*matrix(apply(
                apply(Obj, 2, function(xx) 
                {
                    z1 = xx[,1:p]
                    ppsi = xx[,p+1]
                    out = t(z1)%*%ppsi%*%t(ppsi)%*%z1
                    return(out)
                }), 1, sum), ncol=p)
    }
    
 # the bread part of the sandwitch formula
    inv.S = solve(S)
    COV = 1/N*(inv.S%*%Lambda%*%inv.S)
    se = sqrt(diag(COV))
    t.wald = beta/se
    pval.wald = 2*(1-pnorm(abs(t.wald)))
    lb = beta - qnorm(1-alpha/2)* se
    ub = beta + qnorm(1-alpha/2)* se
    
    stat = cbind(beta, se, t.wald, pval.wald, lb, ub)
    return(list(stat=stat, COV=COV))    
}







############################################################################
########################### Simulation Study ###############################
############################################################################

# Model 1
DGP1 = function(n, m, tau, nsim=2000)
{
#n=30; m=3; tau=0.5; hn=0.1; nsim=5
rho = 0.7
    R = matrix(rho, ncol=m, nrow=m); diag(R)=1  # exchangeable correlation matrix
    subject = rep(1:n, each=m)
    beta = c(1,1)
    set.seed(3333333)
    # Generate the data set
    Data = NULL
    for(i in 1:(nsim))
    {
        x = mvrnorm(n = n,  mu=seq(0.3, m*0.3, length=m), Sigma=0.5*diag(m), tol = 1e-6, empirical = FALSE)
        x = matrix(as.vector(t(x)), ncol=1)
        x = x - mean(x)
        e = mvrnorm(n = n, mu=rep(0,m), Sigma=R, tol = 1e-6, empirical = FALSE)
        e = as.vector(t(e))
        e = e - qnorm(tau)
        x = cbind(1,x)
        y = x%*%beta + e
        tmp = cbind(x[,2], y)
        Data[[i]] = tmp 
    }
    return(Data)
}
DGP2 = function(n, m, tau, nsim=2000)
{
# AR(1) correlation structure
rho = 0.7
    R = diag(m)
    for(j in 1:(m-1))
    {
        for(k in (j+1):m)
        {
            R[j,k] = R[k,j]= rho^(k-j)
        }
    }
    subject = rep(1:n, each=m)
    beta = c(1,1)
    set.seed(12345678)
    # Generate the data set
    Data = NULL
    for(i in 1:(nsim))
    {
        x = mvrnorm(n = n,  mu=seq(0.3, m*0.3, length=m), Sigma=0.5*diag(m), tol = 1e-6, empirical = FALSE)
        x = matrix(as.vector(t(x)), ncol=1)
        x = x - mean(x)
        e = mvrnorm(n = n, mu=rep(0,m), Sigma=R, tol = 1e-6, empirical = FALSE)
        e = as.vector(t(e))
        e = e - qnorm(tau)
        e = 0.25*(1+abs(x))*as.vector(t(e))
        x = cbind(1,x)
        y = x%*%beta + e
        tmp = cbind(x[,2], y)
        Data[[i]] = tmp 
    }
    return(Data)
}

M1.sim.CI.EL = function(n, m, tau, nsim, hn, alpha=0.05)
{
# n=30; m=5; tau=0.5; nsim=5; hn=n^(-0.8); alpha=0.05
    set.seed(1234567)
    Data = DGP1(n, m, tau, nsim)
    SEL = SELc = array(0, dim=c(nsim, 2, 3))
    for(i in 1:nsim)
    {
        dat = Data[[i]]
        y = dat[,2]
        x = cbind(1,dat[,1])
        tmp.EL = CI.ELLDrq(x, y, n, m, tau, hn, alpha)
        SEL[i,,] = tmp.EL$CI.SEL
        SELc[i,,] = tmp.EL$CI.SELc
    }
    # coverage probability and average length
    CP.SEL = apply(SEL, 2, function(x) mean(x[,2]<1 & x[,3]>1))
    CP.SELc = apply(SELc, 2, function(x) mean(x[,2]<1 & x[,3]>1))
    Length.SEL = apply(SEL, 2, function(x) mean(x[,3]-x[,2]))
    Length.SELc = apply(SELc, 2, function(x) mean(x[,3]-x[,2]))
    out = list(SEL=SEL, SELc=SELc, CP.SEL=CP.SEL, CP.SELc=CP.SELc, Length.SEL = Length.SEL, Length.SELc = Length.SELc)
    
    cat("CP of SEL.b0, SEL.b1, SELc.b0, SELc.b1", "\n")
    print(c(CP.SEL, CP.SELc))
    cat("average length of SEL.b0, SEL.b1, SELc.b0, SELc.b1", "\n")
    print(c(Length.SEL, Length.SELc))
    
    return(out)  
}


M2.sim.CI.EL = function(n, m, tau, nsim, hn, alpha=0.05)
{
# n=30; m=5; tau=0.5; nsim=5; hn=n^(-0.8); alpha=0.05
    set.seed(1234567)
    Data = DGP2(n, m, tau, nsim)
    SEL = SELc = array(0, dim=c(nsim, 2, 3))
    for(i in 1:nsim)
    {
        dat = Data[[i]]
        y = dat[,2]
        x = cbind(1,dat[,1])
        tmp.EL = CI.ELLDrq(x, y, n, m, tau, hn, alpha)
        SEL[i,,] = tmp.EL$CI.SEL
        SELc[i,,] = tmp.EL$CI.SELc
    }
    # coverage probability and average length
    CP.SEL = apply(SEL, 2, function(x) mean(x[,2]<1 & x[,3]>1))
    CP.SELc = apply(SELc, 2, function(x) mean(x[,2]<1 & x[,3]>1))
    Length.SEL = apply(SEL, 2, function(x) mean(x[,3]-x[,2]))
    Length.SELc = apply(SELc, 2, function(x) mean(x[,3]-x[,2]))
    out = list(SEL=SEL, SELc=SELc, CP.SEL=CP.SEL, CP.SELc=CP.SELc, Length.SEL = Length.SEL, Length.SELc = Length.SELc)
    
    cat("CP of SEL.b0, SEL.b1, SELc.b0, SELc.b1", "\n")
    print(c(CP.SEL, CP.SELc))
    cat("average length of SEL.b0, SEL.b1, SELc.b0, SELc.b1", "\n")
    print(c(Length.SEL, Length.SELc))
    
    return(out)  
}


Boot = function(x, y, n, m, nboot=500, alpha=0.05)
{
    b.hat = rq(y~x-1, tau=tau)$coef
    index = matrix(1:(n*m), ncol=n)
    beta.b = NULL
    for(j in 1:nboot)
    {
        sub.b = sample(1:n, replace=T)
        ind.b = as.vector(index[,sub.b])
        y.b = y[ind.b]
        x.b = x[ind.b,]
        beta.b = rbind(beta.b, rq(y.b~x.b-1,tau=tau)$coef)
        
    }
    tmp = t(apply(beta.b, 2, quantile, c(0.05, 0.95, 0.05/2, 1-0.05/2, 0.05/4, 1-0.05/4)))
    
    # estimate the bootstrap cov of beta.hat
    Cov = cov(beta.b)
    # obtain the confidence interval for intercept and slope
    se = sqrt(diag(Cov))
    lb = b.hat - qnorm(1-alpha/2)* se
    ub = b.hat + qnorm(1-alpha/2)* se
    # test if (interp, slope)=(1,1)
    tval = (b.hat - c(1,1))^T %*% solve(Cov) %*% (b.hat - c(1,1))
    pval = 1-pchisq(tval, 2)
    tval = rep(tval, 2)
    pval = rep(pval, 2)
    Cov = rep(Cov[1,2],2)
    out = cbind(b.hat, lb, ub, pval, tval, se, Cov, tmp)
    colnames(out) = c("coef", "lb", "ub", "pval(comp test)", "tval(comp test)", "se", "Cov(b0,b1)", "Q05", "Q95", "Q025", "Q975", "Q0125", "Q9875")    
    
    return(out)
}

M1.sim.Boot = function(Data, n, m, tau, nsim, nboot=500, alpha=0.05)
{
# nomnal 95\%
# n=30; m=5; tau=0.5; nsim=4; nboot=20
    set.seed(1234567)
    Data = DGP1(n, m, tau, nsim)
    out = array(0, dim=c(nsim, 2, 13))
    for(i in 1:nsim)
    {
        dat = Data[[i]]
        y = dat[,2]
        x = cbind(1,dat[,1])
        tmp = Boot(x, y, n, m, nboot=nboot, alpha)
        out[i,,] = tmp
        if(i %%200==0) print(i)
    }
    # coverage probability and average length; alpha=0.05
    CP = apply(out, 2, function(x) mean(x[,2]<1 & x[,3]>1))
    CP = c(CP, mean(out[,1,4]>alpha))
    
    Length = apply(out, 2, function(x) mean(x[,3]-x[,2]))
    
    Out = list(Boot=out, CP=CP, Length=Length)
    
    cat("CP (nominal 0.95) of boot for b0, b1, b0&b1 are", "\n")
    print(CP)
    cat("average lengths of boot for b0 and b1 are ", "\n")
    print(Length)
    
    return(Out)  
}

M2.sim.Boot = function(Data, n, m, tau, nsim, nboot=500, alpha=0.05)
{
# nomnal 95\%
# n=30; m=5; tau=0.5; nsim=4; nboot=20
    set.seed(1234567)
    Data = DGP2(n, m, tau, nsim)
    out = array(0, dim=c(nsim, 2, 13))
    for(i in 1:nsim)
    {
        dat = Data[[i]]
        y = dat[,2]
        x = cbind(1,dat[,1])
        tmp = Boot(x, y, n, m, nboot=nboot, alpha)
        out[i,,] = tmp
        if(i %%200==0) print(i)
    }
    # coverage probability and average length; alpha=0.05
    CP = apply(out, 2, function(x) mean(x[,2]<1 & x[,3]>1))
    CP = c(CP, mean(out[,1,4]>alpha))
    
    Length = apply(out, 2, function(x) mean(x[,3]-x[,2]))
    
    Out = list(Boot=out, CP=CP, Length=Length)
    
    cat("CP (nominal 0.95) of boot for b0, b1, b0&b1 are", "\n")
    print(CP)
    cat("average lengths of boot for b0 and b1 are ", "\n")
    print(Length)
    
    return(Out)  
}


sim.Wald = function(Data, n, m, tau, nsim, hn, alpha=0.05)
{
# n=30; m=5; tau=0.5; nsim=5; hn=n^(-0.8); alpha=0.05
    r = rep(1:n, each=m)
    subject = model.matrix(~-1+factor(r)) 
    
    wald2 = matrix(0, ncol=2, nrow=nsim)
    wald1 = array(0, dim=c(2,6, nsim))
    for(i in 1:nsim)
    {
        dat = Data[[i]]
        y = dat[,2]
        x = cbind(1,dat[,1])
        tmp = Wald(y, x, subject, n, m, hn, tau, alpha)
        Cov = tmp$COV
        tmp.wald = tmp$stat
        # test if (interp, slope)=(1,1)
        tval = (tmp.wald[,1] - c(1,1))^T %*% solve(Cov) %*% (tmp.wald[,1] - c(1,1))
        pval = 1-pchisq(tval, 2)
        
        wald1[,,i] = tmp.wald
        wald2[i,] = c(tval, pval)
    }
    # coverage probability for testing (1), (2) and (3)
    CP = c(apply(wald1, 1, function(x) mean(x[5,]<1 & x[6,]>1)), mean(wald2[,2] >= alpha))    
    # average length for testing (1) and (2)
    Length = c(mean(wald1[1,6,]-wald1[1,5,]), mean(wald1[2,6,]-wald1[2,5,]))
    out = list(wald1=wald1, wald2=wald2, CP=CP, Length=Length)
    
    cat("CP of wald for testing (1), (2) and (3) are", "\n")
    print(CP)
    cat("Average length of wald for testing (1) and (2) are ", "\n")
    print(Length)
    
    return(out)  
}