library(quantreg)
library(splines)
library(cobs)
library(sn)


source('bivariate quantile functions.R')

# simulate data from a nonlinear model with location-scale errors.

set.seed(666666)
n = 5000
t = runif(n,0,3)#
h = 40*t/(1+2*t)+ (1+(t/5))*rst(n,shape=5)
w = log(1+t) +(1+0.2*t)*h + (1+(h/10))*rnorm(n)


# -----Step 1:  Estimation of stratified quantile models -----------------------------------#


    knots = c(0.5,1,1.5,2,2.5); ord = 4;  # internal knots and order of B-splines   
    boundry.knots=range(t);knots=c(rep(boundry.knots[1],ord),knots,rep(boundry.knots[2],ord))   
    xm = splineDesign(knots=knots,t,ord=ord)  # design matrix for estimating one coefficient function
    ntau = 200                                # number of quantile levels  
    tau = seq(0,1,length = ntau + 1)[-c(1, ntau + 1)]; 
    
    # modelling order 1

      # (a.1) one marginal model of h ---  Q_tau(h) = g_1(t)   --- #
        fit1 =rq(h ~ xm - 1,tau=tau)
      # (b.1) one conditional coefficient-varying model of w given h --- w = alpha(t) + beta(t)*h --- #
         xm2 =  cbind(xm, xm*h)
         fit2 = rq(w~ xm2 - 1,tau=tau)

    # modelling order 2
    
      # (a.2) one marginal model of w --- #    
        fit3 =rq(w ~ xm - 1,tau=tau)
      # (b.2) one conditional coefficient-varying model of h given w.    --- #
        xm2 =  cbind(xm, xm*w)
        fit4 = rq(h~ xm2 - 1,tau=tau)


# ------Step 2:  model based simulation at a target time t0 -----------------------------------#

        t0 = 0.5  # t0 is the target time
        use = 1 
           # use = 1, use model order 1
           # use = 2, use model order 2
           # use = 3, use model order 3, combine the 2 orders
        newxm = as.vector(splineDesign(knots = knots,t0,ord = ord))
        index = 1:length(tau); 
        
        set.seed(123456)
        #simulate data based on the models of order 1 
        simdata1 = NULL
        for(m in 1:10000){
            simh = newxm%*%fit1$coef[,sample(index,1)]
            newxm2 = c(newxm,newxm*simh)
            simw = newxm2%*%fit2$coef[,sample(index,1)]
            simdata1 = rbind(simdata1,c(simh,simw))        }

       #simulate data based on the models of order 2 
        simdata2 = NULL
        for(m in 1:10000){
            simw = newxm%*%fit3$coef[,sample(index,1)]
            newxm2 = c(newxm,newxm*simw)
            simh = newxm2%*%fit4$coef[,sample(index,1)]
            simdata2 = rbind(simdata2,c(simh,simw))        }
      
       
       if(use == 1){x1=simdata1[,1]; x2=simdata1[,2]}
       if(use == 2){x1=simdata2[,1]; x2=simdata2[,2]}
       if(use == 3){x1=c(simdata1[,1],simdata2[,1]); 
                         x2=c(simdata1[,2],simdata2[,2])}

# -------Step 3: Estimation of conditional bivariate quantile contours-----------------------------#

       coverage =  c(0.5,0.75,0.95) 
       
       CM = NULL
       for(i in 1:length(coverage)){
             con = RefCountour(x1,x2,coverage=coverage[i], nknots=10, standardize = T) 
             Cont = as.data.frame(cbind(con$X, con$Y))
             colnames(Cont) = c("x","y")
             CM = rbind(CM, cbind(t0,coverage[i],Cont))
       }  
       
            
      

# -------Step 4: Plot the estimated contours -----------------------------------#
        
    
    
    win.graph(height=4,width=4)
    par(mfcol=c(1,1),mgp=c(2,0.6,0),mar=c(3,3,1,1))   
    
    xlim = range(CM[,3]);
    ylim = range(CM[,4]);
    plot(1,1,type = "n",  xlab="x1",ylab="x2",xlim =xlim, ylim = ylim)
    for(i in 1:length(coverage)){
          Cont = as.data.frame(CM[CM[,2]==coverage[i],3:4])
          lines(Cont[,1],Cont[,2],lty=1)}
     points(con$center[1],con$center[2])