# This code can be used to reproduce the simulation results under the t distribution with 10 degrees of freedom scenario (Figure 6 in the manuscript).
# Adjust line 296 to reproduce the simulation results for the different combinations of \phi_1 and \theta.
# In line 296, the ci.sim(nsim, n, p, s1, s2, r1, r2, g=0, intercept, ar.coef, B, conf.level) function uses the following arguments:
# nsim: number of simulation replicates (nsim=500)
# n: the tree size (n = 31, or n = 63, or n = 127)
# p: the order of the bifurcating autoegressive model; first order is used (p=1)
# s1, s2: the standard deviations elements in the variance-covariance errors matrix (s1= 1, s2= 1)
# r1, r2: the error correlation elements in the variance-covariance errors matrix which represent the values of \theta (r1 = -0.9, r2 = -0.9)
# g: degrees of freedom (g = 10)
# intercept: the value of \phi_0 (intercept = 10 )
# ar.coef: the value of \phi_1 (ar.coef = c(-0.85))
# B: number of simulations of single bootstrap  and fast double bootstrap approachs (B = 499)
# conf.level: confidence interval level (conf.level = 0.95)

library(parallel)
library(dplyr)
library(tidyr)
library(bifurcatingr)
library(fCopulae)

# Generating trees based on t distribution errors with 10 degrees of freedom.
bfa.tree.gen.t <- function(n, p, s1, s2, r1, r2, g, intercept, ar.coef){
  if (sum(ar.coef)==1) {stop("Error: the sum of the coefficients must not equal 1")}
  #Use start-up BAR(1) tree
  y <- rep(NA,n)
  y[1] <- intercept/(1-sum(ar.coef))
  rmvst(n, dim = 2, mu = rep(0, 2), Omega = matrix(c(s1,r1,r2,s2), 2, 2), alpha = rep(0, 2), df = g)
  error <- rep(NA,n)
  for(i in 1:((n-1)/2)){
    error[(i*2)] <- e[i,1]
    error[(i*2)+1] <- e[i,2]
  }
  for(j in 2:n){
    y[j] <- intercept+ar.coef[1]*y[floor(j/(2^1))]+error[j]
  }
  #Generate BAR(p) tree using last subtree with (2^p-1) y-values as starting values
  term.matrix <- matrix(rep(NA,p*n),ncol=p)
  tree <- rep(NA,n)
  tree[1:((2^p)-1)] <- y[bfa.subtree(n,p)]
  e<-rmvsnorm((n-1)/2, dim = 2, mu = rep(0, 2), Omega = matrix(c(s1,r1,r2,s2), 2, 2), alpha = rep(g, 2))
  error <- rep(NA,n)
  for(i in 1:((n-1)/2)){
    error[(i*2)] <- e[i,1]
    error[(i*2)+1] <- e[i,2]
  }
  for(i in (2^p):n){
    for(k in 1:p){
      term.matrix[,k] <- ar.coef[k]*tree[floor(i/(2^k))]
    }
    tree[i] <- intercept+sum(term.matrix[i,])+error[i]
  }
  return(tree)
}

