R code for performing quantile regression using bivariate copulas

Nelsen explained in his 1999 book An introduction to copulas how to fit a (nonlinear) quantile regression model by means of a bivariate copula (pp. 175-176).

nonlinear-quantile-regression-copulas

In short, Nelsen’s method for fitting a (nonlinear) quantile regression model is as follows:

  1. Take the partial derivative of the copula function C(u, v) with respect to u, where u and v are both defined in [0, 1]. Denote this partial derivative by cu(u, v), and note that cu(u, v) = P{V ≤ v | U = u}.
  2. For fitting the quantile regression model, set cu(u, v) = p, where p is defined in [0, 1].
  3. For regressing v on u, evaluate cu(u, v) at u, and subsequently solve cu(u, v) = p for v (or, similarly, solve 0 = cu(u, v) – p for v).
  4. Replace u by Fx-1(u) and v by Fy-1(v), where Fx-1(·) and Fy-1(·) are the quantile functions for x (=independent variable) and y (=dependent variable), respectively.

The following R code implements this copula method proposed by Nelsen for fitting a (nonlinear) quantile regression model. In addition, the R code may also compute confidence intervals for the fitted quantiles using a Monte-Carlo method.

In this blog post I will demonstrate how to fit a multiple (nonlinear) quantile regression model by means of a copula.

Suggestions and/or questions? Please contact Stefan Gelissen (email: info at datall-analyse.nl).
See this page for an overview of all of Stefan’s
R code blog posts.

R code

library(copula)
library(MASS)
library(mvtnorm)

#Function for performing (nonlinear) quantile regression
#by means of a fitted bivariate copula
qrCopula <- function(copula, p=.5, marginal, x) {
  #-copula is a fitted bivariate copula as returned by fitMvdc().
  #-p is a number specifying the quantile to be estimated, where p
  # should be within the range [0, 1].
  #-marginal: either 1 or 2. If the second marginal is the dependent
  # variable in the quantile regression model, then use marginal=2.
  # On the other hand, if the first marginal is treated as the dependent
  # variable, then marginal=1. See Example 1 below for details.
  #-x is the data point at which to estimate the quantile. See Example 1
  # below for details.
  
  mvdcQR <- copula
  
  #copula function
  if (class(mvdcQR@mvdc@copula)[1]=="normalCopula") {NULL}
  else if (class(mvdcQR@mvdc@copula)[1]=="tCopula") {NULL}
  else if (class(mvdcQR@mvdc@copula)[1]=="joeCopula") {
    biCopFun <- expression(cdf = 1 - ( (1-u1)^alpha + (1-u2)^alpha -
                                         ((1-u1)^alpha) * ((1-u2)^alpha) )^1/alpha)}
  else {biCopFun <- mvdcQR@mvdc@copula@exprdist[1]}
  
  #partial derivative of copula function with respect to u1
  if (class(mvdcQR@mvdc@copula)[1]=="normalCopula") {NULL}
  else if (class(mvdcQR@mvdc@copula)[1]=="tCopula") {NULL}
  else {dC.du1 <- deriv(biCopFun, "u1")}
  
  i <- marginal #dependent variable
  j <- ifelse(i==1, 2, 1) #independent variable
  
  #fitted copula parameters
  alpha <- mvdcQR@mvdc@copula@parameters
  
  #construct cdf and quantile function for marginal distributions
  cdf.margin <- copula:::asCall(paste0("p", mvdcQR@mvdc@margins[j]),
                                mvdcQR@mvdc@paramMargins[[j]])
  qdf.margin <- copula:::asCall(paste0("q", mvdcQR@mvdc@margins[i]),
                                mvdcQR@mvdc@paramMargins[[i]])
  u1 <- eval(cdf.margin, list(x=x))
  
  #compute u2
  if (class(mvdcQR@mvdc@copula)[1]=="normalCopula"){
    u2 <- pnorm(alpha*qnorm(u1) + sqrt(1-alpha^2)*qnorm(p))}
  else if (class(mvdcQR@mvdc@copula)[1]=="tCopula"){
    u2 <- cCopula(c(u1, p), copula=mvdcQR@mvdc@copula, inverse=TRUE)[2]}
  else {
    u2 <- NA
    u2fun <- function (u2, u1, alpha, t) {t - attr(eval(dC.du1),"gradient")}
    try(u2 <- uniroot(u2fun, interval=c(1-.99999999999999, .99999999999999),
                      u1=u1, alpha=alpha, t=p, tol=1e-5)$root, silent=FALSE)}
  x2 <- eval(qdf.margin, list(x=u2))
  x2
}

#Function for estimating a confidence interval for the p-th quantile
#using a Monte-Carlo method
qrCopulaMC <- function(copula, p=.5, marginal, x, n) {
  mvdcMC <- copula
  #-copula is a fitted bivariate copula as returned by fitMvdc().
  #-p is a number specifying the quantile to be estimated, where p
  # should be within the range [0, 1].
  #-marginal: either 1 or 2. If the second marginal is the dependent
  # variable in the quantile regression model, then use marginal=2.
  # On the other hand, if the first marginal is treated as the dependent
  # variable, then marginal=1. See Example 1 below for details.
  #-x is the data point at which to estimate the quantile. See Example 1
  # below for details.
  #-n is the number of Monte-Carlo samples
  
  #Draw n samples from the copula estimates.
  #Here I assume that the copula estimates follow a multivariate
  #Gaussian distribution.
  sampleValues <- rmvnorm(n, mean=coef(mvdcMC), sigma=vcov(mvdcMC))
  #However, for smaller numbers of observations, it may be reasonable to
  #doubt the assumption that the copula estimates follow a multivariate
  #distribution, and it might be better to resort to a bootstrap method.
  
  #copula function
  if (class(mvdcMC@mvdc@copula)[1]=="normalCopula") {NULL}
  else if (class(mvdcMC@mvdc@copula)[1]=="tCopula") {
    stop("t copula is not implemented")
  }
  else if (class(mvdcMC@mvdc@copula)[1]=="joeCopula") {
    biCopFun <- expression(cdf = 1 - ( (1-u1)^alpha + (1-u2)^alpha -
                                         ((1-u1)^alpha) * ((1-u2)^alpha) )^1/alpha)}
  else {biCopFun <- mvdcMC@mvdc@copula@exprdist[1]}
  
  #partial derivative of copula function with respect to u1
  if (class(mvdcMC@mvdc@copula)[1]=="normalCopula") {NULL}
  else if (class(mvdcMC@mvdc@copula)[1]=="tCopula") {NULL}
  else {dC.du1 <- deriv(biCopFun, "u1")}
  
  i <- marginal #dependent variable
  j <- ifelse(i==1, 2, 1) #independent variable
  
  sapply(1:n, function(k) {
    #replace fitted marginal estimates with sampled ones
    nmp <- length(unlist(mvdcMC@mvdc@paramMargins))
    temp1 <- as.relistable(mvdcMC@mvdc@paramMargins)
    temp2 <- unlist(temp1)
    temp2[1:nmp] <- sampleValues[k, 1:nmp]
    marginsParam <- relist(temp2)
    attr(marginsParam, "class") <- NULL
    mvdcMC@mvdc@paramMargins <- marginsParam
    
    #replace fitted copula parameters with sampled ones
    alpha <- sampleValues[k, -c(1:nmp)]
    
    #construct cdf and quantile function for marginal distributions
    cdf.margin <- copula:::asCall(paste0("p", mvdcMC@mvdc@margins[j]),
                                  mvdcMC@mvdc@paramMargins[[j]])
    qdf.margin <- copula:::asCall(paste0("q", mvdcMC@mvdc@margins[i]),
                                  mvdcMC@mvdc@paramMargins[[i]])
    u1 <- eval(cdf.margin, list(x=x))
    
    #compute u2
    if (class(mvdcMC@mvdc@copula)[1]=="normalCopula"){
      u2 <- pnorm(alpha*qnorm(u1) + sqrt(1-alpha^2)*qnorm(p))}
    else if (class(mvdcMC@mvdc@copula)[1]=="tCopula"){NULL}
    else {
      u2 <- NA
      u2fun <- function (u2, u1, alpha, t) {t - attr(eval(dC.du1),"gradient")}
      try(u2 <- uniroot(u2fun, interval=c(1-.99999999999999, .99999999999999),
                        u1=u1, alpha=alpha, t=p, tol=1e-5)$root, silent=FALSE)}
    x2 <- eval(qdf.margin, list(x=u2))
    x2})
}