# Single bootstrap procedure for the LS estimates of BAR(1) models
bfa.boot.ls.bc <- function(z, p, method="boot1",burn = 5, B, boot.est=TRUE, boot.data=FALSE){  # burn can be determined by simulation
  n=length(z)
  boot.dat <- matrix(rep(NA, n*B), ncol=n)
  est <- bfa.ls(z,p,resids=TRUE)
  intercept <- est$coef[1,1]
  coef <- t(as.vector(est$coef[1,2:(p+1)]))
  e <- est$resids-mean(est$resids,na.rm=TRUE)
  e.e=matrix(e,ncol=2, byrow=T)
  ee <- e.e[stats::complete.cases(e.e),]
  for (h in 1:B){
    # start-up tree
    n0 = 2^(p + burn) - 1
    term.matrix0 <- matrix(rep(NA,p*n0),ncol=p)
    dat0 <- rep(NA,n0)
    er0 <- matrix(rep(NA,((n0-1)/2)),ncol=2,nrow=((n0-1)/2))
    samp0 <- sample(1:nrow(ee), ((n0-1)/2) ,replace=T)
    for (s in 1: nrow(er0)) {
      l <- samp0[s]
      er0[s,] <- ee[l,]
    }
    error0 <- rep(NA,n0)
    for(i in 1:((n0-1)/2)){
      error0[(i*2)] <- er0[i,1]
      error0[(i*2)+1] <- er0[i,2]
    }
    dat0[1:((2^p)-1)] <- z[1:((2^p)-1)]
    for(j in (2^p) : n0){
      for(k in 1:p){
        term.matrix0[,k] <- coef[1,k]*dat0[floor(j/(2^k))]
      }
      dat0[j] <- intercept+sum(term.matrix0[j,])+error0[j]
    }
    # bootstrap tree
    term.matrix <- matrix(rep(NA,p*n),ncol=p)
    dat <- rep(NA,n)
    er <- matrix(rep(NA,((n-1)/2)),ncol=2,nrow=((n-1)/2))
    samp <- sample(1:nrow(ee), ((n-1)/2) ,replace=T)
    for (s in 1: nrow(er)) {
      l <- samp[s]
      er[s,] <- ee[l,]
    }
    error <- rep(NA,n)
    for(i in 1:((n-1)/2)){
      error[(i*2)] <- er[i,1]
      error[(i*2)+1] <- er[i,2]
    }
    dat[1:((2^p)-1)] <- dat0[bfa.subtree(n0,p)]
    for(j in (2^p) : n){
      for(k in 1:p){
        term.matrix[,k] <- coef[1,k]*dat[floor(j/(2^k))]
      }
      dat[j] <- intercept+sum(term.matrix[j,])+error[j]
    }
    # copying missing values as in the original tree
    for (k in 1:n){
      dat[k] <- ifelse(is.na(z[k])==T, NA, dat[k])
    }
    boot.dat[h,] <- dat
  }
  out <- list()
  if(boot.est)
    out$boot.bcest=matrix(unlist(apply(boot.dat, 1, bfa.ls.bc, p=p, method=method)),nrow=B, byrow=TRUE)
  if(boot.data)
    out$boot.data <- boot.dat
  return(out)
}

##Asymptotic CI
bfa.asy.ci <- function(a1.ls, a1.ls.star, conf.level){
  alpha = 1 - conf.level
  se = sd(a1.ls.star)
  a1.L = a1.ls - qnorm(1-alpha/2)*se
  a1.U = a1.ls + qnorm(1-alpha/2)*se
  return(asy.ci = c(a1.L,a1.U))
}

##Percentile CI
bfa.perc.ci <- function(a1.ls.star, conf.level){
  alpha = 1 - conf.level
  return(perc.ci = quantile(a1.ls.star, prob=c(alpha/2, 1-alpha/2)))
}

##BC Percentile CI
bfa.perc.bc.ci <- function(z, a1.ls, a1.ls.star, conf.level) {
  # generate the bootstrap samples
  n <- length(z)
  B <- length(a1.ls.star)
  # compute z.hat
  z.hat <- qnorm(sum(a1.ls.star < a1.ls)/B)
  # pick out the BC quantiles
  alpha = 1 - conf.level
  p <- c(alpha/2, 1 - alpha/2)
  w <- 2*z.hat + qnorm(p)
  q <- pnorm(w)
  ind <- ceiling(q*B)
  a1.ls.star <- sort(a1.ls.star)
  return(percbc.ci = a1.ls.star[ind])
}

#########################
#BCa Function
# phi1 is the estimated phi
#phis is the vector of bootstrap replicates
##############################
Boot.BCa <- function(z,p,B=99,J=2,conf.level=0.95){
  n=length(z)
  nn=floor((n-1)/J)
  Rs=c(J:((2*J)-1))
  M=length(Rs)
  phi1=bfa.ls(z,1)$coef[1,2]
  dat.matrx = matrix(NA, nrow=J, ncol=nn)
  dat=list()
  # Bootstarp LS estimates  
  phis <- matrix(unlist(apply(matrix(unlist(bfa.boot1.ls(z,1,B=B, boot.est=FALSE,boot.data = TRUE)), ncol=n, nrow = B),1,p=p, bfa.ls)),ncol=3,byrow=TRUE)[,2]
  alpha <- (1 + c(-conf.level, conf.level))/2 
  zalpha <- qnorm(alpha)
  # the bias correction factor
  z0 <- qnorm(sum(as.vector(phis) < phi1) / B)
  # the acceleration factor (jackknife est.) 
  for(R in Rs){
    indx=NULL
    for (i in 0:10){
      ind=((2^i)*R):((2^i)*R+(2^i-1))
      indx=c(indx,if(tail(ind, n=1)<=n){ind})
    }
    dat=append(dat, list(indx))
  }
  for (i in 1:M){
    dat.matrx[i,] = z[dat[[i]]]
  }
  phi.jack <- matrix(unlist(apply(dat.matrx,1,p=p, bfa.ls)),ncol=3,byrow=TRUE)[,2]
  L <- mean(phi.jack) - phi.jack
  a <- sum(L^3)/(6 * sum(L^2)^1.5)
  # BCa conf. limits
  adj.alpha <- pnorm(z0 + (z0+zalpha)/(1-a*(z0+zalpha))) 
  limits <- quantile(phis, adj.alpha, names=F, type=8,na.rm = TRUE)
  return(list("est"=phi1, "BCa"=limits))
  
}