##EXAMPLE 1: performing nonlinear quantile regression

##Generate some artificial data using a Gumbel copula

#Bivariate Gumbel copula
myCop <- gumbelCopula(param=3, dim=2)

#Marginals:
#the first marginal is a normal distribution, and
#the second marginal is a gamma distribution
myMvdc <- mvdc(copula=myCop, c("norm", "gamma"),
               list(list(mean=5, sd=.02), list(shape=1.2, scale=3)))

#Visualize the copula
persp(myMvdc, dMvdc, xlim=c(4.95, 5.1), ylim=c(0, 30), phi=40, theta=140)
contour(myMvdc, dMvdc, xlim=c(4.95, 5.1), ylim=c(0, 30))

#Draw sample (this sample will subsequently be used as the artificial data)
ns <- 200 #sample size
xys <- rMvdc(ns, myMvdc)

#Inspect data
plot(xys[,1], xys[,2], xlab="x", ylab="y")



##Fit a bivariate Gumbel copula to the artificial data

#Obtain start values for fitting the copula
normMLE <- fitdistr(xys[,2], densfun="normal")
normMean1 <- normMLE$estimate[1]
normSD1 <- normMLE$estimate[2]

gammaMLE <- fitdistr(xys[,2], densfun="gamma")
gammaShape2 <- gammaMLE$estimate[1]
gammaScale2 <- 1/gammaMLE$estimate[2]

#Starting value for the copula parameter
tauEst <- cor(xys[,1], xys[,2], method="kendall")
thetaGumbel <- 1 / (1-tauEst)

#Fit a Gumbel copula to the artificial data
myGumbel <- mvdc(copula=gumbelCopula(param=1.5, dim=2), c("norm", "gamma"),
                 list(list(mean=0, sd=1), list(shape=1, scale=1)))

start <- c(normMean1, normSD1, gammaShape2, gammaScale2, thetaGumbel)

fitGumbel <- fitMvdc(xys, myGumbel, start=start)
fitGumbel



##Quantile regression: regress y on x
xpoints <- seq(4.90, 5.10, length.out=100)
#compute the median
ymed <- sapply(xpoints, function (k) qrCopula(copula=fitGumbel, p=.5,
                                              marginal=2, x=k))
#compute the .05th-quantile
y05 <- sapply(xpoints, function (k) qrCopula(copula=fitGumbel, p=.05,
                                             marginal=2, x=k))
#compute the .95th-quantile
y95 <- sapply(xpoints, function (k) qrCopula(copula=fitGumbel, p=.95,
                                             marginal=2, x=k))


#Plot the results
plot(xys[,1], xys[,2], xlab="x", ylab="y", col=gray(.2,.2))
#add the median to the plot
lines(xpoints, ymed, lty=1, col="red")
#add the 90% prediction interval: approximately 20 (out of 200) points
#should lie outside the boundaries of the 90% prediction interval
lines(xpoints, y05, lty=2, col="blue")
lines(xpoints, y95, lty=2, col="blue")
#add a legend
legend("topleft", lty=c(1, 2),
       col=c("red", "blue"),
       legend=c("median", "90% prediction interval"),
       bty="n", y.intersp=1.2, cex=.7)


#Use a Monte-Carlo method for computing a confidence interval
#for the median at some given x-value
x <- 5.01
predsY <- qrCopulaMC(copula=fitGumbel, p=.5, marginal=2, x=x, n=5000)
(qsY <- quantile(x=predsY, probs=c(.05, .95), na.rm=TRUE))
#add the confidence interval to the plot
points(x, qrCopula(copula=fitGumbel, p=.5, marginal=2, x=x), col="darkgreen")
arrows(x, qsY[1], x , qsY[2], code=3, length=0.1, angle=90, col="darkgreen")

#Compute the confidence interval for the median at some other points
xpoints2 <- seq(4.95, 5.04, length.out=5)
multy <- sapply(xpoints2, function (k) {
  preds <- qrCopulaMC(copula=fitGumbel, p=.5, marginal=2, x=k, n=1000)
  quantile(x=preds, probs=c(.05, .95), na.rm=TRUE)})
#add these confidence intervals to the plot
arrows(xpoints2, multy[1, ], xpoints2, multy[2, ],
       code=3, length=0.1, angle=90, col="darkgreen")



##Inspect the Monte-Carlo sample used for constructing the confidence
#interval at x=5.01
plot(density(predsY))
abline(v=qrCopula(copula=fitGumbel, p=.5, marginal=2, x=x), col="blue", lty=2)
abline(v=c(qsY), col="red", lty=2)



##Quantile regression: regress x on y
ypoints <- seq(.0001, 25, length.out=100)
#compute the median
xmed <- sapply(ypoints, function (k) qrCopula(copula=fitGumbel, p=.5,
                                              marginal=1, x=k))
#compute the .05th-quantile
x05 <- sapply(ypoints, function (k) qrCopula(copula=fitGumbel, p=.05,
                                             marginal=1, x=k))
#compute the .95th-quantile
x95 <- sapply(ypoints, function (k) qrCopula(copula=fitGumbel, p=.95,
                                             marginal=1, x=k))


#Plot the results
plot(xys[,1], xys[,2], xlab="x", ylab="y", col=gray(.2,.2))
#add the median
lines(xmed, ypoints, lty=1, col="red")
#add the 90% prediction interval: approximately 20 (out of 200) points
#should lie outside the boundaries of the 90% prediction interval
lines(x05, ypoints, lty=2, col="blue")
lines(x95, ypoints, lty=2, col="blue")
#add a legend
legend("topleft", lty=c(1, 2),
       col=c("red", "blue"),
       legend=c("median", "90% prediction interval"),
       bty="n", y.intersp=1.2, cex=.7)


#Use a Monte-Carlo method for computing a confidence interval
#for the median at some given y-value
y <- qrCopula(copula=fitGumbel, p=.5, marginal=2, x=x)
predsX <- qrCopulaMC(copula=fitGumbel, p=.5, marginal=1, x=y, n=5000)
(qsX <- quantile(x=predsX, probs=c(.05, .95), na.rm=TRUE))
#add the confidence interval to the plot
points(qrCopula(copula=fitGumbel, p=.5, marginal=1, x=y), y, col="darkgreen")
arrows(qsX[1], y, qsX[2], y, code=3, length=0.1, angle=90, col="darkgreen")


#Compute the confidence interval for the median at some other points
ypoints2 <- seq(1, 17, length.out=5)
multx <- sapply(ypoints2, function (k) {
  preds <- qrCopulaMC(copula=fitGumbel, p=.5, marginal=1, x=k, n=1000)
  quantile(x=preds, probs=c(.05, .95), na.rm=TRUE)})
#add these confidence intervals to the plot
arrows(multx[1, ], ypoints2, multx[2, ], ypoints2, 
       code=3, length=0.1, angle=90, col="darkgreen")



##Inspect the Monte-Carlo sample used for constructing the confidence interval
plot(density(predsX))
abline(v=qrCopula(copula=fitGumbel, p=.5, marginal=1, x=y), col="blue", lty=2)
abline(v=c(qsX), col="red", lty=2)





##EXAMPLE 2: comparing linear quantile regression and Ordinary Least Squares 

##Generate some artificial data using a Gaussian copula

#Bivariate Gaussian copula
myCop <- normalCopula(param=.6, dispstr="un", dim=2)

#Marginals:
#the first marginal is a normal distribution, and
#the second marginal is also a normal distribution
myMvdc <- mvdc(copula=myCop, c("norm", "norm"),
               list(list(mean=25, sd=3), list(mean=3, sd=2)))

#Draw sample (this sample will subsequently be used as the artificial data)
ns <- 100 #sample size
xys <- rMvdc(ns, myMvdc)

#Inspect data
plot(xys[,1], xys[,2], xlab="x", ylab="y")

#Visualize the copula
persp(myMvdc, dMvdc, xlim=c(15, 35), ylim=c(-5, 10))
contour(myMvdc, dMvdc, xlim=c(15, 35), ylim=c(-5, 10))



##Fit a bivariate Gaussian copula to the artificial data

#Obtain start values for fitting the copula
normMLE <- fitdistr(xys[,1], densfun="normal")
normMean1 <- normMLE$estimate[1]
normSD1 <- normMLE$estimate[2]

normMLE <- fitdistr(xys[,2], densfun="normal")
normMean2 <- normMLE$estimate[1]
normSD2 <- normMLE$estimate[2]