ci.sim <- function(nsim, n, p, s1, s2, r1, r2, g=0, intercept, ar.coef, B, conf.level){
  final_ans = matrix(NA,ncol=55,nrow=nsim)
  i <- seq(1,nsim)
  ci.itr <- function(itr){
    z = bfa.tree.gen.t(n, p, s1, s2, r1, r2, g, intercept, ar.coef)
    n = length(z)
    ls.fit = bfa.ls(z, p, x.data = FALSE, y.data = FALSE, resids = TRUE, error.cor = TRUE, error.var = FALSE,
                    cov.matrix = FALSE, conf = TRUE, conf.level = conf.level, p.value = FALSE)
    a0.ls = ls.fit$coef[1,1]
    a1.ls = ls.fit$coef[1,2]
    res.ls = ls.fit$resids
    a1.ls.boot1 =bfa.ls.bc(z,p,method="boot1", B1=99, B2=49)
    a1.ls.boot2f = bfa.ls.bc(z,p,method="boot2fast", B1=99, B2=49)
    a1.ls.star = bfa.boot1.ls(z, p, burn = 5, B, boot.est = TRUE, boot.data = FALSE)
    a1.ls.boot1.star = bfa.boot.ls.bc(z, p, method="boot1",burn = 5, B, boot.est=TRUE, boot.data=FALSE)
    a1.ls.boot2f.star = bfa.boot.ls.bc(z, p, method="boot2fast",burn = 5, B, boot.est=TRUE, boot.data=FALSE)

    asy.fm.ci=ls.fit$conf
    asy.fm.ind.l=ifelse(ar.coef<asy.fm.ci[2,1], 1, 0)
    asy.fm.ind.u=ifelse(ar.coef>asy.fm.ci[2,2], 1, 0)
    asy.fm.ind.t = ifelse(ar.coef>=asy.fm.ci[2,1] & ar.coef<=asy.fm.ci[2,2], 1, 0)
    asy.fm.w=asy.fm.ci[2,2]-asy.fm.ci[2,1]

    asy.ci = bfa.asy.ci(a1.ls, a1.ls.star$boot.est[,2], conf.level)
    asy.ind.l = ifelse(ar.coef<asy.ci[1], 1, 0)
    asy.ind.u = ifelse(ar.coef>asy.ci[2], 1, 0)
    asy.ind.t = ifelse(ar.coef>=asy.ci[1] & ar.coef<=asy.ci[2], 1, 0)
    asy.w = asy.ci[2] - asy.ci[1]

    asy.boot1.ci = bfa.asy.ci(a1.ls.boot1, a1.ls.boot1.star$boot.bcest[,1], conf.level)
    asy.boot1.ind.l = ifelse(ar.coef<asy.boot1.ci[1], 1, 0)
    asy.boot1.ind.u = ifelse(ar.coef>asy.boot1.ci[2], 1, 0)
    asy.boot1.ind.t = ifelse(ar.coef>=asy.boot1.ci[1] & ar.coef<=asy.boot1.ci[2], 1, 0)
    asy.boot1.w = asy.boot1.ci[2] - asy.boot1.ci[1]

    asy.boot2f.ci = bfa.asy.ci(a1.ls.boot2f, a1.ls.boot2f.star$boot.bcest[,1], conf.level)
    asy.boot2f.ind.l = ifelse(ar.coef<asy.boot2f.ci[1], 1, 0)
    asy.boot2f.ind.u = ifelse(ar.coef>asy.boot2f.ci[2], 1, 0)
    asy.boot2f.ind.t = ifelse(ar.coef>=asy.boot2f.ci[1] & ar.coef<=asy.boot2f.ci[2], 1, 0)
    asy.boot2f.w = asy.boot2f.ci[2] - asy.boot2f.ci[1]

    perc.ci = bfa.perc.ci(a1.ls.star$boot.est[,2], conf.level)
    perc.ind.l = ifelse(ar.coef<perc.ci[1], 1, 0)
    perc.ind.u = ifelse(ar.coef>perc.ci[2], 1, 0)
    perc.ind.t = ifelse(ar.coef>=perc.ci[1] & ar.coef<=perc.ci[2], 1, 0)
    perc.w = perc.ci[2] - perc.ci[1]

    perc.boot1.ci = bfa.perc.ci(a1.ls.boot1.star$boot.bcest[,1], conf.level)
    perc.boot1.ind.l = ifelse(ar.coef<perc.boot1.ci[1], 1, 0)
    perc.boot1.ind.u = ifelse(ar.coef>perc.boot1.ci[2], 1, 0)
    perc.boot1.ind.t = ifelse(ar.coef>=perc.boot1.ci[1] & ar.coef<=perc.boot1.ci[2], 1, 0)
    perc.boot1.w = perc.boot1.ci[2] - perc.boot1.ci[1]

    perc.boot2f.ci = bfa.perc.ci(a1.ls.boot2f.star$boot.bcest[,1], conf.level)
    perc.boot2f.ind.l = ifelse(ar.coef<perc.boot2f.ci[1], 1, 0)
    perc.boot2f.ind.u = ifelse(ar.coef>perc.boot2f.ci[2], 1, 0)
    perc.boot2f.ind.t = ifelse(ar.coef>=perc.boot2f.ci[1] & ar.coef<=perc.boot2f.ci[2], 1, 0)
    perc.boot2f.w = perc.boot2f.ci[2] - perc.boot2f.ci[1]

    percbc.ci = bfa.perc.bc.ci(z, a1.ls, a1.ls.star$boot.est[,2], conf.level)
    percbc.ind.l = ifelse(ar.coef<percbc.ci[1], 1, 0)
    percbc.ind.u = ifelse(ar.coef>percbc.ci[2], 1, 0)
    percbc.ind.t = ifelse(ar.coef>=percbc.ci[1] & ar.coef<=percbc.ci[2], 1, 0)
    percbc.w = percbc.ci[2] - percbc.ci[1]

    bca.ci=Boot.BCa(z,p,B,J=4,conf.level)
    bca.ind.l=ifelse(ar.coef<bca.ci$BCa[1], 1, 0)
    bca.ind.u=ifelse(ar.coef>bca.ci$BCa[2], 1, 0)
    bca.ind.t = ifelse(ar.coef>=bca.ci$BCa[1] & ar.coef<=bca.ci$BCa[2], 1, 0)
    bca.w=bca.ci$BCa[2]-bca.ci$BCa[1]
    
    return(list(ar.coef, r1, n,
                asy.fm.ind.l, asy.fm.ind.u, asy.fm.ind.t, asy.fm.w,
                asy.ind.l, asy.ind.u, asy.ind.t, asy.w,
                asy.boot1.ind.l, asy.boot1.ind.u, asy.boot1.ind.t, asy.boot1.w,
                asy.boot2f.ind.l, asy.boot2f.ind.u, asy.boot2f.ind.t, asy.boot2f.w,
                perc.ind.l, perc.ind.u, perc.ind.t, perc.w,
                perc.boot1.ind.l, perc.boot1.ind.u, perc.boot1.ind.t, perc.boot1.w,
                perc.boot2f.ind.l, perc.boot2f.ind.u, perc.boot2f.ind.t, perc.boot2f.w,
                percbc.ind.l, percbc.ind.u, percbc.ind.t, percbc.w,
                bca.ind.l, bca.ind.u, bca.ind.t, bca.w))
  }
  numCores <- detectCores()
  print(paste("Cores=",numCores))
  results <- parallel::mclapply(i, ci.itr, mc.cores = numCores)
  final_ans <- do.call(rbind, results)
  colnames(final_ans) <- c("ar.coef","r1","n",
                           "asy.fm.ind.l", "asy.fm.ind.u", "asy.fm.ind.t", "asy.fm.w",
                           "asy.ind.l", "asy.ind.u", "asy.ind.t", "asy.w",
                           "asy.boot1.ind.l", "asy.boot1.ind.u", "asy.boot1.ind.t", "asy.boot1.w",
                           "asy.boot2f.ind.l", "asy.boot2f.ind.u", "asy.boot2f.ind.t", "asy.boot2f.w",
                           "perc.ind.l", "perc.ind.u", "perc.ind.t", "perc.w",
                           "perc.boot1.ind.l", "perc.boot1.ind.u", "perc.boot1.ind.t", "perc.boot1.w",
                           "perc.boot2f.ind.l", "perc.boot2f.ind.u", "perc.boot2f.ind.t", "perc.boot2f.w",
                           "percbc.ind.l", "percbc.ind.u", "percbc.ind.t", "percbc.w",
                           "bca.ind.l", "bca.ind.u", "bca.ind.t", "bca.w")
  final_ans <- as.data.frame(final_ans)
  return(final_ans)
}

dat <- ci.sim(nsim=500,n=63,p=1,s1=1,s2=1,r1=0.9,r2=0.9,g=10,intercept=10,ar.coef=c(-0.85),B=499, conf.level=0.95) # Change (r1, r2, and ar.coef) to simulate the other combinations
dat <- unnest(dat,cols = c(ar.coef, r1, n,
                               asy.fm.ind.l, asy.fm.ind.u, asy.fm.ind.t, asy.fm.w,
                               asy.ind.l, asy.ind.u, asy.ind.t, asy.w,
                               asy.boot1.ind.l, asy.boot1.ind.u, asy.boot1.ind.t, asy.boot1.w,
                               asy.boot2f.ind.l, asy.boot2f.ind.u, asy.boot2f.ind.t, asy.boot2f.w,
                               perc.ind.l, perc.ind.u, perc.ind.t, perc.w,
                               perc.boot1.ind.l, perc.boot1.ind.u, perc.boot1.ind.t, perc.boot1.w,
                               perc.boot2f.ind.l, perc.boot2f.ind.u, perc.boot2f.ind.t, perc.boot2f.w,
                               percbc.ind.l, percbc.ind.u, percbc.ind.t, percbc.w,
                               bca.ind.l, bca.ind.u, bca.ind.t, bca.w))