#Starting value for the copula parameter
tauEst <- sin(cor(xys[,1], xys[,2], method="kendall")*pi/2)

#Fit a Gaussian copula to the artificial data
myGaussian <- mvdc(copula=normalCopula(param=1, dispstr="un"), c("norm", "norm"),
                   paramMargins=list(list(mean=0, sd=1), list(mean=0, sd=1)))

start <- c(normMean1, normSD1, normMean2, normSD2, tauEst)

fitGaussian <- fitMvdc(xys, myGaussian, start=start)
fitGaussian



##Quantile regression: regress y on x
xpoints <- seq(10, 40, length.out=100)
#compute the median
ymed <- sapply(xpoints, function (k) qrCopula(copula=fitGaussian, p=.5,
                                              marginal=2, x=k))
#compute the .05th-quantile
y05 <- sapply(xpoints, function (k) qrCopula(copula=fitGaussian, p=.05,
                                             marginal=2, x=k))
#compute the .95th-quantile
y95 <- sapply(xpoints, function (k) qrCopula(copula=fitGaussian, p=.95,
                                             marginal=2, x=k))


#Plot the results
plot(xys[,1], xys[,2], xlab="x", ylab="y",
     xlim=c(15, 35), col=gray(.2,.2))
#add the median
lines(xpoints, ymed, lty=1, col="red")
#add the 90% prediction interval: approximately 10 (out of 100) points
#should lie outside the boundaries of the 90% prediction interval
lines(xpoints, y05, lty=2, col="blue")
lines(xpoints, y95, lty=2, col="blue")
#add a legend
legend("topleft", lty=c(1, 2),
       col=c("red", "blue"),
       legend=c("median", "90% prediction interval"),
       bty="n", y.intersp=1.2, cex=.7)


#Use a Monte-Carlo method for computing a confidence interval
#for the median at some given x-value
x <- 24.5
predsY <- qrCopulaMC(copula=fitGaussian, p=.5, marginal=2, x=x, n=5000)
(qsY <- quantile(x=predsY, probs=c(.05, .95), na.rm=TRUE))
#add the confidence interval to the plot
points(x, qrCopula(copula=fitGaussian, p=.5, marginal=2, x=x), col="darkgreen")
arrows(x, qsY[1], x , qsY[2], code=3, length=0.1, angle=90, col="darkgreen")


#Compute the confidence interval for the median at some other points
xpoints2 <- seq(16, 33, length.out=4)
multy <- sapply(xpoints2, function (k) {
  preds <- qrCopulaMC(copula=fitGaussian, p=.5, marginal=2, x=k, n=5000)
  quantile(x=preds, probs=c(.05, .95), na.rm=TRUE)})
#add these confidence intervals to the plot
arrows(xpoints2, multy[1, ], xpoints2, multy[2, ],
       code=3, length=0.1, angle=90, col="darkgreen")



##Fit a linear model (using OLS) to the data
regline <- lm(y ~ x, data=data.frame(y=xys[,2], x=xys[,1]))
#Predict the conditional means and compute their 90% confidence intervals
xpoints3 <- c(xpoints2[1:2], x, xpoints2[3:4])
preds <- predict(regline, newdata=data.frame(x=xpoints3), 
                 interval="confidence", level=.90)
#add the regression line to the plot
abline(reg=regline, col="darkgreen", lty=2, lwd=2)
#notice that the fitted line for the median (obtained with quantile regression)
#coincides with the fitted regression line (as obtained by OLS)

#add the confidence interval for the predicted conditional means to the plot
lines(xpoints3, preds[ ,3], lty=2, col="darkgreen")
lines(xpoints3, preds[ ,2], lty=2, col="darkgreen")
#notice that the widths of the Monte-Carlo and OLS
#confidence intervals are similar



##Inspect the Monte-Carlo sample used for constructing the confidence interval
#at x=24.5
plot(density(predsY))
abline(v=qrCopula(copula=fitGaussian, p=.5, marginal=2, x=x), col="blue", lty=2)
abline(v=c(qsY), col="red", lty=2)





##EXAMPLE 3: Koenker's copula method

#In the vignette of the quantreg package, Koenker implements a nonlinear
#quantile regression model using a Frank copula.
#(See section "Nonlinear quantile regression" at:
#https://cran.r-project.org/web/packages/quantreg/vignettes/rq.pdf)
#Readers of the vignette may have noticed some differences between the copula
#method as implemented by Koenker in his vignette and the implementation
#of the copula method in the qrCopula function. For instance, Koenker's
#implementation in the vignette allows the copula parameter (denoted by delta)
#to vary across different p values (where Koenker denotes p by tau).
#However, in the qrCopula function the copula parameter is kept constant
#across all p values.
#The following example will (hopefully) demonstrate what the implications
#are of either varying or keeping the copula parameter constant across
#different p values.

library(quantreg)

#Generate artificial data in way that is identical to Koenker in the vignette
df <- 8
delta <- 8

#Define a Frank copula with t-distributions as marginals
myMvdc <- mvdc(copula=frankCopula(param=8, dim=2), c("t", "t"),
               list(list(df=8), list(df=8)))

ns <- 200 #sample size
xys <- rMvdc(ns, myMvdc) #artificial data sample


#Plot the data
plot(xys[,1], xys[,2], xlab="x", ylab="y", col="blue", cex=.25)


#Fit a bivariate Frank copula to the data
myFrank <- mvdc(copula=frankCopula(param=1, dim=2), c("t", "t"),
                list(list(df=1), list(df=1)))

#Specify some naive start values for fitting the Frank copula
start <- c(5, 5, 5)

fitFrank <- fitMvdc(xys, myFrank, start=start)
fitFrank
#Note that the Frank copula also estimated the degrees of freedom
#for the two marginal t-distributions. 


#Fit a nonlinear quantile regression model using the qrCopula function.
#Remember that the qrCopula function keeps the copula parameter constant
#across different p values.
xpoints <- seq(-8, 8, length.out=200)
ymed <- sapply(xpoints, function (k) qrCopula(copula=fitFrank, p=.5,
                                              marginal=2, x=k))
y05 <- sapply(xpoints, function (k) qrCopula(copula=fitFrank, p=.05,
                                             marginal=2, x=k))
y95 <- sapply(xpoints, function (k) qrCopula(copula=fitFrank, p=.95,
                                             marginal=2, x=k))

#add the estimated quantile curves to the plot
lines(xpoints, ymed, lty=2, col="blue")
lines(xpoints, y05, lty=2, col="blue")
lines(xpoints, y95, lty=2, col="blue")


#Compute the true quantile curves.
#Note that these quantile curves will be identical to the true quantile
#curves as computed by Koenker in the vignette.
fitFrank2 <- fitFrank
fitFrank2@mvdc@copula@parameters <- 8
fitFrank2@mvdc@paramMargins[[1]]$df <- 8
fitFrank2@mvdc@paramMargins[[2]]$df <- 8

ymed <- sapply(xpoints, function (k) qrCopula(copula=fitFrank2, p=.5,
                                              marginal=2, x=k))
y05 <- sapply(xpoints, function (k) qrCopula(copula=fitFrank2, p=.05,
                                             marginal=2, x=k))
y95 <- sapply(xpoints, function (k) qrCopula(copula=fitFrank2, p=.95,
                                             marginal=2, x=k))

#add the true quantile curves to the plot
lines(xpoints, ymed, lty=1, col="black")
lines(xpoints, y05, lty=1, col="black")
lines(xpoints, y95, lty=1, col="black")


#Koenker model 1:
#This is the model as fitted by Koenker in the vignette. Notice that
#Koenker has a different copula parameter estimated at each
#p-value (i.e., the copula parameter is allowed to vary across different
#p-values), but he fixes the degrees of freedom at df=8.
Dat <- NULL
Dat$x <- x <- xys[order(xys[,1]) ,1]
Dat$y <- y <- xys[order(xys[,1]) ,2]
us <- c(.05, .5, .95)
deltas <- matrix(0, 3, length(us))
colnames(deltas) <- c("Copula parameter", "Mu t-distribution",
                      "Sigma t-distribution")

#In the following function, Koenker implements steps 1-4 of Nelsen's
#method for fitting a (nonlinear) quantile regression model.
FrankModel <- function(x, delta, mu, sigma, df, tau){
  z <- qt(-log(1-(1-exp(-delta))/
                 (1+exp(-delta*pt(x, df))*((1/tau)-1)))/delta, df)
  mu + sigma*z}