dat.summary <- dat %>%
  summarize(ar.coef=mean(ar.coef),r1=mean(r1),n=mean(n),
            asy.fm.ind.l=mean(asy.fm.ind.l), asy.fm.ind.u=mean(asy.fm.ind.u), asy.fm.ind.t=mean(asy.fm.ind.t), asy.fm.w=mean(asy.fm.w),
            asy.ind.l=mean(asy.ind.l), asy.ind.u=mean(asy.ind.u), asy.ind.t=mean(asy.ind.t), asy.w=mean(asy.w),
            asy.boot1.ind.l=mean(asy.boot1.ind.l), asy.boot1.ind.u=mean(asy.boot1.ind.u), asy.boot1.ind.t=mean(asy.boot1.ind.t), asy.boot1.w=mean(asy.boot1.w),
            asy.boot2f.ind.l=mean(asy.boot2f.ind.l), asy.boot2f.ind.u=mean(asy.boot2f.ind.u), asy.boot2f.ind.t=mean(asy.boot2f.ind.t), asy.boot2f.w=mean(asy.boot2f.w),
            perc.ind.l=mean(perc.ind.l), perc.ind.u=mean(perc.ind.u), perc.ind.t=mean(perc.ind.t), perc.w=mean(perc.w),
            perc.boot1.ind.l=mean(perc.boot1.ind.l), perc.boot1.ind.u=mean(perc.boot1.ind.u), perc.boot1.ind.t=mean(perc.boot1.ind.t), perc.boot1.w=mean(perc.boot1.w),
            perc.boot2f.ind.l=mean(perc.boot2f.ind.l), perc.boot2f.ind.u=mean(perc.boot2f.ind.u), perc.boot2f.ind.t=mean(perc.boot2f.ind.t), perc.boot2f.w=mean(perc.boot2f.w),
            percbc.ind.l=mean(percbc.ind.l), percbc.ind.u=mean(percbc.ind.u), percbc.ind.t=mean(percbc.ind.t), percbc.w=mean(percbc.w),
            bca.ind.l=mean(bca.ind.l), bca.ind.u=mean(bca.ind.u), bca.ind.t=mean(bca.ind.t), bca.w=mean(bca.w))

outputName=paste("senario-",5,".csv",sep="")
outputPath=file.path("/nas/longleaf/home/tmsy/ProjectFoo/Output2",outputName)
write.csv(dat.summary,outputPath, row.names = FALSE)