#Koenker fits the quantile regression model by estimating the parameters
#delta (=copula parameter), mu, and sigma (where mu and sigma are, respectively,
#the location and scale of the t-distribution). Koenker fits these
#parameters at each specified value of p (=tau). In that way, a parameter
#such as delta (=copula parameter) is allowed to vary across different
#values for p.
for(i in 1:length(us)){
  tau=us[i]
  fit <- nlrq(y ~ FrankModel(x, delta ,mu, sigma, df=8, tau=tau),
              data=Dat, tau=tau, start=list(delta=5,
                                            mu=0, sigma=1), trace=FALSE)
  lines(x, predict(fit, newdata=x), lty=2, col="green")
  deltas[i,] <- coef(fit)}
#Note: for some artificial data samples this model will yield an error.
#If such an error might occur, then proceed with "Koenker model 2" below,
#which is usually more stable.

#The fitted copula parameter indeed varied across the different p values
cbind(us, deltas) 


#Koenker model 2:
#This is a modification of "Koenker model 1", where mu and sigma of the
#t-distribution are not estimated anymore (but fixed at 0 and 1, respectively).
#This modified model should yield a better fit than "Koenker model 1". 
FrankModel2<- function(x, delta, df1, df2, tau){
  z <- qt(-log(1-(1-exp(-delta))/
                 (1+exp(-delta*pt(x, df1))*((1/tau)-1)))/delta, df2)
  z}

deltas <- matrix(0,length(us), 1)

for(i in 1:length(us)){
  tau=us[i]
  fit <- nlrq(y~FrankModel2(x, delta, df1=8, df2=8, tau=tau),
              data=Dat, tau=tau, start=list(delta=5), trace=FALSE)
  lines(x, predict(fit, newdata=x), lty=2, col="darkgreen")
  deltas[i,] <- coef(fit)}

#Again, the fitted copula parameter varied across the different p values
cbind(us, deltas)


#As a final comparison, we will fit the Frank copula again, but similar to
#"Koenker model 2" we will fix df=8, mu=0, and sigma=1 for the
#two marginal t-distributions.
u1 <- pt(x, df=8)
u2 <- pt(y, df=8)

paramCop <- fitCopula(frankCopula(param=5, dim=2), cbind(u1, u2))

fitFrank3 <- fitFrank

fitFrank3@mvdc@copula@parameters <- paramCop@estimate
fitFrank3@mvdc@paramMargins[[1]]$df <- 8
fitFrank3@mvdc@paramMargins[[2]]$df <- 8


#We finally fit a nonlinear quantile regression model using the qrCopula function.
#The input for the model will be the Frank copula above, where the marginal
#t-distributions were fixed at df=8, mu=0, and sigma=1.
#This final model will be closest to "Koenker model 2", except that
#"Koenker model 2" allows the copula parameter to vary, whereas this model
#keeps the copula constant across different p values.
ymed <- sapply(xpoints, function (k) qrCopula(copula=fitFrank3, p=.5,
                                              marginal=2, x=k))
y05 <- sapply(xpoints, function (k) qrCopula(copula=fitFrank3, p=.05,
                                             marginal=2, x=k))
y95 <- sapply(xpoints, function (k) qrCopula(copula=fitFrank3, p=.95,
                                             marginal=2, x=k))

#add the estimated quantile curves to the plot
lines(xpoints, ymed, lty=2, col="purple")
lines(xpoints, y05, lty=2, col="purple")
lines(xpoints, y95, lty=2, col="purple")
#finally, add a legend to the plot
legend("topleft", lty=c(1, 2, 2, 2, 2),
       col=c("black", "blue", "purple", "green", "darkgreen"),
       legend=c("True quantile curve", "qrCopula, df free",
                "qrCopula, df fixed", "Koenker model 1", "Koenker model 2"),
       bty="n", y.intersp=1.2, cex=.7)

#If you repeat the above procedure for fitting qrCopula (constant copula
#parameter) and the "Koenker models" (varying copula parameter) several times,
#but each time for a newly generated artificial data sample, you will get a sense
#of which of these two implementations (constant vs. varying copula parameter)
#fits quantile curves that lie the closest to the true quantile curves.