R/qfa5.0.R

Defines functions sqr_deriv.plot sqr.plot boot.sqr sqdft sqdft.fit tsqr.fit sqr.fit sqr1.fit sqr3.fit sqr.fit.optim sqr.optim_dobj sqr.optim_obj optim.grad optim.adam qspec.lw qspec.ar qser2ar Asmooth Vsmooth qser2qacf qcser per sar.eq.test sar.eq.bootstrap sar.gc.test sar.gc.bootstrap sar.gc.param.validity sar.residual sar.gc.coef qper2 qspec.sar qser2sar ar2qspec qper qser qacf rq.fit.fnb2 qkl.divergence qsmooth.qper qsmooth qspec2qcoh qacf2speclw qdft2qser qdft2qacf qdft2qper qfa.plot qdft tqr.fit

Documented in boot.sqr per qacf qcser qdft qdft2qacf qdft2qper qdft2qser qfa.plot qkl.divergence qper qper2 qser qser2ar qser2qacf qser2sar qspec2qcoh qspec.ar qspec.lw qspec.sar sar.eq.bootstrap sar.eq.test sar.gc.bootstrap sar.gc.coef sar.gc.test sqdft sqdft.fit sqr1.fit sqr3.fit sqr_deriv.plot sqr.fit sqr.fit.optim sqr.plot tqr.fit tsqr.fit

# -- Library of R functions for Quantile-Frequency Analysis (QFA) --
# by Ta-Hsin Li  (thl024@outlook.com)  October 30, 2022; April 8, 2023; August 21, 2023;
# November 26, 2024; December 14, 2024; April 8, 2025; September 10, 2025; March 23, 2026



#' Trigonometric Quantile Regression (TQR)
#'
#' This function computes ordinary trigonometric quantile regression (TQR) for univariate time series at a single frequency.
#' @param y vector of time series
#' @param f0 frequency in [0,1)
#' @param tau sequence of quantile levels in (0,1)
#' @param prepared if \code{TRUE}, intercept is removed and coef of cosine is doubled when \code{f0 = 0.5}
#' @return object of \code{quantreg::rq()} 
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' fit <- tqr.fit(y,f0=0.1,tau=tau)
#' plot(tau,fit$coef[1,],type='o',pch=0.75,xlab='QUANTILE LEVEL',ylab='TQR COEF')
tqr.fit <- function(y,f0,tau,prepared=TRUE) {

  tqr_fix_coef <- function(coef) {
  # prepare coef from tqr for qdft
  # input:  coef = p*ntau tqr coefficient matrix from tqr.fit()
  # output: 2*ntau matrix of tqr coefficients
    ntau <- ncol(coef)
    if(nrow(coef)==1) {   
      # for f = 0
      coef <- rbind(rep(0,ntau),rep(0,ntau))
    } else if(nrow(coef)==2) {  
      # for f = 0.5: rescale coef of cosine by 2 so qdft can be defined as usual
      coef <- rbind(2*coef[2,],rep(0,ntau))
    } else {
      # for f in (0,0.5)
      coef <- coef[-1,]
    }
    coef
  }

  n <- length(y)
  tt <- c(1:n)
  if(f0 != 0.5 & f0 != 0) {
    fit <- suppressWarnings(quantreg::rq(y ~ cos(2*pi*f0*tt)+sin(2*pi*f0*tt),tau=tau))
  }
  if(f0 == 0.5) {
    fit <- suppressWarnings(quantreg::rq(y ~ cos(2*pi*f0*tt),tau=tau))
  }
  if(f0 == 0) {
    fit <- suppressWarnings(quantreg::rq(y ~ 1,tau=tau))
  }
  if(prepared) fit$coefficients <- tqr_fix_coef(fit$coefficients)
  fit
}


#' Quantile Discrete Fourier Transform (QDFT)
#'
#' This function computes quantile discrete Fourier transform (QDFT) for univariate or multivariate time series.
#' @param y vector or matrix of time series (if matrix, \code{nrow(y)} = length of time series)
#' @param tau sequence of quantile levels in (0,1)
#' @param n.cores number of cores for parallel computing (default = 1)
#' @param cl pre-existing cluster for repeated parallel computing (default = \code{NULL})
#' @return matrix or array of quantile discrete Fourier transform of \code{y}
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y.qdft <- qdft(y,tau)
#' # Make a cluster for repeated use
#' n.cores <- 2
#' cl <- parallel::makeCluster(n.cores)
#' parallel::clusterExport(cl, c("tqr.fit"))
#' doParallel::registerDoParallel(cl)
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y.qdft <- qdft(y1,tau,n.cores=n.cores,cl=cl)
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y.qdft <- qdft(y2,tau,n.cores=n.cores,cl=cl)
#' parallel::stopCluster(cl)
qdft <- function(y,tau,n.cores=1,cl=NULL) {

  z <- function(x) { x <- matrix(x,nrow=2); x[1,]-1i*x[2,] }

  extend.qdft <- function(y,tau,result,sel.f) {
  # define qdft from tqr coefficients
  # input: y = ns*nc-matrix or ns-vector of time series
  #        tau = ntau-vector of quantile levels
  #        result = list of qdft in (0,0.5]
  #        sel.f = index of Fourier frequencies in (0,0.5)
  # output: out = nc*ns*ntau-array or ns*ntau matrix of qdft

    if(!is.matrix(y)) y <- matrix(y,ncol=1)
    nc <- ncol(y)
    ns <- nrow(y)
    ntau <- length(tau)
  
    out <- array(NA,dim=c(nc,ns,ntau))
    for(k in c(1:nc)) {
      # define QDFT at freq 0 as ns * quantile
      out[k,1,] <- ns * quantile(y[,k],tau)
      # retrieve qdft for freq in (0,0.5]
      tmp <- matrix(unlist(result[[k]]),ncol=ntau,byrow=TRUE)
      # define remaining values by conjate symmetry (excluding f=0.5) 
      tmp2 <- NULL
      for(j in c(1:ntau)) {
        tmp2 <- cbind(tmp2,rev(Conj(tmp[sel.f,j])))
      }
      # assemble & rescale everything by ns/2 so that periodogram = |dft|^2/ns
      out[k,c(2:ns),] <- rbind(tmp,tmp2) * ns/2
    }
    if(nc == 1) out <- matrix(out[1,,],ncol=ntau)
    out
  }

  if(!is.matrix(y)) y <- matrix(y,ncol=1)
  ns <- nrow(y)
  nc <- ncol(y)
  f2 <- c(0:(ns-1))/ns
  # Fourier frequencies in (0,0.5]
  f <- f2[which(f2 > 0 & f2 <= 0.5)]
  sel.f <- which(f < 0.5)
  nf <- length(f)
  ntau <- length(tau)
  keep.cl <- TRUE
  if(n.cores>1 & is.null(cl)) {
    cl <- parallel::makeCluster(n.cores)
    parallel::clusterExport(cl, c("tqr.fit"))
    doParallel::registerDoParallel(cl)
    keep.cl <- FALSE
  }
  `%dopar%` <- foreach::`%dopar%`
  `%do%` <- foreach::`%do%`
  # compute qdft for f in (0,0.5]
  result <- list()
  i <- 0
  for(k in c(1:nc)) {
    yy <- y[,k] 
    if(n.cores>1) {
      tmp <- foreach::foreach(i=1:nf) %dopar% {
	    tqr.fit(yy,f[i],tau)$coef
      }
    } else {
      tmp <- foreach::foreach(i=1:nf) %do% {
        tqr.fit(yy,f[i],tau)$coef
      }
    }
    # tmp = a list over freq of 2 x ntau coefficiets 
    tmp <- lapply(tmp,FUN=function(x) {apply(x,2,z)})
    result[[k]] <- tmp
  }
  if(n.cores>1 & !keep.cl) {
    parallel::stopCluster(cl)
    cl <-NULL
  }
  # extend qdft to f=0 and f in (0.5,1) 
  out <- extend.qdft(y,tau,result,sel.f)
  return(out) 
}



#' Quantile-Frequency Plot
#'
#' This function creates an image plot of quantile spectrum.
#' @param freq sequence of frequencies in (0,0.5) at which quantile spectrum is evaluated
#' @param tau sequence of quantile levels in (0,1) at which quantile spectrum is evaluated
#' @param rqper real-valued matrix of quantile spectrum evaluated on the \code{freq} x \code{tau} grid
#' @param rg.qper \code{zlim} for \code{qper} (default = \code{range(qper)})
#' @param rg.tau  \code{ylim} for \code{tau} (default = \code{range(tau)})
#' @param rg.freq \code{xlim} for \code{freq} (default = \code{c(0,0.5)})
#' @param color colors (default = \code{colorRamps::matlab.like2(1024)})
#' @param ylab label of y-axis (default = \code{"QUANTILE LEVEL"})
#' @param xlab label of x-axis (default = \code{"FREQUENCY"})
#' @param tlab title of plot (default = \code{NULL})
#' @param set.par if \code{TRUE}, \code{par()} is set internally (single image)
#' @param legend.plot if \code{TRUE}, legend plot is added
#' @param line distance between the image and the axis (default = 0.5)
#' @return no return value
#' @import 'fields'
#' @import 'graphics'
#' @import 'colorRamps'
#' @export
qfa.plot <- function(freq,tau,rqper,rg.qper=range(rqper),rg.tau=range(tau),rg.freq=c(0,0.5),
  color=colorRamps::matlab.like2(1024),ylab="QUANTILE LEVEL",xlab="FREQUENCY",tlab=NULL,
  set.par=TRUE,legend.plot=TRUE,line=0.5) {

  if(set.par) { 
    oldpar <- par(no.readonly = TRUE) 
    on.exit(par(oldpar))
    graphics::par(mfrow=c(1,1),pty="m",lab=c(7,7,7),mar=c(4,4,3,6)+0.1,las=1)
  }
	
  tmp<-rqper
  tmp[tmp<rg.qper[1]]<-rg.qper[1]
  tmp[tmp>rg.qper[2]]<-rg.qper[2]
  graphics::image(freq,tau,tmp,xlab=xlab,ylab=ylab,xlim=rg.freq,ylim=rg.tau,col=color,zlim=rg.qper,frame=F,axes=F)
  graphics::axis(1,line=line)
  graphics::axis(2,line=line,at=seq(0,1,0.1))
  if(line==0) { abline(h=max(tau)); abline(v=0.5) }
  if(legend.plot) {
    fields::image.plot(zlim=rg.qper,legend.only=TRUE,nlevel = 32,legend.mar = NULL, legend.shrink = 1,col=color)
  }
  if(!is.null(tlab)) graphics::title(tlab)
}


#' Quantile Periodogram (QPER)
#'
#' This function computes quantile periodogram (QPER) from QDFT.
#' @param y.qdft matrix or array of QDFT from \code{qdft()}
#' @return matrix or array of quantile periodogram
#' @export
#' @examples
#' # single time series
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y.qdft <- qdft(y1,tau)
#' y.qper <- qdft2qper(y.qdft)
#' n <- length(y1)
#' ff <- c(0:(n-1))/n
#' sel.f <- which(ff > 0 & ff < 0.5)
#' qfa.plot(ff[sel.f],tau,Re(y.qper[sel.f,]))
#' # multiple time series
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' y.qdft <- qdft(cbind(y1,y2),tau)
#' y.qper <- qdft2qper(y.qdft)
#' qfa.plot(ff[sel.f],tau,Re(y.qper[1,1,sel.f,]))
#' qfa.plot(ff[sel.f],tau,Re(y.qper[1,2,sel.f,]))
qdft2qper <- function(y.qdft) {

  if(is.matrix(y.qdft)) {
    y.qdft <- array(y.qdft,dim=c(1,dim(y.qdft)[1],dim(y.qdft)[2]))
  }
  nc <- dim(y.qdft)[1]
  ns  <- dim(y.qdft)[2]
  nf  <- ns
  ntau <- dim(y.qdft)[3]
  y.qper <- array(NA,dim=c(nc,nc,nf,ntau))
  for(k in c(1:nc)) {
     tmp1 <- matrix(y.qdft[k,,],ncol=ntau)
     y.qper[k,k,,] <- (Mod(tmp1)^2) / ns
     if(k < nc) {
       for(kk in c((k+1):nc)) {
            tmp2 <- matrix(y.qdft[kk,,],ncol=ntau)
            y.qper[k,kk,,] <- tmp1 * Conj(tmp2) / ns
            y.qper[kk,k,,] <- Conj(y.qper[k,kk,,])
       }
     }
  }
  if(nc==1) y.qper <- matrix(y.qper[1,1,,],nrow=nf,ncol=ntau)
  return(y.qper)
}



#' Quantile Autocovariance Function (QACF)
#'
#' This function computes quantile autocovariance function (QACF) from QDFT.
#' @param y.qdft matrix or array of QDFT from \code{qdft()}
#' @param return.qser if \code{TRUE}, return quantile series (QSER) along with QACF
#' @return matrix or array of quantile autocovariance function if \code{return.sqer = FALSE} (default), else a list with the following elements:
#'   \item{qacf}{matirx or array of quantile autocovariance function}
#'   \item{qser}{matrix or array of quantile series}
#' @import 'stats'
#' @export
#' @examples
#' # single time series
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y.qdft <- qdft(y1,tau)
#' y.qacf <- qdft2qacf(y.qdft)
#' plot(c(0:9),y.qacf[c(1:10),1],type='h',xlab="LAG",ylab="QACF")
#' y.qser <- qdft2qacf(y.qdft,return.qser=TRUE)$qser
#' plot(y.qser[,1],type='l',xlab="TIME",ylab="QSER")
#' # multiple time series
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' y.qdft <- qdft(cbind(y1,y2),tau)
#' y.qacf <- qdft2qacf(y.qdft)
#' plot(c(0:9),y.qacf[1,2,c(1:10),1],type='h',xlab="LAG",ylab="QACF")
qdft2qacf <- function(y.qdft,return.qser=FALSE) {

  if(is.matrix(y.qdft)) {
    y.qdft <- array(y.qdft,dim=c(1,dim(y.qdft)[1],dim(y.qdft)[2]))
  }
  nc <- dim(y.qdft)[1]
  ns  <- dim(y.qdft)[2]
  nf  <- ns
  ntau <- dim(y.qdft)[3]
  yy <- array(NA,dim=dim(y.qdft))
  for(k in c(1:nc)) {
    yy[k,,] <- Re(matrix(apply(matrix(y.qdft[k,,],ncol=ntau),2,stats::fft,inverse=TRUE),ncol=ntau)) / ns
  }
  y.qacf <- array(NA,dim=c(nc,nc,dim(y.qdft)[-1]))
  if(nc > 1) {
    for(j in c(1:ntau)) {
      # demean = TRUE is needed
      tmp <- stats::acf(t(as.matrix(yy[,,j],ncol=nc)), type = "covariance", lag.max = ns-1, plot = FALSE, demean = TRUE)$acf
      for(k in c(1:nc)) {
        for(kk in c(1:nc)) {
          y.qacf[k,kk,,j] <- tmp[,k,kk] 
	    }
      }
    }
    rm(tmp)
  } else {
    for(j in c(1:ntau)) {
      # demean = TRUE is needed
      tmp <- stats::acf(c(yy[,,j]), type = "covariance", lag.max = ns-1, plot = FALSE, demean = TRUE)$acf 
      y.qacf[1,1,,j] <- tmp[,1,1]
    }	 
  }
  if(nc==1) {
    y.qacf <- matrix(y.qacf[1,1,,],ncol=ntau)
    yy <- matrix(yy[1,,],ncol=ntau)
  }
  if(return.qser) {
    return(list(qacf=y.qacf,qser=yy))
  } else {
    return(y.qacf)
  }
}



#' Quantile Series (QSER)
#'
#' This function computes quantile series (QSER) from QDFT.
#' @param y.qdft matrix or array of QDFT from \code{qdft()}
#' @return matrix or array of quantile series
#' @import 'stats'
#' @export
#' @examples
#' # single time series
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y.qdft <- qdft(y1,tau)
#' y.qser <- qdft2qser(y.qdft)
#' plot(y.qser[,1],type='l',xlab="TIME",ylab="QSER")
#' # multiple time series
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' y.qdft <- qdft(cbind(y1,y2),tau)
#' y.qser <- qdft2qser(y.qdft)
#' plot(y.qser[1,,1],type='l',xlab="TIME",ylab="QSER")
qdft2qser <- function(y.qdft) {

  if(is.matrix(y.qdft)) {
    y.qdft <- array(y.qdft,dim=c(1,dim(y.qdft)[1],dim(y.qdft)[2]))
  }
  nc <- dim(y.qdft)[1]
  ns  <- dim(y.qdft)[2]
  nf  <- ns
  ntau <- dim(y.qdft)[3]
  y.qser <- array(NA,dim=dim(y.qdft))
  for(k in c(1:nc)) {
    y.qser[k,,] <- Re(matrix(apply(matrix(y.qdft[k,,],ncol=ntau),2,stats::fft,inverse=TRUE),ncol=ntau)) / ns
  }
  if(nc==1) {
    y.qser <- matrix(y.qser[1,,],ncol=ntau)
  }
  return(y.qser)
}


# Lag-Window (LW) Estimator of Quantile Spectrum
#
# This function computes lag-window (LW) estimate of quantile spectrum from QACF.
# @param y.qacf matrix or array of pre-calculated QACF from \code{qdft2qacf()}
# @param M bandwidth parameter of lag window (default = \code{NULL}: quantile periodogram)
# @return A list with the following elements:
#   \item{spec}{matrix or array of quantile spectrum}
#   \item{lw}{lag-window sequence}
# @import 'stats'
qacf2speclw <- function(y.qacf,M=NULL) {
  
  if(is.matrix(y.qacf)) y.qacf <- array(y.qacf,dim=c(1,1,dim(y.qacf)))
  
  nc <- dim(y.qacf)[1]
  ns <- dim(y.qacf)[3]
  ntau <- dim(y.qacf)[4]
  nf <- ns
  nf2 <- 2*ns
    
  Hanning <- function(n,M) {
    # Hanning window
    lags<-c(0:(n-1))
    tmp<-0.5*(1+cos(pi*lags/M))
    tmp[lags>M]<-0
    tmp<-tmp+c(0,rev(c(tmp[-1])))
    tmp
  }

  if(is.null(M)) {
    lw <- rep(1,nf2)
  } else {
    lw<-Hanning(nf2,M)
  }
  qper.lw <- array(NA,dim=c(nc,nc,nf,ntau))
  for(k in c(1:nc)) {
      for(j in c(1:ntau)) {
	    gam <- c(y.qacf[k,k,,j],0,rev(y.qacf[k,k,-1,j]))*lw
  	    qper.lw[k,k,,j] <- Mod(stats::fft(gam,inverse=FALSE))[seq(1,nf2,2)]
	  }
      if(k < nc) {
        for(kk in c((k+1):nc)) {
          for(j in c(1:ntau)) {
		    gam <- c(y.qacf[k,kk,,j],0,rev(y.qacf[kk,k,-1,j]))*lw
  	        qper.lw[k,kk,,j] <- stats::fft(gam,inverse=FALSE)[seq(1,nf2,2)]
		  }
          qper.lw[kk,k,,] <- Conj(qper.lw[k,kk,,])
        }
	  }
  }
  
  # return spectrum
  if(nc==1) qper.lw <- matrix(qper.lw[1,1,,],ncol=ntau)
  return(out=list(spec=qper.lw,lw=lw))
}




#' Quantile Coherence Spectrum
#'
#' This function computes quantile coherence spectrum (QCOH) from quantile spectrum of multiple time series.
#' @param qspec array of quantile spectrum 
#' @param k index of first series (default = 1)
#' @param kk index of second series (default = 2)
#' @return matrix of quantile coherence evaluated at Fourier frequencies in (0,0.5)
#' @export
#' @examples
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' n <- length(y1)
#' ff <- c(0:(n-1))/n
#' sel.f <- which(ff > 0 & ff < 0.5)
#' y.qacf <- qacf(cbind(y1,y2),tau)
#' y.qper.lw <- qspec.lw(y.qacf=y.qacf,M=5)$spec
#' y.qcoh <- qspec2qcoh(y.qper.lw,k=1,kk=2)
#' qfa.plot(ff[sel.f],tau,y.qcoh)
qspec2qcoh<-function(qspec,k=1,kk=2) {

  nf <- dim(qspec)[3]
  ff <- c(0:(nf-1))/nf
  sel.f <- which(ff > 0 & ff < 0.5)
  coh <- Mod(qspec[k,kk,sel.f,])^2/(Re(qspec[k,k,sel.f,])*Re(qspec[kk,kk,sel.f,]))
  coh[coh > 1] <- 1
  coh
}


qsmooth <- function(x,link=c('linear','log'),method=c('gamm','sp'),spar=NULL) {
# smooth a vector of real or complex data
# input:   x = vector of real or complex data to be smoothed
#          link = in linear or log domain for smoothing
#          method = gamm: spline smoothing using gamm()
#                   sp:   spline smoothing using smooth.spline()
#          spar = smoothing parameter for smooth.spline() (default: spar=NULL for GCV)
# output:  tmp = smoothed data
# imports: gamm() in 'mgcv'
#          corAR1() in 'nlme' (included in 'mgcv')
  
  smooth.spline.gamm <- function(x,y=NULL) {
  # smoothing splines with correlated data using gamm()
  # input:        x = vector of response (if y=NULL) or independent variable
  #               y = vector of response if x is independent variable
    if(is.null(y)) {
      y <- x
	  x <- c(1:length(y)) / length(y)
    }
    # fitting spline with AR(1) residuals
    fit <- mgcv::gamm(y ~ s(x), data=data.frame(x=x,y=y), correlation=nlme::corAR1())$gam 
    return(list(x=x,y=stats::fitted(fit),fit=fit))
  }

  eps <- 1e-16
  if(is.complex(x)) {
  	tmp<-Re(x)
    if(abs(diff(range(tmp))) > 0) {
      if(method[1] == 'gamm') {
        tmp<-smooth.spline.gamm(tmp)$y
      } else {
        tmp<-stats::smooth.spline(tmp,spar=spar)$y
      }
	}
  	tmp2<-Im(x)	
    if(abs(diff(range(tmp2))) > 0) {
      if(method[1] == 'gamm') {
        tmp2<-smooth.spline.gamm(tmp2)$y
      } else {
        tmp2<-stats::smooth.spline(tmp2,spar=spar)$y
      }
	}
    tmp <- tmp + 1i*tmp2
  } else {
  	tmp<-x
    if(abs(diff(range(tmp))) > 0) {
      if(link[1]=='log') {
	    sig<-mean(tmp)
	    tmp[tmp<eps*sig] <- eps*sig
        if(method[1] == 'gamm') {
          tmp<-exp(smooth.spline.gamm(log(tmp))$y)	
		} else {
          tmp<-exp(stats::smooth.spline(log(tmp),spar=spar)$y)
		}
      } else {
        if(method[1] == 'gamm') {
          tmp<-smooth.spline.gamm(tmp)$y		
		} else {
          tmp<-stats::smooth.spline(tmp,spar=spar)$y
		}
  	  }
	}
  }
  tmp
}



# Quantile Smoothing of Quantile Periodogram or Spectral Estimate
#
# This function computes quantile-smoothed version of quantile periodogram or other quantile spectral estimate.
# @param y.qper matrix or array of quantile periodogram or spectral estimate
# @param method smoothing method: \code{"gamm"} for \code{mgcv::gamm()} (default), 
# \code{"sp"} for \code{stats::smooth.spline()}
# @param spar smoothing parameter in \code{smooth.spline()} if \code{method = "sp"} (default = NULL for GCV)
# @param n.cores number of cores for parallel computing (default = 1)
# @param cl pre-existing cluster for repeated parallel computing (default = \code{NULL})
# @return matrix or array of quantile-smoothed quantile spectral estimate
# @import 'foreach'
# @import 'parallel'
# @import 'doParallel'
# @import 'mgcv'
# @import 'nlme'
# @export
# @examples
# y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
# y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
# tau <- seq(0.1,0.9,0.05)
# n <- length(y1)
# ff <- c(0:(n-1))/n
# sel.f <- which(ff > 0 & ff < 0.5)
# y.qdft <- qdft(cbind(y1,y2),tau)
# y.qacf <- qdft2qacf(y.qdft)
# y.qper.lw <- qspec.lw(y.qacf=y.qacf,M=5)$spec
# qfa.plot(ff[sel.f],tau,Re(y.qper.lw[1,1,sel.f,]))
# y.qper.lwqs <- qsmooth.qper(y.qper.lw,method="sp",spar=0.9)
# qfa.plot(ff[sel.f],tau,Re(y.qper.lwqs[1,1,sel.f,]))
qsmooth.qper <- function(y.qper,method=c("gamm","sp"),spar=NULL,n.cores=1,cl=NULL) {
  
  if(is.matrix(y.qper)) y.qper <- array(y.qper,dim=c(1,1,dim(y.qper)))
  nc<-dim(y.qper)[1]
  nf<-dim(y.qper)[3]  
  ntau<-dim(y.qper)[4]  
  eps <- 1e-16
  method <- method[1]
  f2 <- c(0:(nf-1))/nf
  # half of the Fourier frequencies excluding 0 for smoothing
  sel.f1 <- which(f2 >= 0 & f2 <= 0.5)
  # half Fourier frequencies excluding 0 and 0.5 for freq folding
  sel.f2 <- which(f2 > 0 & f2 < 0.5)  
  sel.f3 <- which(f2 > 0.5)
  keep.cl <- TRUE
  if(n.cores>1 & is.null(cl)) {
    cl <- parallel::makeCluster(n.cores)
    parallel::clusterExport(cl, c("qsmooth"))
    doParallel::registerDoParallel(cl)
    keep.cl <- FALSE
  }
  spar0 <- spar
  spars <- NULL
  qper.sm <-y.qper
  for(k in c(1:nc)) {
    spars <- rbind(spars,c(k,k,spar0))
	if(!is.null(cl)) {
	  qper.sm[k,k,sel.f1,] <- t(parallel::parApply(cl,Re(y.qper[k,k,sel.f1,]),1,qsmooth,link="log",method=method,spar=spar0))
	} else {
	  qper.sm[k,k,sel.f1,] <- t(apply(Re(y.qper[k,k,sel.f1,]),1,qsmooth,link="log",method=method,spar=spar0))
	}
    for(j in c(1:ntau)) {
	  qper.sm[k,k,sel.f3,j] <- rev(qper.sm[k,k,sel.f2,j])
	}    	
    if(k < nc) {
      for(kk in c((k+1):nc)) {
        spars <- rbind(spars,c(k,kk,spar0))
        if(!is.null(cl)) {
          tmp1 <- t(parallel::parApply(cl,Re(y.qper[k,kk,,]),1,qsmooth,method=method,spar=spar0))
          tmp2 <- t(parallel::parApply(cl,Im(y.qper[k,kk,,]),1,qsmooth,method=method,spar=spar0))
        } else {
          tmp1 <- t(apply(Re(y.qper[k,kk,,]),1,qsmooth,method=method,spar=spar0))
          tmp2 <- t(apply(Im(y.qper[k,kk,,]),1,qsmooth,method=method,spar=spar0))
        }
  	    qper.sm[k,kk,,] <- tmp1 + 1i*tmp2
        qper.sm[kk,k,,] <- Conj(qper.sm[k,kk,,])
	  }
	}
  }
  if(n.cores>1 & !keep.cl) {
    parallel::stopCluster(cl)
    cl <- NULL
  }
  if(nc==1) qper.sm <- matrix(qper.sm[1,1,,],ncol=ntau)
  return(qper.sm)
}



#' Kullback-Leibler Divergence of Quantile Spectral Estimate
#'
#' This function computes Kullback-Leibler divergence (KLD) of quantile spectral estimate.
#' @param y.qper matrix or array of quantile spectral estimate from, e.g., \code{qspec.lw()}
#' @param qspec matrix of array of true quantile spectrum (same dimension as \code{y.qper})
#' @param sel.f index of selected frequencies for computation (default = \code{NULL}: all frequencies)
#' @param sel.tau index of selected quantile levels for computation (default = \code{NULL}: all quantile levels)
#' @return real number of Kullback-Leibler divergence
#' @export
qkl.divergence <- function(y.qper,qspec,sel.f=NULL,sel.tau=NULL) {

  if(is.matrix(y.qper)) {
      y.qper <- array(y.qper,dim=c(1,1,dim(y.qper)))
      qspec <- array(qspec,dim=c(1,1,dim(qspec)))
  }
  nc <- dim(y.qper)[1]
  if(is.null(sel.f)) sel.f <- c(2:dim(y.qper)[3])
  if(is.null(sel.tau)) sel.tau <- c(2:dim(y.qper)[4])
  if(nc > 1) {
    out <- 0
    for(j in c(1:length(sel.tau))) {
      for(i in c(1:length(sel.f))) {
         S <- matrix(y.qper[,,sel.f[i],sel.tau[j]],ncol=nc)
         S0 <- matrix(qspec[,,sel.f[i],sel.tau[j]],ncol=nc)
         tmp1 <- Mod(sum(diag(S %*% solve(S0))))
         tmp2 <- log(prod(abs(Re(diag(qr(S)$qr)))) / prod(abs(Re(diag(qr(S0)$qr)))))
         out <- out + tmp1 - tmp2 - nc
      }
    }
    out <- out / (length(sel.f) * length(sel.tau))
  } else {
    S <- Re(y.qper[1,1,sel.f,sel.tau])
    S0 <- Re(qspec[,1,sel.f,sel.tau])
    out <- mean(S/S0 - log(S/S0) - 1)    
  }
  out
}


#' @useDynLib qfa, .registration=TRUE
rq.fit.fnb2 <- function(x, y, zeta0, rhs, beta=0.99995,eps=1e-06) {
# modified version of rq.fit.fnb() in 'quantreg' to accept user-provided zeta0 and rhs
# for the Primal-Dual Interior Point algorithm (rqfnb) of Portnoy and Koenker (1997) 
# input:   x = n*p design matrix
#          y = n-vector of response
#      zeta0 = n-vector of initial values of dual variables (x in rqfnb) [ zeta0 = 1-tau in rq.fit.fnb ]
#        rhs = right hand side in the dual formulation [rhs = (1 - tau) * apply(x, 2, sum) in rq.fit.fnb ]
# output: coefficients = p-vector of regression coefficients (primal variables)
#             dualvars = n-vector of dual variables
#                  nit = number of iterations
    n <- length(y)
    p <- ncol(x)
    if (n != nrow(x)) 
        stop("x and y don't match n")
    d <- rep(1, n)
    u <- rep(1, n)
    wn <- rep(0,9*n)
    wp <- rep(0,(p+3)*p)
    wn[1:n] <- zeta0  # initial value of dual variables 
    z <- .Fortran("rqfnb", as.integer(n), as.integer(p), a = as.double(t(as.matrix(x))),
        c = as.double(-y), rhs = as.double(rhs), d = as.double(d), 
        u = as.double(u), beta = as.double(beta), eps = as.double(eps), 
        wn = as.double(wn), wp = as.double(wp), 
        nit = integer(3), info = integer(1))
    if (z$info != 0) 
        warning(paste("Error info = ", z$info, "in stepy: possibly singular design"))
    coefficients <- -z$wp[1:p]
    names(coefficients) <- dimnames(x)[[2]]
    list(coefficients = coefficients, dualvars = z$wn[(2*n+1):(3*n)], nit = z$nit, info = z$info)
}



#' Quantile Autocovariance Function (QACF)
#'
#' This function computes quantile autocovariance function (QACF) from time series 
#' or quantile discrete Fourier transform (QDFT).
#' @param y vector or matrix of time series (if matrix, \code{nrow(y)} = length of time series)
#' @param tau sequence of quantile levels in (0,1)
#' @param y.qdft matrix or array of pre-calculated QDFT (default = \code{NULL}: compute from \code{y} and \code{tau});
#' if \code{y.qdft} is supplied, \code{y} and \code{tau} can be left unspecified
#' @param n.cores number of cores for parallel computing of QDFT if \code{y.qdft = NULL} (default = 1)
#' @param cl pre-existing cluster for repeated parallel computing of QDFT (default = \code{NULL})
#' @return matrix or array of quantile autocovariance function
#' @import 'stats'
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' # compute from time series
#' y.qacf <- qacf(y,tau)
#' # compute from QDFT 
#' y.qdft <- qdft(y,tau) 
#' y.qacf <- qacf(y.qdft=y.qdft)
qacf <- function(y,tau,y.qdft=NULL,n.cores=1,cl=NULL) {
  if(is.null(y.qdft)) y.qdft <- qdft(y=y,tau=tau,n.cores=n.cores,cl=cl)
  return(qdft2qacf(y.qdft))
}


#' Quantile Series (QSER)
#'
#' This function computes quantile series (QSER) from time series or quantile discrete Fourier transform (QDFT).
#' @param y vector or matrix of time series (if matrix, \code{nrow(y)} = length of time series)
#' @param tau sequence of quantile levels in (0,1)
#' @param y.qdft matrix or array of pre-calculated QDFT (default = \code{NULL}: compute from \code{y} and \code{tau});
#' if \code{y.qdft} is supplied, \code{y} and \code{tau} can be left unspecified
#' @param n.cores number of cores for parallel computing of QDFT if \code{y.qdft = NULL} (default = 1)
#' @param cl pre-existing cluster for repeated parallel computing of QDFT (default = \code{NULL})
#' @return matrix or array of quantile series
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' # compute from time series
#' y.qser <- qser(y,tau)  
#' # compute from QDFT
#' y.qdft <- qdft(y,tau)  
#' y.qser <- qser(y.qdft=y.qdft) 
qser <- function(y,tau,y.qdft=NULL,n.cores=1,cl=NULL) {
  if(is.null(y.qdft)) y.qdft <- qdft(y=y,tau=tau,n.cores=n.cores,cl=cl)
  return(qdft2qser(y.qdft))
}


#' Quantile Periodogram (QPER)
#'
#' This function computes quantile periodogram (QPER) from time series 
#' or quantile discrete Fourier transform (QDFT).
#' @param y vector or matrix of time series (if matrix, \code{nrow(y)} = length of time series)
#' @param tau sequence of quantile levels in (0,1)
#' @param y.qdft matrix or array of pre-calculated QDFT (default = \code{NULL}: compute from \code{y} and \code{tau});
#' if \code{y.qdft} is supplied, \code{y} and \code{tau} can be left unspecified
#' @param n.cores number of cores for parallel computing of QDFT if \code{y.qdft = NULL} (default = 1)
#' @param cl pre-existing cluster for repeated parallel computing of QDFT (default = \code{NULL})
#' @return matrix or array of quantile periodogram
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' # compute from time series
#' y.qper <- qper(y,tau)  
#' # compute from QDFT 
#' y.qdft <- qdft(y,tau) 
#' y.qper <- qper(y.qdft=y.qdft)  
qper <- function(y,tau,y.qdft=NULL,n.cores=1,cl=NULL) {
  if(is.null(y.qdft)) y.qdft <- qdft(y=y,tau=tau,n.cores=n.cores,cl=cl)
  return(qdft2qper(y.qdft))
}



# Quantile Spectrum from AR Model of Quantile Series
#
# This function computes quantile spectrum (QSPEC) from an AR model of Quantile Series (QSER).
# @param fit object of AR model from \code{qser2sar()} or \code{qser2ar()}
# @param freq sequence of frequencies in [0,1) (default = \code{NULL}: all Fourier frequencies)
# @return a list with the following elements:
#   \item{spec}{matrix or array of quantile spectrum}
#   \item{freq}{sequence of frequencies}
# @export
ar2qspec<-function(fit,freq=NULL) {
  V <- fit$V
  A <- fit$A
  n <- fit$n
  p <- fit$p
  if(is.vector(V)) {
    V <- array(V,dim=c(1,1,length(V)))
	if(p > 0) A <- array(A,dim=c(1,1,dim(A)))
  }
  nc <- dim(V)[1]
  ntau <- dim(V)[3]
  if(is.null(freq)) freq<-c(0:(n-1))/n
  nf <- length(freq)
  spec <- array(0,dim=c(nc,nc,nf,ntau))
  for(ell in c(1:ntau)) {
    for(jj in c(1:nf)) {
	  if(p > 0) {
        U <- matrix(0,nc,nc)
        for(j in c(1:p)) {
          U <- U + A[,,j,ell] * exp(-1i*2*pi*j*freq[jj])
        }
        U <- matrix(solve(diag(nc) - U),ncol=nc)
        spec[,,jj,ell] = U %*% (matrix(V[,,ell],nrow=nc,ncol=nc) %*% t(Conj(U)))
      } else {
	    spec[,,jj,ell] = matrix(V[,,ell],nrow=nc,ncol=nc)
	  }
    }
  }
  if(nc==1) spec <- spec[1,1,,]
  return(list(spec=spec, freq=freq))
}


#' Spline Autoregression (SAR) Model of Quantile Series
#'
#' This function fits spline autoregression (SAR) model to quantile series (QSER).
#' @param y.qser matrix or array of pre-calculated QSER at \code{tau0} using, e.g., \code{qser()}
#' @param tau0 quantile levels used to compute \code{y.qser}
#' @param tau  quantile levels for evaluation (\code{min(tau0)} <= \code{tau} <= \code{max(tau0)}; 
#' default = \code{tau0})
#' @param p order of SAR model (default = \code{NULL}: automatically selected by AIC)
#' @param order.max maximum order for AIC if \code{p = NULL} (default = \code{NULL}: determined by \code{stats::ar()})
#' @param spar penalty parameter alla \code{stats::smooth.spline()} (default = \code{NULL}: automatically selected)
#' @param method criterion for penalty parameter selection: \code{"GCV"} (default), \code{"AIC"}, or \code{"BIC"} 
#' @param type type of quantile smoothing for residual variance: 1 = direct (default) or 2 = square root
#' @param weights sequence of weights (default = \code{rep(1,length(tau0))})
#' @param interval interval for \code{spar} optimization (default = \code{c(-1.5,1.5)})
#' @return a list with the following elements:
#'   \item{A}{matrix or array of SAR coefficients}
#'   \item{V}{vector or matrix of SAR residual covariance}
#'   \item{p}{order of SAR model}
#'   \item{spar}{penalty parameter}
#'   \item{n}{length of time series}
#'   \item{tau}{quantile levels for evaluation}
#'   \item{tau0}{quantile levels for fitting}
#;   \item{type}{type of quantile smoothing for residual variance}
#'   \item{weights}{weights in penalty function}
#'   \item{sdm}{spline design matrices}
#'   \item{fit}{object containing details of SAR fit}
#' @import 'stats'
#' @import 'splines'
#' @export
qser2sar <- function(y.qser,tau0,tau=tau0,p=NULL,order.max=NULL,spar=NULL,method=c("GCV","AIC","BIC"),
            type=c(1,2),weights=rep(1,length(tau0)),interval=c(-1.5,1.5)) {

  create_spline_matrix <- function(tau0,tau) {
  # create spline matrix and its second derivative for SQR
  # input:    tau0 = quantiles for fitting
  #            tau = quantiles for evaluation
  # output: bsMat  = matrix of basis functions
  #         dbsMat = matrix of second derivatives
  #         bsMat2 = maxtrix of basis function for interpolation to tau
  # imports: 'splines'
    knots <- stats::smooth.spline(tau0,tau0)$fit$knot 
    # rescale tau0 and tau to [0,1] for spline matrix to be consistent with smooth.spline
    bsMat <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)))
    dbsMat <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)),derivs=2)
    bsMat2 <- splines::splineDesign(knots,(tau-min(tau))/(max(tau)-min(tau)))
    return(list(bsMat=bsMat,dbsMat=dbsMat,bsMat2=bsMat2,knots=knots))
  }

  sar_rescale_penalty <- function(Z2,D2) {
  # rescalng factor for the SAR penalty parameter alla smooth.spline
    r <- sum(diag(Z2))
    r <- r / sum(diag(D2))
	r
  }

  sar_obj <- function(spar,Z2,D2,YZ,Y,Z,r,nc,L,n,p,method=c("GCV","AIC","BIC"),return.solution=FALSE) {
  # this function computes the criterion for penalty parameter selection in SAR using spar
  # with the option of returning the details of fit
    lam <- r * 256^(3*spar-1)
    # compute least-squares solution
    Z2 <- Z2 + lam * D2
    Omega <- solve(Z2)
    Theta <- YZ %*% Omega
    # compute residual covariance matrix
    df <- 0
    rss <- 0
    V0 <- array(0,dim=c(nc,nc,L))
    for(l in c(1:L)) {
      df <- df + sum(diag(crossprod(Z[[l]],Omega) %*% Z[[l]])) * nc
      resid <- Y[[l]] - Theta %*% Z[[l]]
      if(method[1]=="GCV") {
        rss <- rss + sum(resid^2)
      }
	  if(method[1] %in% c("AIC","BIC") | return.solution) {
        V0[,,l] <- matrix(tcrossprod(resid),ncol=nc) / n
      }
    }
    # compute parameter selection criterion
    crit <- 0
    if(method[1]=="AIC") {
      for(l in c(1:L)) {
        if(nc==1) {
          crit <- crit + c(abs(V0[,,l])) 
        } else {
          crit <- crit + c(determinant(V0[,,l])$modulus)
        }
      }
      crit <- n*crit + 2*df
    }
    if(method[1]=="BIC") {
      for(l in c(1:L)) {
        if(nc==1) {
          crit <- crit + c(abs(V0[,,l]))
        } else {
          crit <- crit + c(determinant(V0[,,l])$modulus)
        }
      }
      crit <- n*crit + log(n)*df
    }
	if(method[1]=="GCV") {
      crit <- (rss/(L*(n-p))) / (1-df/(L*(n-p)))^2
    }
    if(return.solution) {
      return(list(crit=crit,Theta=Theta,V0=V0,Omega=Omega,lam=lam,df=df,method=method[1]))
    } else { 
      return(crit)
    }
  }
  
  sar_var_smooth <- function(fit,tau0,tau,type=c(1,2)) {
  # smooth and interpolate residual covariance matrix
    V0 <- fit$V0
	p <- fit$p
    nc <- dim(V0)[1]
	nq <- dim(V0)[3]
    ntau <- length(tau)
    V <- array(0,dim=c(nc,nc,ntau))
    U <- V0
    if(type[1]==2) {
      if(nc==1) {
		U[1,1,] <- sqrt(U[1,1,]) 
	  } else {
	    for(ell in c(1:nq)) {
	      tmp <- svd(U[,,ell])
		  U[,,ell] <- tmp$u %*% diag(sqrt(tmp$d))
        }
	  }
	}
    if(p==0) {
      lam <- NULL
      for(k in c(1:nc)) {
	    for(kk in c(1:nc)) {   
	      lam <- c(lam,stats::smooth.spline(tau0,U[k,kk,])$lambda)
	    }
	  }
      fit$lam <- mean(lam)
	}
    if(type[1]==2) {
      for(k in c(1:nc)) {
  	    for(kk in c(1:nc)) {
  	      fit.ss <- stats::smooth.spline(tau0,U[k,kk,],lambda=fit$lam)
	      V[k,kk,] <- stats::predict(fit.ss,tau)$y
	    }
	  } 
	  for(ell in c(1:ntau)) {
		V[,,ell] <- V[,,ell] %*% t(V[,,ell])
      }
	} else {
      for(k in c(1:nc)) {
        fit.ss <- stats::smooth.spline(tau0,U[k,k,],lambda=fit$lam)
        V[k,k,] <- stats::predict(fit.ss,tau)$y
        if(k < nc) {
          for(kk in c((k+1):nc)) {
            fit.ss <- stats::smooth.spline(tau0,U[k,kk,],lambda=fit$lam)
            V[k,kk,] <- stats::predict(fit.ss,tau)$y
	      }
          V[kk,k,] <- V[k,kk,]
        }
	  }
    }
	# fix possible nonpositive-definite matrix
    eps <- 1e-8
	lam0 <- rep(0,ntau)
    for(ell in c(1:ntau)) {
	  lam0[ell] <- min(eigen(V[,,ell])$values)
	}
	lam0 <- min(lam0)
    if(lam0 < 0.0) {
	  for(ell in c(1:ntau)) {
        V[,,ell] <- V[,,ell] + diag(eps+abs(lam0),nrow=nc,ncol=nc)
      }
    }
    return(V)
  }
  

  if(is.matrix(y.qser)) y.qser <- array(y.qser,dim=c(1,dim(y.qser)))
  nc <- dim(y.qser)[1]
  n <- dim(y.qser)[2]
  L <- dim(y.qser)[3]
  ntau <- length(tau)
  type <- type[1]
  
  if( length(tau0) != L ) {
    stop("mismatched dimension: y.qser should be computed at tau0!")
  }
  if( length(tau0) < 10 ) {
    stop("too few quantile levels for fitting (must be 10 or more)!")
  }
  if( min(tau0) > min(tau) | max(tau0) < max(tau) ) {
      stop("tau is not in the range of tau0!")
  }
  if( is.null(p) ) {
    if( is.null(order.max) ) {
       stop("order.max unspecified!")	
	} else {
      if( order.max >= n/2 ) {
        stop("order.max too large (must be less than n/2)!")
      }
    }
  }
  if( !is.null(p) ) {
    if(p >= n) {
      stop("order p too large (must be less than n)!")
    }
  }
  if( length(weights) != length(tau0) ) {  
      stop("length of weights not equal length of tau0!")
  }
  if( is.null(spar) & !(method[1] %in% c("GCV","AIC","BIC")) ) {
    stop("method of smoothing parameter selection not GCV, AIC, or BIC!")
  }

  # create L-by-K spline design matrix with knots provided by smooth.spline
  sdm <- create_spline_matrix(tau0,tau)
  K <- ncol(sdm$bsMat)
  # create weights
  w <- weights
  sel.tau0 <- c(1:L)
  
  # fit VAR to obtain AIC for order selection
  if(is.null(p)) {
    aic <- NULL
    for(l in c(1:L)) {
      mu <- apply(matrix(y.qser[,,sel.tau0[l]],nrow=nc),1,mean)
      aic <- cbind(aic,stats::ar(t(matrix(y.qser[,,sel.tau0[l]]-mu,nrow=nc)),
                   order.max = order.max, method=c("yule-walker"))$aic)
    }
    aic <- apply(aic,1,mean)
    # optimal order minimizes average aic
    p <- as.numeric(which.min(aic)) - 1
  }

  # SAR model of order 0
  if(p == 0) {
	 # compute residual variance
     V0 <- array(0,dim=c(nc,nc,L)) 
     for(l in c(1:L)) {
       mu <- apply(matrix(y.qser[,,sel.tau0[l]],nrow=nc),1,mean)
       V0[,,l] <- matrix(tcrossprod(matrix(y.qser[,,sel.tau0[l]]-mu,nrow=nc)),ncol=nc) / n
     }
	 fit <- NULL
	 fit$V0 <- V0
	 fit$p <- 0
	 # smooth residual variance
	 V <- sar_var_smooth(fit,tau0,tau,type)
     if(nc==1) {
       V <- V[1,1,]
	   fit$V0 <- fit$V0[1,1,]
     }
     return(list(A=NULL,V=V,p=p,spar=NULL,tau=tau,n=n,tau0=tau0,type=type,weights=weights,fit=fit))
  }
  
  # SAR model of order p > 0 
  Z2 <- matrix(0,nrow=K*nc*p,ncol=K*nc*p)
  D2 <- matrix(0,nrow=K*nc*p,ncol=K*nc*p)
  YZ <- matrix(0,nrow=nc,ncol=K*nc*p)
  Z <- list()
  Y <- list()
  for(l in c(1:L)) {
    Zt <- matrix(0,nrow=n-p,ncol=K*nc*p)
    mu <- apply(matrix(y.qser[,,sel.tau0[l]],nrow=1),1,mean)
    for(j in c(1:p)) {
      Yt <- t(matrix(y.qser[,c((p-j+1):(n-j)),sel.tau0[l]]-mu,nrow=nc))
      j0 <- (j-1)*K*nc
      for(k in c(1:K)) {
        Zt[,(j0+(k-1)*nc+1):(j0+k*nc)] <- sdm$bsMat[l,k] * Yt
      }
    }
    D0 <- matrix(0,nrow=nc,ncol=K*nc)
    for(k in c(1:K)) {
      D0[,((k-1)*nc+1):(k*nc)] <- diag(sdm$dbsMat[l,k],nrow=nc,ncol=nc)
    }
	Dt <- matrix(0,nrow=nc*p,ncol=K*nc*p)	
    for(j in c(1:p)) {
      Dt[((j-1)*nc+1):(j*nc),((j-1)*K*nc+1):(j*K*nc)] <- D0
    }
    Yt <- t(matrix(y.qser[,c((p+1):n),sel.tau0[l]]-mu,nrow=nc))
    Z2 <- Z2 + crossprod(Zt)
    YZ <- YZ + crossprod(Yt,Zt)
    D2 <- D2 + w[l] * crossprod(Dt)
    Z[[l]] <- t(Zt)
    Y[[l]] <- t(Yt)
  }
  # clean up
  rm(Zt,Dt,Yt,D0)
  # compute scaling factor of penalty parameter alla smooth.spline
  r <- sar_rescale_penalty(Z2,D2)
  # select optimal penalty parameter
  if(is.null(spar)) {
    spar <- stats::optimize(sar_obj,interval=interval,
	        Z2=Z2,D2=D2,YZ=YZ,Y=Y,Z=Z,r=r,nc=nc,L=L,n=n,p=p,method=method[1])$minimum
  }
  # compute SAR solution
  fit <- sar_obj(spar,Z2,D2,YZ,Y,Z,r,nc,L,n,p,method=method[1],return.solution=TRUE)
  # clean up    
  rm(Z,Y,YZ,Z2,D2)
  # compute SAR coefficient functions
  A <- array(0,dim=c(nc,nc,p,ntau))
  for(j in c(1:p)) {
    theta <- matrix(fit$Theta[,((j-1)*K*nc+1):(j*K*nc)],nrow=nc)
    for(ell in c(1:ntau)) {
	  Phi <- matrix(0,nrow=K*nc,ncol=nc)
      for(k in c(1:K)) {
        Phi[((k-1)*nc+1):(k*nc),] <- diag(sdm$bsMat2[ell,k],nrow=nc,ncol=nc)
      }
      A[,,j,ell] <- theta %*% Phi
    }
  }
  fit$p <- p
  # smooth residual covariance matrix
  V <- sar_var_smooth(fit,tau0,tau,type)
  
  # special case of nc=1
  if(nc==1) {
    A <- matrix(A[1,1,,],nrow=p)
    V <- V[1,1,]
    fit$V0 <- fit$V0[1,1,]
  }
  return(list(A=A,V=V,p=p,spar=spar,tau=tau,n=n,tau0=tau0,type=type,weights=weights,sdm=sdm,fit=fit))
}


#' Spline Autoregression (SAR) Estimator of Quantile Spectrum
#'
#' This function computes spline autoregression (SAR) estimate of quantile spectrum.
#' @param y vector or matrix of time series (if matrix, \code{nrow(y)} = length of time series)
#' @param y.qser matrix or array of pre-calculated QSER at \code{tau0} (default = \code{NULL}: compute from \code{y} and \code{tau0});
#' if supplied, \code{y} can be left unspecified
#' @param tau0 quantile levels used to computer QSER (must be supplied with or without \code{y.qser})
#' @param tau  quantile levels for evaluation (\code{min(tau0)} <= \code{tau} <= \code{max(tau0)}; 
#' default = \code{tau0})
#' @param p order of SAR model (default = \code{NULL}: automatically selected by AIC)
#' @param order.max maximum order for AIC if \code{p = NULL} (default = \code{NULL}: determined by \code{stats::ar()})
#' @param spar penalty parameter alla \code{stats::smooth.spline()} (default = \code{NULL}: automatically selected)
#' @param method criterion for penalty parameter selection: \code{"GCV"} (default), \code{"AIC"}, or \code{"BIC"}
#' @param type type of quantile smoothing for residual variance: 1 = direct (default) or 2 = square root
#' @param weights sequence of weights (\code{length(weights) = length(tau0)}; default: \code{weights = rep(1,length(tau0))})
#' @param interval interval for spar optimization (default: \code{interval = c(-1.5,1.5)})
#' @param freq sequence of frequencies in [0,1) (default = \code{NULL}: all Fourier frequencies)
#' @param n.cores number of cores for parallel computing of QDFT if \code{y.qser = NULL} (default = 1)
#' @param cl pre-existing cluster for repeated parallel computing of QDFT (default = \code{NULL})
#' @return a list with the following elements:
#'   \item{spec}{matrix or array of SAR quantile spectrum}
#'   \item{freq}{sequence of frequencies}
#'   \item{fit}{object of SAR model}
#'   \item{qser}{matrix or array of quantile series if \code{y.qser = NULL}}
#' @import 'stats'
#' @import 'splines'
#' @export
#' @examples
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' n <- length(y1)
#' ff <- c(0:(n-1))/n
#' sel.f <- which(ff > 0 & ff < 0.5)
#' # compute from time series
#' y.sar <- qspec.sar(cbind(y1,y2),tau0=tau,p=1)
#' qfa.plot(ff[sel.f],tau,Re(y.sar$spec[1,1,sel.f,]))
#' # compute from quantile series
#' y.qser <- qser(cbind(y1,y2),tau)
#' y.sar <- qspec.sar(y.qser=y.qser,tau0=tau,p=1)
#' qfa.plot(ff[sel.f],tau,Re(y.sar$spec[1,1,sel.f,]))
qspec.sar <- function(y,y.qser=NULL,tau0,tau=tau0,p=NULL,order.max=NULL,spar=NULL,method=c("GCV","AIC","BIC"),
type=c(1,2),weights=rep(1,length(tau0)),interval=c(-1.5,1.5),freq=NULL,n.cores=1,cl=NULL) {
  return.qser <- FALSE
  if(is.null(y.qser)) { 
    y.qser <- qser(y,tau0,n.cores=n.cores,cl=cl)
    return.qser <- TRUE
  }
  fit <- qser2sar(y.qser,tau0=tau0,tau=tau,p=p,order.max=order.max,spar=spar,
                  method=method[1],type=type[1],weights=weights,interval=interval)
  tmp <- ar2qspec(fit)
  if(return.qser) {
    return(list(spec=tmp$spec,freq=tmp$freq,fit=fit,qser=y.qser)) 
  } else {
    return(list(spec=tmp$spec,freq=tmp$freq,fit=fit))
  }
}



#' Quantile Periodogram Type II (QPER2)
#'
#' This function computes type-II quantile periodogram for univariate time series.
#' @param y univariate time series
#' @param freq sequence of frequencies in [0,1)
#' @param tau sequence of quantile levels in (0,1)
#' @param weights sequence of weights in quantile regression (default = \code{NULL}: weights equal to 1)
#' @param n.cores number of cores for parallel computing (default = 1)
#' @param cl pre-existing cluster for repeated parallel computing (default = \code{NULL})
#' @return matrix of quantile periodogram evaluated on \code{freq * tau} grid
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @import 'quantreg'
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' n <- length(y)
#' ff <- c(0:(n-1))/n
#' sel.f <- which(ff > 0 & ff < 0.5)
#' y.qper2 <- qper2(y,ff,tau)
#' qfa.plot(ff[sel.f],tau,Re(y.qper2[sel.f,]))
qper2 <- function(y,freq,tau,weights=NULL,n.cores=1,cl=NULL) {

 rsoid2 <- function(tt,ff,a,b) {
 # this function computes trigonometric regression function
  one <- rep(1,length(tt))
  tmp <- (one %o% a) * cos(2*pi*(tt %o% ff)) + (one %o% b) * sin(2*pi*(tt %o% ff))
  tmp <- apply(tmp,1,sum)
  return(tmp)
 }

 qr.cost <- function(y,tau,weights=NULL) {
 # this function computes the cost of quantile regression
   n <- length(y)
   if(is.null(weights)) weights<-rep(1,n)
   tmp <- tau*y
   sel <- which(y < 0)
   if(length(sel) > 0) tmp[sel] <- (tau-1)*y[sel]
   return(sum(tmp*weights,na.rm=T))
 }
 
 tqr.cost <- function(yy,ff,tau,weights=NULL,method='fn') {
 # this function computes the cost of trigonometric quantile regression
 # for a vector of quantile levels tau at a single frequency ff
    n <- length(yy)
    tt <- c(1:n)
	ntau <- length(tau)
    if(ff == 0.5) {
      fit <- suppressWarnings(quantreg::rq(yy ~ cos(2*pi*ff*tt),method=method,tau=tau,weights=weights))
      fit$coefficients <- rbind(fit$coefficients,rep(0,ntau))
    }
    if(ff == 0) {
	  fit <- suppressWarnings(quantreg::rq(yy ~ 1,method=method,tau=tau,weights=weights))
      fit$coefficients <- rbind(fit$coefficients,rep(0,ntau),rep(0,ntau))
    }
    if(ff != 0.5 & ff != 0) {
      fit <- suppressWarnings(quantreg::rq(yy ~ cos(2*pi*ff*tt)+sin(2*pi*ff*tt),method=method,tau=tau,weights=weights))
    }
    fit$coefficients <- matrix(fit$coefficients,ncol=ntau)
    if(ntau == 1) fit$coefficients <- matrix(fit$coefficients,ncol=1)
    cost <- rep(NA,ntau)
    for(i.tau in c(1:ntau)) {
       resid <- yy - fit$coefficients[1,i.tau]-rsoid2(tt,ff,fit$coefficients[2,i.tau],fit$coefficients[3,i.tau])
       cost[i.tau] <- qr.cost(resid,tau[i.tau],weights)
    }
    return(cost)
 }
 
 n <- length(y)
 nf <- length(freq)
 ntau <- length(tau)
 if(is.null(weights)) weights<-rep(1,n)
 keep.cl <- TRUE
 if(n.cores>1 & is.null(cl)) {
    cl <- parallel::makeCluster(n.cores)
    parallel::clusterExport(cl, c("rq"))
    doParallel::registerDoParallel(cl)
    keep.cl <- FALSE
 }
 `%dopar%` <- foreach::`%dopar%`
 `%do%` <- foreach::`%do%`
 # cost without trigonometric regressor
 cost0 <- tqr.cost(y,0,tau,weights)
 # cost with trigonometric regressor
 i <- 0  
 if(n.cores>1) {
   out <- foreach::foreach(i=1:nf) %dopar% { tqr.cost(y,freq[i],tau,weights) }
 } else {
   out <- foreach::foreach(i=1:nf) %do% { tqr.cost(y,freq[i],tau,weights) }
 }
 out <- matrix(unlist(out),ncol=ntau,byrow=T)
 out <- matrix(rep(cost0,nf),ncol=ntau,byrow=T) - out
 if(n.cores>1 & !keep.cl) {
    parallel::stopCluster(cl)
    cl <-NULL
 }
 out[out<0]<-0
 if(ntau==1) out<-c(out)
 return(out)
}

  
#' Extraction of SAR Coefficients for Granger-Causality Analysis
#'
#' This function extracts the spline autoregression (SAR) coefficients from an SAR model for Granger-causality analysis.
#' See \code{sar.gc.bootstrap} for more details regarding the use of \code{index}.
#' @param fit object of SAR model from \code{qser2sar()} or \code{qspec.sar()$fit}
#' @param index a pair of component indices for multiple time series 
#' or a sequence of lags for single time series (default = \code{c(1,2)})
#' @return matrix of selected SAR coefficients (number of lags by number of quantiles)
#' @export
#' @examples
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y.sar <- qspec.sar(cbind(y1,y2),tau0=tau,p=1)
#' A <- sar.gc.coef(y.sar$fit,index=c(1,2))
sar.gc.coef <- function(fit,index=c(1,2)) {

  test_for_validity <- function(index,nc,p) {
    if(p==0) {
      stop("method not applicable to order 0 models!")	  
    }
    if(nc==1) {
      if(sum(index < 1 | index > p) > 0) {
        stop(paste0("index outside the range of lags 1-",p,"!"))   
      }
    } else {
      if(length(index) != 2) {
        stop(paste0("length of index not equal 2!"))	  	
      }
      if(sum(index < 1 | index > nc) > 0) {
        stop(paste0("index outside the range of time series components 1-",nc,"!"))	   
      }
    }
  }

  if(is.vector(fit$V)) {
	nc <- 1
  } else {
    nc <- dim(fit$V)[1]	  
  }
  p <- fit$p
  test_for_validity(index,nc,p)
  if(nc == 1) {
    return(matrix(fit$A[index,],nrow=length(index)))
  } else {
    return(matrix(fit$A[index[1],index[2],,],nrow=p))
  }
}


# This function computes the residual of SAR fit
sar.residual <- function(y.qser,fit) {
    p <- fit$p
    n <- fit$n
    ntau <- dim(y.qser)[3]
    if(is.matrix(y.qser)) y.qser <- array(y.qser,dim=c(1,dim(y.qser)))  
    nc <- dim(y.qser)[1] 
    mu <- apply(y.qser,c(1,3),mean)
	y.qser0 <- y.qser
	if(nc==1) mu <- matrix(mu,nrow=1)
    for(l in c(1:ntau)) y.qser0[,,l] <- y.qser0[,,l] - mu[,l]
	if(p > 0) {
      tmp <- array(0,dim=c(nc,n,ntau))
      for(l in c(1:ntau)) {
	    if(nc > 1) {
          for(tt in c((p+1):n)) {
            for(j in c(1:p)) tmp[,tt,l] <- tmp[,tt,l] + fit$A[,,j,l] %*% y.qser0[,tt-j,l]
          }
		} else {
          for(tt in c((p+1):n)) {
            for(j in c(1:p)) tmp[1,tt,l] <- tmp[1,tt,l] + fit$A[j,l] %*% y.qser0[1,tt-j,l]
          }		
		}
      }
	  resid <- y.qser0 - tmp
      resid <- resid[,-c(1:p),]
	} else {
	  resid <- y.qser0
	}
    resid
}


# This function check the validity of parameters for SAR-based causality analysis
sar.gc.param.validity <- function(index,nc,p) {
    if(p==0) {
      stop("method not applicable to order 0 models!")	  
    }
    if(nc==1) {
      if(sum(index < 1 | index > p) > 0) {
        stop(paste0("index outside the range of lags 1-",p,"!"))   
      }
    } else {
      if(length(index) != 2) {
        stop(paste0("length of index not equal 2!"))	  	
      }
      if(sum(index < 1 | index > nc) > 0) {
        stop(paste0("index outside the range of time series components 1-",nc,"!"))	   
      }
    }
}



#' Bootstrap Simulation of SAR Coefficients for Granger-Causality Analysis
#'
#' This function simulates bootstrap samples of selected spline autoregression (SAR) coefficients 
#' for Granger-causality analysis based on the SAR model of quantile series (QSER) under H0: 
#' (a) for multiple time series, the second series specified in \code{index} is not causal 
#' for the first series specified in \code{index};
#' (b) for single time series, the series is not causal at the lags specified in \code{index}.
#' @param y.qser matrix or array of QSER from \code{qser()} or \code{qspec.sar()$qser}
#' @param fit object of SAR model from \code{qser2sar()} or \code{qspec.sar()$fit}
#' @param index a pair of component indices for multiple time series 
#' or a sequence of lags for single time series (default = \code{c(1,2)})
#' @param nsim number of bootstrap samples (default = 1000)
#' @param method method of residual calculation: \code{"ar"} (default) or \code{"sar"}
#' @param refit if \code{TRUE} the SAR model is refit under H0 (default = \code{FALSE})
#' @param n.cores number of cores for parallel computing (default = 1)
#' @param mthreads if \code{FALSE} (default), set \code{RhpcBLASctl::blas_set_num_threads(1)}
#' @param seed seed for random sampling (default = \code{1234567})
#' @return array of simulated bootstrap samples of selected SAR coefficients
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @import 'RhpcBLASctl'
#' @export
#' @examples
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y.sar <- qspec.sar(cbind(y1,y2),tau0=tau,p=1)
#' A.sim <- sar.gc.bootstrap(y.sar$qser,y.sar$fit,index=c(1,2),nsim=5)
sar.gc.bootstrap <- function(y.qser,fit,index=c(1,2),nsim=1000,method=c("ar","sar"),
                    refit=FALSE,n.cores=1,mthreads=FALSE,seed=1234567) {

  sarc_refit <- function(fit.sar,index=c(1,2)) {
  # compute constrainted SAR coefficients from unconstraint fit.sar
  # where a constraint is imposed on the (index[1],index[2]) entry of the SAR coefficients
    Omega <- fit.sar$fit$Omega
    Theta <- fit.sar$fit$Theta
    lam <- fit.sar$fit$lam
    p <- fit.sar$fit$p
    sdm <- fit.sar$sdm
    nc <- dim(Theta)[1]
    K <- dim(sdm$bsMat)[2]
    ntau <- dim(sdm$bsMat2)[1]
    e1 <- rep(0,nc)
    e1[index[1]] <- 1
    e2 <- rep(0,nc)
    e2[index[2]] <- 1
    Q <- kronecker(diag(K*p),t(e2))
    R <- Q %*% Omega
    # compute constraint vector
    theta <- rep(0,K*p)
    # compute Lagrange multipler
    xi <- solve(R %*% t(Q))
    xi <- xi %*% (Q %*% t(Theta) %*% e1 - theta)
    # compute SARC 
    Theta <- Theta - kronecker(e1,t(xi) %*% R)
    # compute SAR coefficient functions
    A <- fit.sar$A
    for(j in c(1:p)) {
      theta <- matrix(Theta[,((j-1)*K*nc+1):(j*K*nc)],nrow=nc)
      for(ell in c(1:ntau)) {
	    Phi <- matrix(0,nrow=K*nc,ncol=nc)
        for(k in c(1:K)) {
          Phi[((k-1)*nc+1):(k*nc),] <- diag(sdm$bsMat2[ell,k],nrow=nc,ncol=nc)
        }
        A[,,j,ell] <- theta %*% Phi
      }
    }
    return(A)
  }

  simulate_qser_array <- function(resid,sample.idx,fit,idx0,refit=FALSE) {
  # this function simulates quantile series under H0 from bootstrap samples of residuals
  # input: resid = nc*(n-p)*ntau-array of residuals from qser2ar()
  #        sample.idx = nn-vector of sampled time indices
  #        fit = SAR model from qser2sar()  
  #        idx0 = indices of SAR coefficients to set to zero under H0
  #        refit = if TRUE, SAR model is refit under the constraint H0
  # output: nc*n*ntau array of simulated quantile series
      p <- fit$p
      nc <- dim(resid)[1]
	  n <- dim(resid)[2] + p
      ntau <- dim(resid)[3]
      nn <- length(sample.idx)
	  A <- fit$A
	  if(refit) A <- sarc_refit(fit,idx0)
      resid.sim <- array(resid[,sample.idx,],dim=c(nc,nn,ntau))
      qser.sim <- array(0,dim=c(nc,nn,ntau))
      for(ell in c(1:ntau)) {
        mu <- apply(resid.sim[,,ell],1,mean) 
        for(jj in c((p+1):nn)) {
          if(nc == 1) {
            for(j in c(1:p)) {
              A0 <- A[j,ell]
              if(j %in% idx0) A0 <- 0
              A0 <- matrix(A0,ncol=nc,nrow=nc)	  
              qser.sim[,jj,ell] <- qser.sim[,jj,ell] + A0 %*% matrix(qser.sim[,jj-j,ell],nrow=nc)
            }		  
          } else {
            for(j in c(1:p)) {
              A0 <- A[,,j,ell]
              A0[idx0[1],idx0[2]] <- 0
              qser.sim[,jj,ell] <- qser.sim[,jj,ell] + A0 %*% matrix(qser.sim[,jj-j,ell],nrow=nc)
            } 
          }
          qser.sim[,jj,ell] <- qser.sim[,jj,ell] + (resid.sim[,jj,ell]-mu)
        }
      }
      qser.sim <- qser.sim[,c((nn-n+1):nn),]
      return(qser.sim)
  }

  if(is.matrix(y.qser)) y.qser <- array(y.qser,dim=c(1,dim(y.qser)))  
  nc <- dim(y.qser)[1]
  n <- dim(y.qser)[2]
  ntau <- dim(y.qser)[3]
  p <- fit$p
  
  sar.gc.param.validity(index,nc,p) 
  
  if(method[1]=="ar") {
    # compute residuals from AR model fitted to QSER for each quantile
    resid <- qser2ar(y.qser,p=p)$resid
  } else {
    # compute residuals from SAR model fitted to QSER
    resid <- sar.residual(y.qser,fit)
  }
  if(nc==1) resid <- array(resid,dim=c(1,dim(resid)))
   
  # generate bootstrap sampling of time stamps
  nn <- 2*n
  set.seed(seed)
  idx <- list()
  for(sim.idx in c(1:nsim)) {
    idx[[sim.idx]] <- sample(c(1:dim(resid)[2]),nn,replace=TRUE)
  }
  
  `%dopar%` <- foreach::`%dopar%`
  `%do%` <- foreach::`%do%`
  if(n.cores > 1) {
  	if(!mthreads) {
	  ##sink(tempfile())
      RhpcBLASctl::blas_set_num_threads(1)
	  ##sink()
    }
    cl <- parallel::makeCluster(n.cores)
    parallel::clusterExport(cl, c("qser2sar","sar.gc.coef"))
    doParallel::registerDoParallel(cl)
    out <- foreach::foreach(sim.idx=c(1:nsim)) %dopar% {
	  qser.sim <- simulate_qser_array(resid,idx[[sim.idx]],fit,index,refit)
      fit.sim <- qser2sar(qser.sim,tau0=fit$tau0,tau=fit$tau,p=fit$p,spar=fit$spar,type=fit$type,weights=fit$weights)
      A.sim <- sar.gc.coef(fit.sim,index)
    }
    parallel::stopCluster(cl)
	if(!mthreads) {
      ##sink(tempfile())
      RhpcBLASctl::blas_set_num_threads(Sys.getenv('OPENBLAS_NUM_THREADS'))
      ##sink()
    }
  } else {
    out <- foreach::foreach(sim.idx=c(1:nsim)) %do% {
	  qser.sim <- simulate_qser_array(resid,idx[[sim.idx]],fit,index,refit)
      fit.sim <- qser2sar(qser.sim,tau0=fit$tau0,tau=fit$tau,p=fit$p,spar=fit$spar,type=fit$type,weights=fit$weights)
      A.sim <- sar.gc.coef(fit.sim,index)
    }
  }
  A.sim <- array(0,dim=c(nsim,dim(out[[1]])))
  for(sim.idx in c(1:nsim)) {
    A.sim[sim.idx,,] <- out[[sim.idx]]
  }
  return(A.sim)
}


#' Wald Test and Confidence Band for Granger-Causality Analysis
#'
#' This function computes Wald test and confidence band for Granger-causality 
#' using bootstrap samples generated by \code{sar.gc.bootstrap()} 
#' based the spline autoregression (SAR) model of quantile series (QSER).
#' @param A matrix of selected SAR coefficients 
#' @param A.sim simulated bootstrap samples from \code{sar.gc.bootstrap()}
#' @param sel.lag indices of time lags for Wald test (default = \code{NULL}: all lags)
#' @param sel.tau indices of quantile levels for Wald test (default = \code{NULL}: all quantiles)
#' @return a list with the following elements:
#'   \item{test}{list of Wald test result containing \code{wald} and \code{p.value}}
#'   \item{A.u}{matrix of upper limits of 95\% confidence band of \code{A}}
#'   \item{A.l}{matrix of lower limits of 95\% confidence band of \code{A}}
#' @import 'MASS'
#' @export
#' @examples
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y.sar <- qspec.sar(cbind(y1,y2),tau0=tau,p=1)
#' A <- sar.gc.coef(y.sar$fit,index=c(1,2))
#' A.sim <- sar.gc.bootstrap(y.sar$qser,y.sar$fit,index=c(1,2),nsim=5)
#' y.gc <- sar.gc.test(A,A.sim)
sar.gc.test <- function(A,A.sim,sel.lag=NULL,sel.tau=NULL) {
   
  wald_test <- function(A,A.sim,sel.lag,sel.tau) {
  # this function computes the bootstrap Wald statistic and p-value
  # of selected elements in A using bootstrap samples in A.sim for Granger Causality test
  # input:   A     = p*ntau-matrix of selected SAR(p) coefficients for all lags and quantiles
  #          A.sim = nsim*p*ntau-array of bootstrap samples associated with A
  #          sel.lap = selected lags for testing (NULL for all lags )
  #          sel.tau = selected quantile levels for testing (NULL for all quantiles)
  # imports: MASS
    if(is.null(sel.lag)) {
      sel.lag <- c(1:dim(A)[1])
    }
    if(is.null(sel.tau)) {
      sel.tau <- c(1:dim(A)[2])
    }
    if(sum(sel.lag < 1 | sel.lag > dim(A)[1]) > 0) {
      stop(paste0("sel.lag outside the range of lag index 1-",dim(A)[1]))	
    }
    if(sum(sel.tau < 1 | sel.tau > dim(A)[2]) > 0) {
      stop(paste0("sel.tau outside the range of quantile index 1-",dim(A)[2]))	
    }
    nlag <- length(sel.lag)
    ntau <- length(sel.tau)
    B <- dim(A.sim)[1]
    Sig <- matrix(0,nrow=ntau*nlag,ncol=ntau*nlag)
    for(b in c(1:B)) {
      Sig <- Sig + tcrossprod(c(A.sim[b,sel.lag,sel.tau]))
    } 
    Sig <- Sig / B
    Sig <- MASS::ginv(Sig)
    chi2 <- c(crossprod(c(A[sel.lag,sel.tau]),Sig) %*% as.vector(c(A[sel.lag,sel.tau])))
    chi2.H0 <- rep(0,B)
    for(b in c(1:B)) {
      chi2.H0[b] <- crossprod(c(A.sim[b,sel.lag,sel.tau]),Sig) %*% as.vector(c(A.sim[b,sel.lag,sel.tau]))
    }
    p.value <- length(which(chi2.H0 >= chi2))/B
    return(list(wald=chi2,p.value=p.value,Sig.inv=Sig,a=c(A[sel.lag,sel.tau])))
  }

  test <- wald_test(A,A.sim,sel.lag,sel.tau)
  A.u <- apply(A.sim,c(2,3),quantile,prob=0.975)
  A.l <- apply(A.sim,c(2,3),quantile,prob=0.025)
  return(list(test=test,A.u=A.u,A.l=A.l))
}



#' Bootstrap Simulation of SAR Coefficients for Testing Equality of Granger-Causality in Two Samples
#'
#' This function simulates bootstrap samples of selected spline autoregression (SAR) coefficients 
#' for testing equality of Granger-causality in two samples based on their SAR models
#' under H0: effect in each sample equals the average effect.
#' @param y.qser matrix or array of QSER from \code{qser()} or \code{qspec.sar()$qser}
#' @param fit object of SAR model from \code{qser2sar()} or \code{qspec.sar()$fit}
#' @param fit2 object of SAR model for the other sample
#' @param index a pair of component indices for multiple time series 
#' or a sequence of lags for single time series (default = \code{c(1,2)})
#' @param nsim number of bootstrap samples (default = 1000)
#' @param method method of residual calculation: \code{"ar"} (default) or \code{"sar"}
#' @param refit if \code{TRUE} the models are refit under H0 (default=\code{FALSE})
#' @param n.cores number of cores for parallel computing (default = 1)
#' @param mthreads if \code{FALSE} (default), set \code{RhpcBLASctl::blas_set_num_threads(1)} 
#' @param seed seed for random sampling (default = \code{1234567})
#' @return array of simulated bootstrap samples of selected SAR coefficients
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @import 'RhpcBLASctl'
#' @export
#' @examples
#' y11 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y21 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' y12 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y22 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y1.sar <- qspec.sar(cbind(y11,y21),tau0=tau,p=1)
#' y2.sar <- qspec.sar(cbind(y12,y22),tau0=tau,p=1)
#' A1.sim <- sar.eq.bootstrap(y1.sar$qser,y1.sar$fit,y2.sar$fit,index=c(1,2),nsim=5)
#' A2.sim <- sar.eq.bootstrap(y2.sar$qser,y2.sar$fit,y1.sar$fit,index=c(1,2),nsim=5)
sar.eq.bootstrap <- function(y.qser,fit,fit2,index=c(1,2),nsim=1000,method=c("ar","sar"),
                    refit=FALSE,n.cores=1,mthreads=FALSE,seed=1234567) {

  sarc2_refit <- function(fit1.sar,fit2.sar,index=c(1,2)) {
  # compute two-sample constrainted SAR coefficients from unconstraint fit1.sar and fit2.sar
  # where a constraint is imposed on the (index[1],index[2]) entry of the SAR coefficients
    Omega1 <- fit1.sar$fit$Omega
    Theta1 <- fit1.sar$fit$Theta
    Omega2 <- fit2.sar$fit$Omega
    Theta2 <- fit2.sar$fit$Theta
    p <- fit1.sar$fit$p
    sdm <- fit1.sar$sdm
    nc <- dim(Theta1)[1] 
    K <- dim(sdm$bsMat)[2]
    ntau <- dim(sdm$bsMat2)[1]
    e1 <- rep(0,nc)
    e1[index[1]] <- 1
    e2 <- rep(0,nc)
    e2[index[2]] <- 1
    Q <- kronecker(diag(K*p),t(e2))
    R1 <- Q %*% Omega1
    R2 <- Q %*% Omega2
    # compute Lagrange multipler
    xi <- solve((R1+R2) %*% t(Q))
    xi <- xi %*% (Q %*% t(Theta1-Theta2) %*% e1)
    # compute SARC 
    Theta1 <- Theta1 - kronecker(e1,t(xi) %*% R1)
    Theta2 <- Theta2 + kronecker(e1,t(xi) %*% R2)
    # compute SAR coefficient functions
    A1 <- fit1.sar$A
    A2 <- fit2.sar$A
    for(j in c(1:p)) {
      theta1 <- matrix(Theta1[,((j-1)*K*nc+1):(j*K*nc)],nrow=nc)
	  theta2 <- matrix(Theta2[,((j-1)*K*nc+1):(j*K*nc)],nrow=nc)
      for(ell in c(1:ntau)) {
	    Phi <- matrix(0,nrow=K*nc,ncol=nc)
        for(k in c(1:K)) {
          Phi[((k-1)*nc+1):(k*nc),] <- diag(sdm$bsMat2[ell,k],nrow=nc,ncol=nc)
        }
        A1[,,j,ell] <- theta1 %*% Phi
	    A2[,,j,ell] <- theta2 %*% Phi
      }
    }
    return(list(A.1=A1,A.2=A2))
  }

  simulate_qser_array_eq <- function(resid,sample.idx,fit,fit2,idx0,refit=FALSE) {
  # this function simulates quantile series under H0 from bootstrap samples of residuals
  # input: resid = nc*(n-p)*ntau-array of residuals from qser2ar()
  #        sample.idx = nn-vector of sampled time indices
  #        fit, fit2 = SAR models from qser2sar() for sample 1 and sample 2
  #        idx0 = indices of SAR coefficients to set to average under H0
  #        refit = if TRUE, SAR model is refit under the constraint H0
  # output: nc*n*ntau array of simulated quantile series
      p <- fit$p
      nc <- dim(resid)[1]
	  n <- dim(resid)[2] + p
      ntau <- dim(resid)[3]
      nn <- length(sample.idx)
	  A <- fit$A
      A2 <- fit2$A
	  if(refit) { 
	    fit.sarc <- sarc2_refit(fit,fit2,idx0)
		A <- fit.sarc$A.1
		A2 <- fit.sarc$A.2
	  }
      resid.sim <- array(resid[,sample.idx,],dim=c(nc,nn,ntau))
      qser.sim <- array(0,dim=c(nc,nn,ntau))
      for(ell in c(1:ntau)) {
	    mu <- apply(resid.sim[,,ell],1,mean) 
        for(jj in c((p+1):nn)) {
          if(nc == 1) {
            for(j in c(1:p)) {
              A0 <- A[j,ell]
              if(j %in% idx0) A0 <- 0.5*(A0 + A2[j,ell])
              A0 <- matrix(A0,ncol=nc,nrow=nc)			  
              qser.sim[,jj,ell] <- qser.sim[,jj,ell] + A0 %*% matrix(qser.sim[,jj-j,ell],nrow=nc)
            }		  
          } else {
            for(j in c(1:p)) {
              A0 <- A[,,j,ell]
              A0[idx0[1],idx0[2]] <- 0.5*(A0[idx0[1],idx0[2]] + A2[idx0[1],idx0[2],j,ell])
              qser.sim[,jj,ell] <- qser.sim[,jj,ell] + A0 %*% matrix(qser.sim[,jj-j,ell],nrow=nc)
            } 
          }
          qser.sim[,jj,ell] <- qser.sim[,jj,ell] + (resid.sim[,jj,ell] - mu)
        }
      }
      qser.sim <- qser.sim[,c((nn-n+1):nn),]
      return(qser.sim)
  }

  if(is.matrix(y.qser)) y.qser <- array(y.qser,dim=c(1,dim(y.qser)))  
  nc <- dim(y.qser)[1]
  n <- dim(y.qser)[2]
  ntau <- dim(y.qser)[3]
  p <- fit$p
  
  sar.gc.param.validity(index,nc,p) 
  
  # compute residauls from sample 1
  if(method[1]=="ar") {
    # from AR model fitted to QSER for each quantile
    resid <- qser2ar(y.qser,p=p)$resid
  } else {
    # from SAR model fitted to QSER
    resid <- sar.residual(y.qser,fit)
  }
  
  # generate bootstrap sampling of time stamps
  nn <- 2*n
  set.seed(seed)
  idx <- list()
  for(sim.idx in c(1:nsim)) {
    idx[[sim.idx]] <- sample(c(1:dim(resid)[2]),nn,replace=TRUE)
  }
  
  `%dopar%` <- foreach::`%dopar%`
  `%do%` <- foreach::`%do%`
  if(n.cores > 1) {
  	if(!mthreads) {
	  ##sink(tempfile())
      RhpcBLASctl::blas_set_num_threads(1)
	  ##sink()
    }
    cl <- parallel::makeCluster(n.cores)
    parallel::clusterExport(cl, c("qser2sar","sar.gc.coef"))
    doParallel::registerDoParallel(cl)
    out <- foreach::foreach(sim.idx=c(1:nsim)) %dopar% {
	  qser.sim <- simulate_qser_array_eq(resid,idx[[sim.idx]],fit,fit2,index,refit)
      fit.sim <- qser2sar(qser.sim,tau0=fit$tau0,tau=fit$tau,p=fit$p,spar=fit$spar,type=fit$type,weights=fit$weights)
      A.sim <- sar.gc.coef(fit.sim,index)
    }
    parallel::stopCluster(cl)
	if(!mthreads) {
      ##sink(tempfile())
      RhpcBLASctl::blas_set_num_threads(Sys.getenv('OPENBLAS_NUM_THREADS'))
      ##sink()
    }
  } else {
    out <- foreach::foreach(sim.idx=c(1:nsim)) %do% {
	  qser.sim <- simulate_qser_array_eq(resid,idx[[sim.idx]],fit,fit2,index,refit)
      fit.sim <- qser2sar(qser.sim,tau0=fit$tau0,tau=fit$tau,p=fit$p,spar=fit$spar,type=fit$type,weights=fit$weights)
      A.sim <- sar.gc.coef(fit.sim,index)
    }
  }
  A.sim <- array(0,dim=c(nsim,dim(out[[1]])))
  for(sim.idx in c(1:nsim)) {
    A.sim[sim.idx,,] <- out[[sim.idx]]
  }
  return(A.sim)
}
  

#' Wald Test and Confidence Band for Equality of Granger-Causality in Two Samples
#'
#' This function computes Wald test and confidence band for equality of Granger-causality in two samples
#' using bootstrap samples generated by \code{sar.eq.bootstrap()} based on the spline autoregression (SAR) models
#' of quantile series (QSER).
#' @param A1 matrix of selected SAR coefficients for sample 1
#' @param A1.sim simulated bootstrap samples from \code{sar.eq.bootstrap()} for sample 1
#' @param A2 matrix of selected SAR coefficients for sample 2
#' @param A2.sim simulated bootstrap samples from \code{sar.eq.bootstrap()} for sample 2
#' @param sel.lag indices of time lags for Wald test (default = \code{NULL}: all lags)
#' @param sel.tau indices of quantile levels for Wald test (default = \code{NULL}: all quantiles)
#' @return a list with the following elements:
#'   \item{test}{list of Wald test result containing \code{wald} and \code{p.value}}
#'   \item{D.u}{matrix of upper limits of 95\% confidence band for \code{A1 - A2}}
#'   \item{D.l}{matrix of lower limits of 95\% confidence band for \code{A1 - A2}}
#' @import 'MASS'
#' @export
#' @examples
#' y11 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y21 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' y12 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y22 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y1.sar <- qspec.sar(cbind(y11,y21),tau0=tau,p=1)
#' y2.sar <- qspec.sar(cbind(y12,y22),tau0=tau,p=1)
#' A1.sim <- sar.eq.bootstrap(y1.sar$qser,y1.sar$fit,y2.sar$fit,index=c(1,2),nsim=5)
#' A2.sim <- sar.eq.bootstrap(y2.sar$qser,y2.sar$fit,y1.sar$fit,index=c(1,2),nsim=5)
#' A1 <- sar.gc.coef(y1.sar$fit,index=c(1,2))
#' A2 <- sar.gc.coef(y2.sar$fit,index=c(1,2))
#' test <- sar.eq.test(A1,A1.sim,A2,A2.sim,sel.lag=NULL,sel.tau=NULL)
sar.eq.test <- function(A1,A1.sim,A2,A2.sim,sel.lag=NULL,sel.tau=NULL) {
   
  wald_eq_test <- function(A1,A1.sim,A2,A2.sim,sel.lag,sel.tau) {
  # this function computes the two-sample bootstrap Wald statistic and p-value
  # of selected elements in A and A2 using bootstrap samples in A.sim and A2.sim 
  # for testing the equality of A1 and A2, i.e., H0: (A1-A2)[selected elements] = 0
  # input:   A1, A2 = p*ntau-matrix of selected SAR(p) coefficients for all lags and quantiles
  #          A1.sim. A2.sim = nsim*p*ntau-array of bootstrap samples associated with A1 and A2 under H0
  #          sel.lap = selected lags for testing (NULL for all lags )
  #          sel.tau = selected quantile levels for testing (NULL for all quantiles)
  # imports: MASS
    if(is.null(sel.lag)) {
      sel.lag <- c(1:dim(A1)[1])
    }
    if(is.null(sel.tau)) {
      sel.tau <- c(1:dim(A1)[2])
    }
    if(sum(sel.lag < 1 | sel.lag > dim(A1)[1]) > 0) {
      stop(paste0("sel.lag outside the range of lag index 1-",dim(A1)[1]))	
    }
    if(sum(sel.tau < 1 | sel.tau > dim(A1)[2]) > 0) {
      stop(paste0("sel.tau outside the range of quantile index 1-",dim(A1)[2]))	
    }
    nlag <- length(sel.lag)
    ntau <- length(sel.tau)
    B <- dim(A1.sim)[1]
    Sig1 <- matrix(0,nrow=ntau*nlag,ncol=ntau*nlag)
    Sig2 <- matrix(0,nrow=ntau*nlag,ncol=ntau*nlag)
    for(b in c(1:B)) {
      Sig1 <- Sig1 + tcrossprod(c(A1.sim[b,sel.lag,sel.tau]))
      Sig2 <- Sig2 + tcrossprod(c(A2.sim[b,sel.lag,sel.tau]))	  
    } 
    Sig1 <- Sig1 / B
    Sig2 <- Sig2 / B	
    Sig <- MASS::ginv(Sig1+Sig2)
    delta <- c(A1[sel.lag,sel.tau]-A2[sel.lag,sel.tau])
    chi2 <- c(crossprod(delta,Sig) %*% as.vector(delta))
    chi2.H0 <- rep(0,B)
    for(b in c(1:B)) {
	  tmp <- c(A1.sim[b,sel.lag,sel.tau] - A2.sim[b,sel.lag,sel.tau])
      chi2.H0[b] <- crossprod(tmp,Sig) %*% as.vector(tmp)
    }
    p.value <- length(which(chi2.H0 >= chi2))/B
    return(list(wald=chi2,p.value=p.value,Sig.inv=Sig,D=delta))
  }
  
  test <- wald_eq_test(A1,A1.sim,A2,A2.sim,sel.lag,sel.tau)
  D.u <- apply(A1.sim-A2.sim,c(2,3),quantile,prob=0.975)
  D.l <- apply(A1.sim-A2.sim,c(2,3),quantile,prob=0.025)
  return(list(test=test,D.u=D.u,D.l=D.l))
}


#' Periodogram (PER)
#'
#' This function computes the periodogram or periodogram matrix for univariate or multivariate time series.
#' @param y vector or matrix of time series s (if matrix, nrow(y) = length of time series)

#' @return vector or array of periodogram
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y.per <- per(y)
#' plot(y.per)
per <- function(y) {
  if(!is.matrix(y)) y <- matrix(y,ncol=1,nrow=length(y))
  nc <- ncol(y)
  n <- nrow(y)
  y.dft <- matrix(0,nrow=n,ncol=nc)
  for(k in c(1:nc)) {
    y.dft[,k] <- fft(y[,k])
  }  
  y.per <- array(0,dim=c(nc,nc,n))
  for(k in c(1:nc)) {
    for(kk in c(k:nc)) {
	  y.per[k,kk,] <- y.dft[,k] * Conj(y.dft[,kk]) / n
      if(k != kk) y.per[kk,k,] = y.per[k,kk,]
	}
  }
  if(nc == 1) y.per <- Re(y.per[1,1,])
  return(y.per)
}


#' Quantile-Crossing Series (QCSER)
#'
#' This function creates the quantile-crossing series (QCSER) for univariate or multivariate time series.
#' @param y vector or matrix of time series
#' @param tau sequence of quantile levels in (0,1)
#' @param normalize \code{TRUE} or \code{FALSE} (default): normalize QCSER to have unit variance
#' @return A matrix or array of quantile-crossing series
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y.qser <- qcser(y,tau)
#' dim(y.qser)
qcser <- function(y,tau,normalize=FALSE) {
  if(!is.matrix(y)) y <- matrix(y,ncol=1,nrow=length(y))
  n <- nrow(y)
  nc <- ncol(y)
  nq <- length(tau)
  y.qcser <- array(0,dim=c(nc,n,nq))
  for(k in c(1:nc)) {
    q <- quantile(y[,k],tau)
    for(j in c(1:nq)) {
	  y.qcser[k,,j] <- 0
	  y.qcser[k,which(y[,k] <= q[j]),j] <- 1
	  y.qcser[k,,j] <- tau[j] - y.qcser[k,,j]
	  if(normalize) y.qcser[k,,j] <- y.qcser[k,,j] / sqrt(tau[j]*(1-tau[j]))
	}
  }
  if(nc==1) y.qcser <- matrix(y.qcser,ncol=nq)
  return(y.qcser)
}


#' ACF of Quantile Series (QSER) or Quantile-Crossing Series (QCACF)
#'
#' This function creates the ACF of quantile series or quantile-crossing series
#' @param y.qser matrix or array of quantile-crossing series
#' @return A matrix or array of ACF
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y.qser <- qcser(y,tau)
#' y.qacf <- qser2qacf(y.qser)
#' dim(y.qacf)
qser2qacf <- function(y.qser) {
  if(is.matrix(y.qser)) y.qser <- array(y.qser,dim=c(1,nrow(y.qser),ncol(y.qser)))
  nc <- dim(y.qser)[1]
  n <- dim(y.qser)[2]
  nq <- dim(y.qser)[3]
  y.qacf <- array(0,dim=c(nc,nc,n,nq))
  if(nc > 1) {
    for(j in c(1:nq)) {
      tmp <- stats::acf(t(as.matrix(y.qser[,,j],ncol=nc)),
	                    type="covariance",lag.max=n-1,plot=FALSE,demean=TRUE)$acf
      for(k in c(1:nc)) {
        for(kk in c(1:nc)) {
          y.qacf[k,kk,,j] <- tmp[,k,kk] 
	    }
      }
    }
    rm(tmp)
  } else {
    for(j in c(1:nq)) {
      tmp <- stats::acf(c(y.qser[,,j]),
	                    type="covariance",lag.max=n-1,plot=FALSE,demean=TRUE)$acf 
      y.qacf[1,1,,j] <- tmp[,1,1]
    }	 
  }
  if(nc==1) y.qacf <- matrix(y.qacf[1,1,,],ncol=nq)
  return(y.qacf)
}


# Smoothing of AR Residual Covariance Matrix Across Quantiles
#
# This function perfoms quantile smoothing of Residual Covariance matrix
# @param V array of functional residual covariance to be smoothed
# @param method quantile smoothing method: \code{"gamm"}, \code{"sp"}, or \code{"none"} (default)
# @param spar smoothing parameter for "sp" (default = NULL)
# @param type smoothing type: 1 (smooth V, default) or 2 (smooth square-root of V)
Vsmooth <- function(V,method=c("none","gamm","sp"),spar=NULL,type=c(1,2)) {
    nc <- dim(V)[1]
	ntau <- dim(V)[3]
    if(method[1] %in% c("gamm","sp")) {
  	  if(nc == 1) {
	    if(type[1]==2) {
		  V[1,1,] <- qsmooth(sqrt(V[1,1,]),method=method[1],spar=spar)
		  V[1,1,] <- V[1,1,]*V[1,1,]
		} else {
	      V[1,1,] <- qsmooth(V[1,1,],method=method[1],spar=spar)
		}
	  } else {
	    if(type[1]==2) {
  	      U <- V
	      for(ell in c(1:ntau)) {
	        tmp <- svd(V[,,ell])
		    U[,,ell] <- tmp$u %*% diag(sqrt(tmp$d))
          }
          for(k in c(1:nc)) {
	        for(kk in c(1:nc)) {   
	          U[k,kk,] <- qsmooth(U[k,kk,],method=method[1],spar=spar)
	        }
	      }
	      for(ell in c(1:ntau)) {
		    V[,,ell] <- U[,,ell] %*% t(U[,,ell])
          }
	    } else {
          for(k in c(1:nc)) {
	        for(kk in c(k:nc)) {   
	          V[k,kk,] <- qsmooth(V[k,kk,],method=method[1],spar=spar)
		      V[kk,k,] <- V[k,kk,]
	        }
	      }	   
	    }
	  }
    }
    return(V)
}	



# Smoothing of AR Coefficient Matrix Across Quantiles
#
# This function perfoms quantile smoothing of AR Coefficients
# @param A array of functional AR coefficients to be smoothed
# @param method quantile smoothing method: \code{"gamm"}, \code{"sp"}, or \code{"none"} (default)
# @param spar smoothing parameter for "sp" (default = "GCV")
Asmooth <- function(A,method=c("none","gamm","sp"),spar="GCV") {
    nc <- dim(A)[1]
    p <- dim(A)[3]
    if(method[1] %in% c("gamm","sp")) {
      for(k in c(1:nc)) {
	    for(kk in c(1:nc)) {
	      for(j in c(1:p)) {
	        A[k,kk,j,] <- qsmooth(A[k,kk,j,],method=method[1],spar=spar)
	      }
        }
	  }
    }
    return(A)
}



#' Autoregression (AR) Model of Quantile Series
#'
#' This function fits an autoregression (AR) model to quantile series (QSER) separately for each quantile level using \code{stats::ar()}.
#' @param y.qser matrix or array of pre-calculated QSER, e.g., using \code{qser()}
#' @param p order of AR model (default = \code{NULL}: selected by AIC)
#' @param order.max maximum order for AIC if \code{p = NULL} (default = \code{NULL}: determined by \code{stats::ar()})
#' @param method quantile smoothing method: \code{"gamm"} for \code{mgcv::gamm()}, 
#' \code{"sp"} for \code{stats::smooth.spline()}, or \code{"none"} (default)
#' @param spar smoothing parameter for \code{stats::smooth.spline()} (default = NULL for GCV)
#' @param type type of smoothing for covariance matrix: 1 (direct, default) or 2 (square root)
#' @return a list with the following elements:
#'   \item{A}{matrix or array of AR coefficients}
#'   \item{V}{vector or matrix of residual covariance}
#'   \item{p}{order of AR model}
#'   \item{n}{length of time series}
#'   \item{residuals}{matrix or array of residuals}
#' @import 'stats'
#' @import 'mgcv'
#' @export
qser2ar <- function(y.qser,p=NULL,order.max=NULL,method=c("none","gamm","sp"),spar=NULL,type=c(1,2)) {
  if(is.matrix(y.qser)) y.qser <- array(y.qser,dim=c(1,dim(y.qser)))
  nc <- dim(y.qser)[1]
  n <- dim(y.qser)[2]
  ntau <- dim(y.qser)[3]
  # order selection by aic
  if(is.null(p)) {
    aic <- NULL
    for(ell in c(1:ntau)) {
      mu <- apply(matrix(y.qser[,,ell],nrow=nc),1,mean)
      aic <- cbind(aic,stats::ar(t(matrix(y.qser[,,ell]-mu,nrow=nc)),order.max=order.max,method=c("yule-walker"))$aic)
    }
    aic <- apply(aic,1,mean)
    # optimal order minimizes aic
    p <- as.numeric(which.min(aic)) - 1
  }
  # p = 0
  if(p == 0) {
     A <- NULL	
     V <- array(0,dim=c(nc,nc,ntau))
     resid <- array(0,dim=c(nc,n,ntau))
     for(ell in c(1:ntau)) {
       mu <- apply(matrix(y.qser[,,ell],nrow=nc),1,mean)
       V[,,ell] <- matrix(tcrossprod(matrix(y.qser[,,ell]-mu,nrow=nc)),ncol=nc)/n
       resid[,,ell] <- matrix(y.qser[,,ell]-mu,nrow=nc)
     }
     V <- Vsmooth(V,method[1],spar=spar,type=type[1])	 
     if(nc==1) {
       V <- V[1,1,]
       resid <- resid[1,,]
     }
     return(list(A=A,V=V,p=p,n=n,residuals=y.qser))
  }
  # p > 0
  V <- array(0,dim=c(nc,nc,ntau))
  A <- array(0,dim=c(nc,nc,p,ntau))
  resid <- array(0,dim=c(nc,n-p,ntau))
  for(ell in c(1:ntau)) {
    mu <- apply(matrix(y.qser[,,ell],nrow=nc),1,mean)
    fit <- stats::ar(t(matrix(y.qser[,,ell]-mu,nrow=nc)), aic=FALSE, order.max = p, method=c("yule-walker"))
    if(nc==1) {
      V[,,ell] <- sum(fit$resid[-c(1:p)]^2)/n
      for(j in c(1:p)) A[,,j,ell] <- fit$ar[j]	
      resid[,,ell] <- t(fit$resid[-c(1:p)])
    } else {
      V[,,ell] <- matrix(crossprod(fit$resid[-c(1:p),]),ncol=nc)/n
      for(j in c(1:p)) A[,,j,ell] <- fit$ar[j,,]
      resid[,,ell] <- t(fit$resid[-c(1:p),])	  
    }	
  }
  A <- Asmooth(A,method[1],spar=spar)
  V <- Vsmooth(V,method[1],spar=spar,type=type[1])
  if(nc==1) {
    V <- V[1,1,]
    A <- matrix(A[1,1,,],nrow=p)
    resid <- resid[1,,]
  }
  return(list(A=A,V=V,p=p,n=n,residuals=resid))
}



#' Autoregression (AR) Estimator of Quantile Spectrum
#'
#' This function computes autoregression (AR) estimate of quantile spectrum from time series or quantile series (QSER).
#' @param y vector or matrix of time series (if matrix, \code{nrow(y)} = length of time series)
#' @param tau sequence of quantile levels in (0,1)
#' @param y.qser matrix or array of pre-calculated QSER (default = \code{NULL}: compute from \code{y} and \code{tau});
#' if \code{y.qser} is supplied, \code{y} and \code{tau} can be left unspecified
#' @param p order of AR model (default = \code{NULL}: automatically selected by AIC)
#' @param order.max maximum order for AIC if \code{p = NULL} (default = \code{NULL}: determined by \code{stats::ar()})
#' @param freq sequence of frequencies in [0,1) (default = \code{NULL}: all Fourier frequencies)
#' @param method method of quantile smoothing: \code{"gamm"} for \code{mgcv::gamm()}, 
#' \code{"sp"} for \code{stats::smooth.spline()}, or \code{"none"} (default)
#' @param spar smoothing parameter for \code{stats::smooth.spline()} (default = NULL for GCV)
#' @param type type of smoothing for covariance matrix: 1 (direct, default) or 2 (square root)
#' @param n.cores number of cores for parallel computing of QDFT if \code{y.qser = NULL} (default = 1)
#' @param cl pre-existing cluster for repeated parallel computing of QDFT (default = \code{NULL})
#' @return a list with the following elements:
#'   \item{spec}{matrix or array of AR quantile spectrum}
#'   \item{freq}{sequence of frequencies}
#'   \item{fit}{object of AR model}
#'   \item{qser}{matrix or array of quantile series if \code{y.qser = NULL}}
#' @import 'stats'
#' @import 'mgcv'
#' @export
#' @examples
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' y <- cbind(y1,y2)
#' tau <- seq(0.1,0.9,0.05)
#' n <- length(y1)
#' ff <- c(0:(n-1))/n
#' sel.f <- which(ff > 0 & ff < 0.5)
#' y.qspec.ar <- qspec.ar(y,tau,p=1)$spec
#' qfa.plot(ff[sel.f],tau,Re(y.qspec.ar[1,1,sel.f,]))
#' y.qser <- qcser(y1,tau)
#' y.qspec.ar <- qspec.ar(y.qser=y.qser,p=1)$spec
#' qfa.plot(ff[sel.f],tau,Re(y.qspec.ar[sel.f,]))
#' y.qspec.arqs <- qspec.ar(y.qser=y.qser,p=1,method="sp")$spec
#' qfa.plot(ff[sel.f],tau,Re(y.qspec.arqs[sel.f,]))
qspec.ar <- function(y,tau,y.qser=NULL,p=NULL,order.max=NULL,freq=NULL,
                method=c("none","gamm","sp"),spar=NULL,type=c(1,2),n.cores=1,cl=NULL) {
  return.qser <- FALSE
  if(is.null(y.qser)) {
    y.qser <- qser(y,tau,n.cores=n.cores,cl=cl)
	return.qser <- TRUE
  }
  fit <- qser2ar(y.qser,p=p,order.max=order.max,method=method[1],spar=spar,type=type[1])
  tmp <- ar2qspec(fit,freq)
  if(return.qser) {
    return(list(spec=tmp$spec,freq=tmp$freq,fit=fit,qser=y.qser)) 
  } else {
    return(list(spec=tmp$spec,freq=tmp$freq,fit=fit))
  }
}




#' Lag-Window (LW) Estimator of Quantile Spectrum
#'
#' This function computes lag-window (LW) estimate of quantile spectrum
#' with or without quantile smoothing from time series or quantile autocovariance function (QACF).
#' @param y vector or matrix of time series (if matrix, \code{nrow(y)} = length of time series)
#' @param tau sequence of quantile levels in (0,1)
#' @param y.qacf matrix or array of pre-calculated QACF (default = \code{NULL}: compute from \code{y} and \code{tau});
#' if \code{y.qacf} is supplied, \code{y} and \code{tau} can be left unspecified
#' @param M bandwidth parameter of lag window (default = \code{NULL}: quantile periodogram)
#' @param method quantile smoothing method:  \code{"gamm"} for \code{mgcv::gamm()}, 
#' \code{"sp"} for \code{stats::smooth.spline()}, or \code{"none"} (default)
#' @param spar smoothing parameter in \code{smooth.spline()} if \code{method = "sp"} (default = \code{"GCV"})
#' @param n.cores number of cores for parallel computing (default = 1)
#' @param cl pre-existing cluster for repeated parallel computing (default = \code{NULL})
#' @return A list with the following elements:
#'   \item{spec}{matrix or array of spectral estimate}
#'   \item{spec.lw}{matrix or array of spectral estimate without quantile smoothing}
#'   \item{lw}{lag-window sequence}
#'   \item{qacf}{matrix or array of quantile autocovariance function if \code{y.qacf = NULL}}
#' @import 'stats'
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @import 'mgcv'
#' @import 'nlme'
#' @export
#' @examples
#' y1 <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' y2 <- stats::arima.sim(list(order=c(1,0,0), ar=-0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' n <- length(y1)
#' ff <- c(0:(n-1))/n
#' sel.f <- which(ff > 0 & ff < 0.5)
#' y.qacf <- qacf(cbind(y1,y2),tau)
#' y.qper.lw <- qspec.lw(y.qacf=y.qacf,M=5)$spec
#' qfa.plot(ff[sel.f],tau,Re(y.qper.lw[1,1,sel.f,]))
#' y.qper.lwqs <- qspec.lw(y.qacf=y.qacf,M=5,method="sp",spar=0.9)$spec
#' qfa.plot(ff[sel.f],tau,Re(y.qper.lwqs[1,1,sel.f,]))
qspec.lw <- function(y,tau,y.qacf=NULL,M=NULL,method=c("none","gamm","sp"),spar="GCV",n.cores=1,cl=NULL) {
  return.qacf <- FALSE
  if(is.null(y.qacf)) {
    y.qacf <- qacf(y,tau,n.cores=n.cores,cl=cl)
	return.qacf <- TRUE
  }
  y.lw <- qacf2speclw(y.qacf,M)
  if(method[1] %in% c("sp","gamm")) {
    y.qper.lwqs <- qsmooth.qper(y.lw$spec,method=method[1],spar=spar,n.cores=n.cores,cl=cl)
  } else {
    y.qper.lwqs <- y.lw$spec
  }
  if(return.qacf) {
    return(list(spec=y.qper.lwqs,spec.lw=y.lw$spec,lw=y.lw$lw,qacf=y.qacf))
  } else {
    return(list(spec=y.qper.lwqs,spec.lw=y.lw$spec,lw=y.lw$lw))
  }  
}




optim.adam <- function(theta0,fn,gr,...,sg.rate=c(1,1),gr.m=NULL,gr.v=NULL,
  control=list()) {
# The ADAM method to find the minimizer fn - allow stochastic gradient and scheduled stepsize reduction
# input:  theta0 = vector of initial values of parameter
#             fn = objective function
#             gr = gradient function
#            ... = additional parameters of fn and gr
#           gr.m = vector of initial values of mean gradient (same dimension as theta0)
#           gr.v = vector of initial values of mean squared gradient (same dimension as theta0)
#        sg.rate = vector of sampling rates for quantiles and observations in stochastic gradient calculation
# output:    par = optimum parameter
#    convergence = convergence indicator (0 = converged; 1 = maxit is reseached before convergence)
#          value = minimum value of the objective function

  con <- list(maxit = 100, stepsize=0.01, warmup=70, stepupdate=20, stepredn=0.2, seed=1000,trace=0,
  reltol=sqrt(.Machine$double.eps),converge.test=F)
  con[(namc <- names(control))] <- control
  
  maxit <- con$maxit
  reltol <- con$reltol
  stepsize <- con$stepsize
  warmup <- con$warmup
  stepupdate <- con$stepupdate
  stepredn <- con$stepredn
  converge.test <- con$converge.test
  trace <- con$trace
  seed <- con$seed
  
  fn2 <- function(par) fn(par, ...)
  gr2 <- function(par,sg.rate) gr(par, ..., sg.rate)
  
  ptm <- proc.time()

  set.seed(seed)
  
  b1 <- 0.9
  b2 <- 0.999
  eps <- 1e-8
  if(is.null(gr.m)) {
    m <- rep(0,length(theta0))
  } else {
	m <- gr.m
  }
  if(is.null(gr.v)) {
    v <- rep(0,length(theta0))
  } else {
	v <- gr.v
  }
  theta <- theta0
  if(converge.test) cost <- fn2(theta)
  if(trace == -1) { 
    pars <- matrix(NA,nrow=length(theta),ncol=maxit)
    times <- rep(NA,maxit)
    vals <- rep(NA,maxit)
    steps <- rep(NA,maxit)
  }
  convergence <- 1
  i <- 1
  s <- stepsize
  while(i <= maxit & convergence == 1) {
    if( i > warmup ) {
      if(stepupdate > 0) {
        if(i %% stepupdate == 0 ) {
          s <- s * stepredn
        }
      }
    }
    g <- gr2(theta,sg.rate)
    m <- b1 * m	+ (1-b1) * g
    v <- b2 * v + (1-b2) * g^2
    d <- m / (1 - b1^i)
    d <- d / (sqrt(v / (1 - b2^i)) + eps)
    theta <- theta - s * d
    if(trace == -1) { 
      pars[,i] <- theta
      times[i] <- (proc.time()-ptm)[3]
      vals[i] <- fn2(theta)
	  steps[i] <- s
    }
    if(converge.test) {
      cost2 <- fn2(theta)	
      if(abs(cost2-cost) < reltol*(abs(cost)+reltol)) convergence <- 0
      cost <- cost2
    }
    i <- i + 1
  }
  if(!converge.test) cost <- fn2(theta)
  if(trace == -1) {
    return(list(par=theta,convergence=convergence,value=cost,counts=i-1,gr.m=m,gr.v=v,
                vals=vals,all.par=pars,all.time=times,steps=steps)) 
  } else {
    return(list(par=theta,convergence=convergence,value=cost,counts=i-1,gr.m=m,gr.v=v))
  }
}



optim.grad <- function(theta0,fn,gr,...,sg.rate=c(1,1),gr.m=NULL,gr.v=NULL,
  control=list()) {
# enhanced ADAM with limited line search to minimize fn: allow stochastic gradient
# input:  theta0 = vector of initial values of parameter
#             fn = objective function
#             gr = gradient function
#            ... = additional parameters of fn and gr
#           gr.m = vector of initial values of mean gradient (same dimension as theta0)
#           gr.v = vector of initial values of mean squared gradient (same dimension as theta0)
#        sg.rate = vector of sampling rates for quantiles and observations in stochastic gradient calculation
# output:    par = optimum parameter
#    convergence = convergence indicator (0 = converged; 1 = maxit is reseached before convergence)
#          value = minimum value of the objective function

  con <- list(maxit = 100, stepsize=0.01, warmup=70, stepupdate=20, stepredn=0.2, line.search.type=1,
  seed=1000,trace=0,reltol=sqrt(.Machine$double.eps),acctol=1e-04, line.search.max=5, converge.test=F)
  con[(namc <- names(control))] <- control
  
  maxit <- con$maxit
  warmup <- con$warmup
  reltol <- con$reltol
  stepsize <- con$stepsize
  converge.test <- con$converge.test
  trace <- con$trace
  acctol <- con$acctol
  stepredn <- con$stepredn
  stepupdate <- con$stepupdate
  count.max <- con$line.search.max
  line.search.type <- con$line.search.type[1]
  seed <- con$seed
  
  fn2 <- function(par) fn(par, ...)
  gr2 <- function(par,sg.rate) gr(par, ..., sg.rate)
  
  ptm <- proc.time()
  
  set.seed(seed)
  
  b1 <- 0.9
  b2 <- 0.999
  eps <- 1e-8
  if(is.null(gr.m)) {
    m <- rep(0,length(theta0))
  } else {
	m <- gr.m
  }
  if(is.null(gr.v)) {
    v <- rep(0,length(theta0))
  } else {
    v <- gr.v  
  }
  theta <- theta0
  if(converge.test) cost <- fn2(theta)
  if(trace == -1) {
    pars <- matrix(NA,nrow=length(theta),ncol=maxit)
    times <- rep(NA,maxit)
    vals <- rep(NA,maxit)
    steps <- rep(NA,maxit)
  }
  convergence <- 1
  i <- 1
  s <- stepsize
  while(i <= maxit & convergence == 1) {
    g <- gr2(theta,sg.rate)
    m <- b1 * m	+ (1-b1) * g
    v <- b2 * v + (1-b2) * g^2
    d <- m / (1 - b1^i)
    d <- d / (sqrt(v / (1 - b2^i)) + eps)
    if(i > warmup) {
      if(stepupdate > 0) {
	    if(i %% stepupdate == 0) {
          gradproj <- sum(d * d)
          if(line.search.type %in% c(1,2)) {
            # type=1 or 2: start with default step size
			steplength <- min(1, stepsize / stepredn^(floor(count.max/2)))
            ##steplength <- stepsize
          } else {
            # type=3 or 4: start with current step size
			steplength <- min(1, s / stepredn^(floor(count.max/2)))
            ##steplength <- s
          }
          accpoint <- F
          fmin <- fn2(theta)
          count <- 0
          while(!accpoint & count < count.max) {
            f <- fn2(theta - steplength * d)
            count <- count + 1
            accpoint <- (f <= fmin - gradproj * steplength * acctol)  
            if(!accpoint) {
              steplength <- steplength * stepredn
            }
          }
          if(accpoint) {
            s <- steplength
          } else {
            if(line.search.type %in% c(1,4)) {
              # type=1 or 4: fallback to default
              s <- stepsize		
            } else {
              # type=2 or 3: fallback to reduced default 
              s <- stepsize * stepredn
            }
          }
        }
      }
    }
    theta <- theta - s * d
    if(trace == -1) {
      pars[,i] <- theta
      times[i] <- (proc.time()-ptm)[3]
      vals[i] <- fn2(theta)
      steps[i] <- s
    }
    if(converge.test) {
      cost2 <- fn2(theta)	
      if(abs(cost2-cost) < reltol*(abs(cost)+reltol)) convergence <- 0
      cost <- cost2
    }
    i <- i + 1
  }
  if(!converge.test) cost <- fn2(theta)
  if(trace == -1) {
    return(list(par=theta,convergence=convergence,value=cost,counts=i-1,gr.m=m,gr.v=v,stepsize=s,
                steps=steps,vals=vals,all.par=pars,all.time=times))
  } else {
    return(list(par=theta,convergence=convergence,value=cost,counts=i-1,gr.m=m,gr.v=v,stepsize=s))
  }
}



sqr.optim_obj <- function(theta,y,X,tau0,sdm,spar=0,weighted=FALSE) {
  # objective function of spline quantile regression
  # input:  theta = p*K-vector of spline coefficients
  #                 where p = number of regression coefficients
  #                     K = number of spline basis functions
  #             y = n-vector of time series
  #             X = n*p regression design matrix
  #          tau0 = L-vector of selected quantile levels
  #           sdm = objec from create_spline_matrix
  #          spar = user-specified smoothing parameter
  #      weighted = T if penalty term gets quantile-dependent weights
  # output: cost = value of the objecitive function


  create_weights <- function(tau) {
  # quantile-dependent weights of second derivatives in SQR penalty
    0.25 / (tau*(1-tau))
  }

  rescale_penalty <- function(spar,n,tau0,X,sdm,weighted=FALSE) {
  # rescale penalty par by weights and weighted l1 norm of dbsMat
  # input: spar = user-specified penalty parameter
  #           n = number of observations
  #        tau0 = L vector of quantile levels
  #           X = design matrix
  #         sdm = object from create_spline_matrix()
  #    weighted = TRUE if penalty is weighted
  # output:  c2 = n * c * w
  # dependencies: create_weights(), create_Phi_matrix()
    L <- length(tau0)
	p <- ncol(X)
	n <- nrow(X)
    if(weighted) {
      w <- create_weights(tau0)
    } else {
      w <- rep(1,L)  
    }
	r <- 0
	for(ell in c(1:L)) {
      # design matrix for spline coefficients
      X0 <- create_Phi_matrix(sdm$bsMat[ell,],p)
      X0 <- X %*% X0
      r <- r + sum(abs(X0))
    }
	r <- r / sum(abs(w * sdm$dbsMat))
    c2 <- w * r * 1000.0**(spar - 1.0)
    c2
  }

  create_Phi_matrix <- function(bsVec,p) {
  # create spline basis matrix for LP and dual LP
  # input: bsVec = K-vector of bs functions (K=number of basis functions)
  #            p = number of parameters
  # output:  Phi = p*pK matrix of basis functions
    K <- length(bsVec)
    Phi <- matrix(0,nrow=p,ncol=p*K)
    for(i in c(1:p)) Phi[i,((i-1)*K+1):(i*K)] <- bsVec
    Phi
  }
  
  Rho <- function(u, tau) u * (tau - (u < 0))
 
  n <- length(y)
  L <- length(tau0)
  p <- ncol(X)

  # rescale penalty parameter
  c2 <- rescale_penalty(spar,n,tau0,X,sdm,weighted)
  
  theta <- matrix(theta,ncol=1)
  
  cost <- 0
  for(ell in c(1:L)) {
    # regression design matrix for spline coefficients
    Phi <- create_Phi_matrix(sdm$bsMat[ell,],p)
    cost <- cost + sum( Rho(y - c(X %*% (Phi %*% theta)),tau0[ell]) )
    # design matrix for penalty function
    Z0 <- create_Phi_matrix(sdm$dbsMat[ell,],p)
    cost <- cost + c2[ell] * sum( abs( Z0 %*% theta ) )
  }
  cost
}



sqr.optim_dobj <- function(theta,y,X,tau0,sdm,spar=0,weighted=FALSE,sg.rate=c(1,1),eps=1e-8) {
  # derivative of the objective function of spline quantile regression
  # input:  theta = p*K-vector of spline coefficients
  #                 where p = number of regression coefficients
  #                     K = number of spline basis functions
  #            y = n-vector of time series
  #            X = n*p regression design matrix
  #         tau0 = L-vector of selected quantile levels
  #          sdm = object from create_spline_matrix()
  #         spar = user-specified smoothing parameter
  #     weighted = T if penalty term gets quantile-dependent weights
  #      sg.rate = vector of sampling rates for quantiles and observations
  # output: dcost = value of the derivative of the objecitive function
  
  drho.rq <- function(y,tau,eps) {
    # derivative of check function for quantile regression
    # input: y = n-vector of dependent variables
    #        tau = quantile level
    # output: drho = n-vector of drho(y) evaluated at tau
    drho <- y
	drho[y>0] <- tau
    drho[y<0] <- -(1-tau)
	drho[abs(y) < eps] <- 0
	drho
  }

  dabs <- function(y,eps) {
    # derivative of absolute value
    dabs <- y
    dabs[y>0] <- 1
    dabs[y<0] <- -1
    dabs[abs(y) < eps] <- 0
    dabs
  }
  
  create_weights <- function(tau) {
  # quantile-dependent weights of second derivatives in SQR penalty
    0.25 / (tau*(1-tau))
  }

  rescale_penalty <- function(spar,n,tau0,X,sdm,weighted=FALSE) {
  # rescale penalty par by weights and weighted l1 norm of dbsMat
  # input: spar = user-specified penalty parameter
  #           n = number of observations
  #        tau0 = L vector of quantile levels
  #           X = design matrix
  #         sdm = object from create_spline_matrix()
  #    weighted = TRUE if penalty is weighted
  # output:  c2 = n * c * w
  # dependencies: create_weights(), create_Phi_matrix()
    L <- length(tau0)
	p <- ncol(X)
	n <- nrow(X)
    if(weighted) {
      w <- create_weights(tau0)
    } else {
      w <- rep(1,L)  
    }
	r <- 0
	for(ell in c(1:L)) {
      # design matrix for spline coefficients
      X0 <- create_Phi_matrix(sdm$bsMat[ell,],p)
      X0 <- X %*% X0
      r <- r + sum(abs(X0))
    }
	r <- r / sum(abs(w * sdm$dbsMat))
    c2 <- w * r * 1000.0**(spar - 1.0)
    c2
  }

  create_Phi_matrix <- function(bsVec,p) {
  # create spline basis matrix for LP and dual LP
  # input: bsVec = K-vector of bs functions (K=number of basis functions)
  #            p = number of parameters
  # output:  Phi = p*pK matrix of basis functions
    K <- length(bsVec)
    Phi <- matrix(0,nrow=p,ncol=p*K)
    for(i in c(1:p)) Phi[i,((i-1)*K+1):(i*K)] <- bsVec
    Phi
  }
  
  n <- length(y)
  L <- length(tau0)
  p <- ncol(X)

  # rescale penalty parameter
  c2 <- rescale_penalty(spar,n,tau0,X,sdm,weighted)
	
  theta <- matrix(theta,ncol=1)
  
  # sample the quantiles (min=2)
  sel.L <- c(1:L)
  if(sg.rate[1] < 1 && sg.rate[1] > 0) {
    sel.L <- sample(c(1:L),max(2,floor(sg.rate[1] * L)))
  }
  # sample the observations (min=10)
  sel.n <- c(1:n)
  if(sg.rate[2] < 1 && sg.rate[2] > 0) {
    sel.n <- sample(c(1:n),max(10,floor(sg.rate[2] * n)))
  }
  
  dcost <- rep(0,nrow(theta))
  for(ell in sel.L) {
    # regression design matrix for spline coefficients
    X1 <- create_Phi_matrix(sdm$bsMat[ell,],p)
	X1 <- X[sel.n,] %*% X1
    tmp <- matrix(drho.rq(y[sel.n] - c(X1 %*% theta),tau0[ell],eps),ncol=1)
    tmp <- -crossprod(X1,tmp)
    # rescale for observations
    tmp <- tmp * (n / length(sel.n))
    dcost <- dcost + tmp 
    # design matrix for penalty function
    Z0 <- create_Phi_matrix(sdm$dbsMat[ell,],p)
    tmp <- matrix(dabs( c(Z0 %*% theta),eps ),ncol=1)
    tmp <- crossprod(Z0,tmp)
    dcost <- dcost + c2[ell] * tmp
  }
  # rescale to offset the effect of sampling
  dcost <- dcost * (L / length(sel.L))
  dcost
}



#' Cubic Spline Quantile Regression with L1-Norm Roughness Penalty (SQR) Computed by Gradient Algorithms
#'
#' This function computes spline quantile regression with cubic splines and L1-norm roughness penalty by a gradient algorithm BFGS, ADAM, or GRAD.
#' @param X vecor or matrix of explanatory variables (including intercept)
#' @param y vector of dependent variable
#' @param tau sequence of quantile levels in (0,1)
#' @param spar smoothing parameter
#' @param d subsampling rate of quantile levels (default = 1)
#' @param weighted if \code{TRUE}, penalty function is weighted (default = \code{FALSE})
#' @param method optimization method: \code{"BFGS"} (default), \code{"ADAM"}, or \code{"GRAD"}
#' @param beta.rq matrix of regression coefficients from \code{quantreg::rq(y~X)} for initialization (default = \code{NULL})
#' @param theta0 initial value of spline coefficients (default = \code{NULL})
#' @param spar0 smoothing parameter for \code{stats::smooth.spline()} to smooth \code{beta.rq} for initilaiztion (default = \code{NULL})
#' @param sg.rate vector of sampling rates for quantiles and observations in stochastic gradient version of GRAD and ADAM
#' @param mthreads if \code{FALSE}, set \code{RhpcBLASctl::blas_set_num_threads(1)} (default = \code{TRUE})
#' @param control list of control parameters
#' \describe{
#'   \item{\code{maxit}:}{max number of iterations (default = 100)}
#'   \item{\code{stepsize}:}{stepsize for ADAM and GRAD (default = 0.01)}
#'   \item{\code{warmup}:}{length of warmup phase for ADAM and GRAD (default = 70)}
#'   \item{\code{stepupdate}:}{frequency of update for ADAM and GRAD (default = 20)}
#'   \item{\code{stepredn}:}{stepsize discount factor for ADAM and GRAD (default = 0.2)}
#'   \item{\code{line.search.type}:}{line search option (1,2,3,4) for GRAD (default = 1)}
#'   \item{\code{line.search.max}:}{max number of line search trials for GRAD (default = 1)}
#'   \item{\code{seed}:}{seed for stochastic version of ADAM and GRAD (default = 1000)}
#'   \item{\code{trace}:}{-1 return results from all iterations, 0 (default) return final result}
#' }
#' @return A list with the following elements:
#'   \item{beta}{matrix of regression coefficients}
#'   \item{all.beta}{coefficients from all iterations for GRAD and ADAM}
#'   \item{spars}{smoothing parameters from \code{stats::smooth.spline()} for initialization}
#'   \item{fit}{object from the optimization algorithm}
#' @import 'stats'
#' @import 'splines'
#' @import 'RhpcBLASctl'
#' @export
#' @examples
#' data(engel)
#' y <- engel$foodexp
#' X <- cbind(rep(1,length(y)),engel$income-mean(engel$income))
#' tau <- seq(0.1,0.9,0.05)
#' fit.rq <- quantreg::rq(y ~ X[,2],tau)
#' fit.sqr <- sqr(y ~ X[,2],tau,spar=0.2)
#' fit <- sqr.fit.optim(X,y,tau,spar=0.2,d=2,method="BFSG",beta.rq=fit.rq$coef)
#' fit <- sqr.fit.optim(X,y,tau,spar=0.2,d=2,method="BFSG",beta.rq=fit.rq$coef)
#' par(mfrow=c(1,2),pty="m",lab=c(10,10,2),mar=c(4,4,2,1)+0.1,las=1)
#' for(j in c(1:2)) {
#'   plot(tau,fit.rq$coef[j,],type="n",xlab="QUANTILE LEVEL",ylab=paste0("COEFF",j))
#'   points(tau,fit.rq$coef[j,],pch=1,cex=0.5)
#'   lines(tau,fit.sqr$coef[j,],lty=1); lines(tau,fit$beta[j,],lty=2,col=2)
#' }
sqr.fit.optim <- function(X,y,tau,spar=0,d=1,weighted=FALSE,method=c("BFGS","ADAM","GRAD"),
  beta.rq=NULL,theta0=NULL,spar0=NULL,sg.rate=c(1,1),mthreads=TRUE,control=list(trace=0)) {

  create_coarse_grid <- function(tau,d=1) {
  # create index of a coarse quantile grid for SQR
    ntau <- length(tau)
    sel.tau0 <- seq(1,ntau,d)
    if(sel.tau0[length(sel.tau0)] < ntau) sel.tau0 <- c(sel.tau0,ntau)
    sel.tau0
  }

  create_spline_matrix <- function(tau0,tau) {
  # create spline matrix and its second derivative for SQR
  # input:    tau0 = subset of tau
  #            tau = complete set of quantiles for interpolation
  # output: bsMat  = matrix of basis functions
  #         dbsMat = matrix of second derivatives
  #         bsMat2 = maxtrix of basis function for interpolation to tau
  # imports: 'splines'
    knots <- stats::smooth.spline(tau0,tau0)$fit$knot 
    # rescale tau0 and tau to [0,1] for spline matrix to be consistent with smooth.spline
    bsMat <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)))
    dbsMat <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)),derivs=2)
    bsMat2 <- splines::splineDesign(knots,(tau-min(tau))/(max(tau)-min(tau)))
    return(list(bsMat=bsMat,dbsMat=dbsMat,bsMat2=bsMat2,knots=knots))
  }

  create_Phi_matrix <- function(bsVec,p) {
  # create spline basis matrix for LP and dual LP
  # input: bsVec = K-vector of bs functions (K=number of basis functions)
  #            p = number of parameters
  # output:  Phi = p*pK matrix of basis functions
    K <- length(bsVec)
    Phi <- matrix(0,nrow=p,ncol=p*K)
    for(i in c(1:p)) Phi[i,((i-1)*K+1):(i*K)] <- bsVec
    Phi
  }
  
  # gr.m = vector of initial value of mean gradient for 'grad' and 'adam'
  # gr.v = vector of initial value of mean squared gradient for 'grad' and 'adam'  
  gr.m <- NULL
  gr.v <- NULL
  return.beta0 <- FALSE

  con <- list()
  con[(namc <- names(control))] <- control
  method <- method[1]
  if(!(method %in% c('GRAD'))) {
    con$warmup <- NULL
    con$line.search <- NULL
    con$line.search.max <- NULL
    con$line.search.type <- NULL
  }
  if(!(method %in% c('ADAM','GRAD'))) {
    con$seed <- NULL
  }
  
  trace <- con$trace
  if(is.null(trace)) trace <- 0
  
  # turn-off blas multithread to run in parallel
  if(!mthreads) {
    ##sink(tempfile())
    RhpcBLASctl::blas_set_num_threads(1)
    ##sink()
  }
  
  n <- length(y)
  ntau <-length(tau)
  p <- ncol(X)
  # create a coarse quantile grid for LP
  sel.tau0 <- create_coarse_grid(tau,d)
  tau0 <- tau[sel.tau0]
  L <- length(tau0)
  # create spline design matrix with knots provided by smooth.spline
  sdm <- create_spline_matrix(tau0,tau)
  K <- ncol(sdm$bsMat)
  
  # compute intial value of theta from spline smoothing of rq solution
  spars <- NULL  
  if(is.null(theta0)) {
    if(is.null(beta.rq)) {
      theta0 <- rep(0,p*K)
    } else {
      theta0 <- NULL
      for(i in c(1:p)) {
	    fit0 <- stats::smooth.spline(tau0,beta.rq[i,sel.tau0],spar=spar0)
        theta0 <- c(theta0, fit0$fit$coef)
	    spars <- c(spars,fit0$spar)
      }
      rm(fit0)
    }
  }
  
  # compute optimal solution theta
  if( !(method %in% c("BFGS","ADAM","GRAD")) ) method <- "BFGS"
  if(method=="ADAM") {
    fit <- optim.adam(theta0,sqr.optim_obj,sqr.optim_dobj,y=y,X=X,tau0=tau0,sdm=sdm,spar=spar,
               weighted=weighted,sg.rate=sg.rate,gr.m=gr.m,gr.v=gr.v,control=con)
  } 
  if (method=="GRAD") {
    fit <- optim.grad(theta0,sqr.optim_obj,sqr.optim_dobj,y=y,X=X,tau0=tau0,sdm=sdm,spar=spar,
               weighted=weighted,sg.rate=sg.rate,gr.m=gr.m,gr.v=gr.v,control=con)  
  }
  if (method=="BFGS") {
    con$trace <- max(0,con$trace)
    fit <- stats::optim(theta0,sqr.optim_obj,gr=sqr.optim_dobj,method="BFGS",y=y,X=X,tau0=tau0,sdm=sdm,spar=spar,
               weighted=weighted,control=con)
  }
  
  theta <- fit$par
  
  # map to interpolated regression solution
  beta <- matrix(0,nrow=p,ncol=ntau)
  for(ell in c(1:ntau)) {
    Phi <- create_Phi_matrix(sdm$bsMat2[ell,],p)
    beta[,ell] <- Phi %*% matrix(theta,ncol=1)
  }

  all.beta <- list()
  if( (trace == -1) && (method %in% c('ADAM','GRAD')) ) {
    for(i in c(1:length(fit$all.time))) {
      all.beta[[i]] <- matrix(0,nrow=p,ncol=ntau)
      for(ell in c(1:ntau)) {
        Phi <- create_Phi_matrix(sdm$bsMat2[ell,],p)
        all.beta[[i]][,ell] <- Phi %*% matrix(fit$all.par[,i],ncol=1)
      }
    }	  
  }
  
  beta0<-NULL
  if(return.beta0) {
    # map to interpolated regression solution from theta0
    beta0 <- matrix(0,nrow=p,ncol=ntau)
    for(ell in c(1:ntau)) {
      Phi <- create_Phi_matrix(sdm$bsMat2[ell,],p)
      beta0[,ell] <- Phi %*% matrix(theta0,ncol=1)
    } 
  }
  ##return(list(beta=beta,beta0=beta0,sel.tau0=sel.tau0,all.beta=all.beta,theta=theta,sdm=sdm,spars=spars,fit=fit))
  return(list(beta=beta,all.beta=all.beta,spars=spars,fit=fit))
}




#'  Cubic Spline Quantile Regression with L2-Norm Roughness Penalty (SQR3 or Cubic SQR)
#'
#' This function computes spline quantile regression with cubic splines and L2-norm roughness penalty
#' (SQR3 or cubic SQR) from the response vector and the design matrix on a given set of quantile levels.
#' It uses \code{solve_osqp()}  in the "osqp" package or \code{solve_piqp()} in the "piqp" package.
#' Both are general-purpose quadratic program (QP) solvers in the sparse-matrix form.
#' @param X design matrix (requirement: \code{nrow(X) = length(y)})
#' @param y response vector
#' @param tau  sequence of quantile levels for evaluation
#' @param tau0 sequence of quantile levels for fitting (\code{min(tau0)} <= \code{tau} <= \code{max(tau0)}; 
#' default = \code{tau})
#' @param spar smoothing parameter (default = 1)
#' @param w  weight sequence in penalty (default = \code{rep(1,length(tau0))})
#' @param mthreads if \code{FALSE} (default), set \code{RhpcBLASctl::blas_set_num_threads(1)}
#' @param ztol zero-tolerance parameter to determine the model complexity (default = \code{1e-04})
#' @param type type of QP formulation: \code{'dual'} (default) or \code{'primal'}
#' @param solver QP solver: \code{'piqp'} (default) or \code{'osqp'}
#' @param npar experimental parameter (default = 1)
#' @param all.knots \code{TRUE} or \code{FALSE} (default) as in \code{stats::smooth.spline()}
#' @param control list of control parameters for the QP solver (default = \code{list()})
#' @return A list with the following elements:
#'   \item{coefficients}{matrix of regression coefficients}
#'   \item{derivatives}{matrix of derivatives of regression coefficients}
#'   \item{crit}{sequence critera for smoothing parameter select: (AIC,BIC,GIC)}
#'   \item{np}{sequence of complexity measure as the number of effective parameters}
#'   \item{fid}{sequence of fidelity measure as the quasi-likelihood}
#'   \item{info}{convergence status}
#'   \item{nit}{number of iterations}
#'   \item{K}{number of spline basis functions}
#' @importFrom Matrix Matrix bdiag Diagonal crossprod
#' @importFrom methods as
#' @import 'stats'
#' @import 'splines'
#' @import 'RhpcBLASctl'
#' @import 'quantreg'
#' @import 'piqp'
#' @import 'osqp'
#' @export
sqr3.fit <- function(X,y,tau,tau0=tau,spar=1,w=rep(1,length(tau0)),mthreads=FALSE,ztol=1e-04,
  type=c("dual","primal"),solver=c("piqp","osqp"),npar=c(1,2),all.knots=FALSE,control=list()) {  

  create_spline_matrix <- function(tau0,tau,all.knots=FALSE) {
  # create cubic spline matrix and its second derivative for SQR3
  # input:    tau0 = quantiles for fitting
  #            tau = quantiles for evaluation
  # output: bsMat  = matrix of basis functions
  #         dbsMat = matrix of second derivatives
  #         dbsMatb= matrix of first derivatives
  #         bsMat2 = maxtrix of basis function for interpolation to tau
  # imports: 'splines'
    knots <- stats::smooth.spline(tau0,tau0,all.knots=all.knots)$fit$knot 
    # rescale tau0 and tau to [0,1] for spline matrix to be consistent with smooth.spline
    bsMat <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)))
    dbsMat <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)),derivs=2)
	dbsMatb <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)),derivs=1)
    bsMat2 <- splines::splineDesign(knots,(tau-min(tau))/(max(tau)-min(tau)))
	dbsMat2 <- splines::splineDesign(knots,(tau-min(tau))/(max(tau)-min(tau)),derivs=1)
    return(list(bsMat=bsMat,dbsMat=dbsMat,dbsMatb=dbsMatb,bsMat2=bsMat2,dbsMat2=dbsMat2,knots=knots))
  }

  create_Phi_matrix <- function(bsVec,p) {
  # create spline basis matrix
  # input: bsVec = K-vector of bs functions (K=number of basis functions)
  #            p = number of parameters
  # output:  Phi = p*pK matrix of basis functions
    K <- length(bsVec)
    Phi <- matrix(0,nrow=p,ncol=p*K)
    for(i in c(1:p)) Phi[i,((i-1)*K+1):(i*K)] <- bsVec
    Phi
  }
  
  rescale_penalty2 <- function(tau0,X,sdm,w) {
  # rescale penalty par by weights and weighted l2 norm of second derivatives sdm$dbsMat
  # input: tau0 = sequence of quantile levels
  #           X = design matrix (n-by-p)
  #         sdm = object from create_spline_matrix()
  #           w = weight sequence, same length as tau0
  # output:  c2 = n * r * w
  # dependencies: create_Phi_matrix()
    L <- length(tau0)
    p <- ncol(X)
    n <- nrow(X)
    r <- 0
    r2 <- 0
    for(ell in c(1:L)) {
      # design matrix for spline coefficients
      X0 <- create_Phi_matrix(sdm$bsMat[ell,],p)
      X0 <- X %*% X0
      r <- r + sum(abs(X0))
      X0 <- create_Phi_matrix(sdm$dbsMat[ell,],p)
      r2 <- r2 + w[ell] * sum(diag(crossprod(X0)))
    }
    r <- r / r2
    c2 <- w * r
    c2
  }

  
  Rho <- function(u, tau) u * (tau - (u < 0))
  
  # turn-off blas multithread to run in parallel
  if(!mthreads) {
    ##sink(tempfile())
    RhpcBLASctl::blas_set_num_threads(1)
    ##sink()
  }
  
  # set default values
  if(is.null(ztol)) {
    ztol <- 1e-04
  } else {
    if(ztol <= 0) ztol <- 1e-04
  }
  if(is.null(solver)) {
    solver <- "piqp"
  } else {
    if( !(solver %in% c("piqp","osqp")) ) solver <- "piqp"
  } 
  n <- length(y)
  ntau <-length(tau)
  p <- ncol(X)
  L <- length(tau0)
  # create spline design matrix with knots provided by smooth.spline
  sdm <- create_spline_matrix(tau0,tau,all.knots)
  K <- ncol(sdm$bsMat)
  # rescale penalty parameter
  c2 <- rescale_penalty2(tau0,X,sdm,w)
  c0 <- c2*1000.0**(spar-1.0)
  
  if(solver[1]=="osqp") {
    con <- osqp::osqpSettings(verbose=FALSE,max_iter=50000L,eps_abs=1e-06,eps_rel=1e-06)
  } else {
    con <- piqp::piqp_settings(verbose=FALSE,max_iter=500L,eps_abs=1e-04,eps_rel=1e-05)
  }
  con[(namc <- names(control))] <- control
  accept <- c("solved","solved inaccurate","maximum iterations reached","max iterations reached")
  
  if(type[1]=="primal") {
  
  if(solver[1]=="osqp") {
  
    # set up the primal QP for osqp
    ptm <- proc.time()
    ##B <- Matrix::Matrix(0,p*K,p*K,sparse=TRUE)
	B <- Matrix(0,p*K,p*K,sparse=TRUE)
	A <- NULL
	a <- NULL
    for(ell in c(1:L)) {
      X0 <- as(create_Phi_matrix(sdm$dbsMat[ell,],p),"sparseMatrix")
	  B <- B + 2 * c0[ell] * crossprod(X0)
      X0 <- as(create_Phi_matrix(sdm$bsMat[ell,],p),"sparseMatrix")
      X0 <- as(X %*% X0,"sparseMatrix")
	  A <- rbind(A,X0)
	  a <- c(a,rep(tau0[ell],n))
    }
    rm(X0)	
	A <- cbind(A,Diagonal(n*L,1),-Diagonal(n*L,1))
    A <- rbind(A,cbind(Matrix(0,2*n*L,p*K,sparse=TRUE),Diagonal(2*n*L,1)))
    B <- bdiag(B,Matrix(0,2*n*L,2*n*L,sparse=TRUE))
    a <- c(rep(0,p*K),a,1-a)
    lb <- c(rep(y,L),rep(0,2*n*L)) 
    ub <- c(rep(y,L),rep(Inf,2*n*L))
	setup.time <- (proc.time()-ptm)[3]
	ptm <- proc.time()
    # compute the primal QP solution by osqp
    ##sink(tempfile())
    fit <- osqp::solve_osqp(B, a, A, lb, ub, con)
    ##sink()
    rm(B,a,A,lb,ub)
    theta <- fit$x[c(1:(p*K))]
	solve.time <- (proc.time()-ptm)[3]	
  } else {
    # set up the primal QP for piqp
	ptm <- proc.time()
	B <- Matrix(0,p*K,p*K,sparse=TRUE)
	A <- NULL
	a <- NULL
    for(ell in c(1:L)) {
	  X0 <- as(create_Phi_matrix(sdm$dbsMat[ell,],p),"sparseMatrix")
      B <- B + 2 * c0[ell] * crossprod(X0)
      X0 <- as(create_Phi_matrix(sdm$bsMat[ell,],p),"sparseMatrix")
      X0 <- as(X %*% X0,"sparseMatrix")
	  A <- rbind(A,X0)
      a <- c(a,rep(tau0[ell],n))
    }
    rm(X0) 
	A <- cbind(A,Diagonal(n*L,1),-Diagonal(n*L,1))
    B <- bdiag(B,Matrix(0,2*n*L,2*n*L,sparse=TRUE))
    a <- c(rep(0,p*K),a,1-a)
    b <- rep(y,L)
    lb <- c(rep(-Inf,p*K),rep(0,2*n*L)) 
    ub <- c(rep(Inf,p*K),rep(Inf,2*n*L))
	setup.time <- (proc.time()-ptm)[3]
	ptm <- proc.time()	
    # compute the primal QP solution by piqp
    ##sink(tempfile())
    fit <- piqp::solve_piqp(P=B, c=a, A=A, b=b, x_l=lb, x_u=ub, settings=con)
    ##sink()
    rm(B,a,A,b,lb,ub)
    theta <- fit$x[c(1:(p*K))]
	solve.time <- (proc.time()-ptm)[3]
  }
  
  } else {  # dual
  
  if(solver[1]=="osqp") {
    # set up the dual QP for osqp
    ptm <- proc.time()
	B <- Matrix(0, p*K,p*K,sparse=TRUE)
	A <- NULL
    Xp <- t(X)
    a <- rep(0,p*K)
    for(ell in c(1:L)) {
	  X0 <- as(create_Phi_matrix(sdm$dbsMat[ell,],p),"sparseMatrix")
      B <- B + 2 * c0[ell] * crossprod(X0)
      X0 <- as(create_Phi_matrix(sdm$bsMat[ell,],p),"sparseMatrix")
      X0 <- crossprod(X0,Xp)
      A <- cbind(A,X0)
      a <- a + (1-tau0[ell]) * X0 %*% rep(1,n)
    }
    rm(X0,Xp)
    A <- cbind(-B,A)
	A <- rbind(A,cbind(Matrix(0,n*L,p*K,sparse=TRUE),Diagonal(n*L,1)))
    B <- bdiag(B,Matrix(0,n*L,n*L,sparse=TRUE))
    lb <- c(as.vector(a),rep(0,n*L)) 
    ub <- c(as.vector(a),rep(1,n*L))
    b <- c(rep(0,p*K),-rep(y,L))
	setup.time <- (proc.time()-ptm)[3]
    ptm <- proc.time()
    # compute the dual QP solution by osqp
    ##sink(tempfile())
    fit <- osqp::solve_osqp(B, b, A, lb, ub, con)
    ##sink()
    rm(B,A,b,lb,ub)
    # take the primal QP solution (dual of the dual QP)
    theta <- fit$y[c(1:(p*K))] 
	solve.time <- (proc.time()-ptm)[3]
  } else { 
    # set up the dual QP for piqp
    ptm <- proc.time()
	B <- Matrix(0, p*K,p*K,sparse=TRUE)
	A <- NULL
    Xp <- t(X)
    a <- rep(0,p*K)
    for(ell in c(1:L)) {
	  X0 <- as(create_Phi_matrix(sdm$dbsMat[ell,],p),"sparseMatrix")
      B <- B + 2 * c0[ell] * crossprod(X0)
      X0 <- as(create_Phi_matrix(sdm$bsMat[ell,],p),"sparseMatrix")
      X0 <- crossprod(X0,Xp)
      A <- cbind(A,X0)
      a <- a + (1-tau0[ell]) * X0 %*% rep(1,n)
    }
    rm(X0,Xp)
    a <- as.vector(a)
    A <- cbind(-B,A)
	B <- bdiag(B,Matrix::Matrix(0,n*L,n*L,sparse=TRUE))
    lb <- c(rep(-Inf,p*K),rep(0,n*L)) 
    ub <- c(rep(Inf,p*K),rep(1,n*L))
    b <- c(rep(0,p*K),-rep(y,L))
	setup.time <- (proc.time()-ptm)[3]
    ptm <- proc.time()
    # compute the dual QP solution by piqp
    ##sink(tempfile())
    fit <- piqp::solve_piqp(P=B, c=b, A=A, b=a, x_l=lb, x_u=ub, settings=con)
    ##sink()
    rm(B,a,A,b,lb,ub)
    # take the primal QP solution (dual of the dual QP)
    theta <- fit$y[c(1:(p*K))]
	solve.time <- (proc.time()-ptm)[3]	
  } 
  
  }
  info <- fit$info$status
  nit <- fit$info$iter
  if(!(info %in% accept)) theta <- rep(0,p*K)
  # clean up
  rm(fit)
  solve.time <- (proc.time()-ptm)[3]
  # map to interpolated regression solution
  theta <- matrix(theta,ncol=1)
  beta <- matrix(0,nrow=p,ncol=ntau)
  dbeta <- matrix(0,nrow=p,ncol=ntau)
  for(ell in c(1:ntau)) {
    Phi <- create_Phi_matrix(sdm$bsMat2[ell,],p)
    beta[,ell] <- Phi %*% theta
	Phi <- create_Phi_matrix(sdm$dbsMat2[ell,],p)
    dbeta[,ell] <- Phi %*% theta
  }
  # np = number of points interpolated by the fitted values: invalide for QP
  # bic: aka sic (Koenker 2005, p. 234)
  fid <- rep(0,L)
  np <- rep(0,L)
  beta0 <- matrix(0,nrow=p,ncol=L)  # fitted values
  beta1 <- matrix(0,nrow=p,ncol=L)  # 1st derivatives
  beta2 <- matrix(0,nrow=p,ncol=L)  # 2nd derivatives
  for(ell in c(1:L)) {
    Phi <- create_Phi_matrix(sdm$bsMat[ell,],p)
    beta0[,ell] <- Phi %*% theta
    resid <- y - X %*% beta0[,ell]
    fid[ell] <- mean(Rho(resid,tau0[ell]))
    np[ell] <- sum( abs(resid) < ztol * sd(y) )
	Phi <- create_Phi_matrix(sdm$dbsMatb[ell,],p)
    beta1[,ell] <- Phi %*% theta
    Phi <- create_Phi_matrix(sdm$dbsMat[ell,],p)
    beta2[,ell] <- Phi %*% theta
  }
  np2 <- rep(0,p)
  del <- mean(abs(diff(tau0)))
  for(j in c(1:p)) {
    np2[j] <- sum( abs(beta2[j,]) > ztol ) / p
  }
  if(npar[1]==2) np <- np2
  hqc <- 2 * log( mean(fid) ) + log(log(n)) * mean(np) /n
  bic <- 2 * log( mean(fid) ) + log(n) * mean(np) /n
  aic <- 2 * log( mean(fid) ) + 2 * mean(np) /n
  gic <- aic*0.4 + bic*0.6
  crit <- c(aic,bic,gic)
  names(crit) <- c("AIC","BIC","GIC")
  if(!mthreads) {
    ##sink(tempfile())
    RhpcBLASctl::blas_set_num_threads(Sys.getenv('OPENBLAS_NUM_THREADS'))
    ##sink()
  }
  return(list(coefficients=beta,crit=crit,np=np,fid=fid,info=info,derivatives=dbeta,npar=npar,nit=nit,control=control,
              beta0=beta0,beta1=beta1,beta2=beta2,tau0=tau0,K=K,setup.time=setup.time,solve.time=solve.time))
}




#' Linear Spline Quantile Regression with Total-Variation Roughness Penalty (SQR1 or Linear SQR)
#'
#' This function computes spline quantile regression with linear splines and total-variation roughness penalty 
#' (SQR1 or linear SQR) from the response vector and the design matrix on a given set of quantile levels.
#' It uses the FORTRAN code \code{rqfnb.f} in the "quantreg" package with the kind permission of Dr. R. Koenker
#' or the R  function \code{rq.fit.sfn()} in the same package as a sparse-matrix alternative. Both solve the SQR1 problem as a linear program (LP).
#' @param X design matrix (\code{nrow(X) = length(y)})
#' @param y response vector
#' @param tau  sequence of quantile levels for evaluation
#' @param tau0 sequence of quantile levels for fitting (\code{min(tau0)} <= \code{tau} <= \code{max(tau0)}; 
#' default = \code{tau})
#' @param spar smoothing parameter (default = 1)
#' @param w weight sequence in penalty (default = \code{rep(1,length(tau0)-1)})
#' @param mthreads if \code{FALSE} (default), set \code{RhpcBLASctl::blas_set_num_threads(1)}
#' @param ztol zero-tolerance parameter to determine the model complexity (default = \code{1e-05})
#' @param solver LP solver: \code{'fnb'} (defaut) or \code{'sfn'}
#' @param npar experimental parameter (default = 1)
#' @param all.knots \code{TRUE} or \code{FALSE} (default), same as in \code{stats::smooth.spline()}
#' @return A list with the following elements:
#'   \item{coefficients}{matrix of regression coefficients}
#'   \item{derivatives}{matrix of derivatives of regression coefficients}
#'   \item{crit}{sequence critera for smoothing parameter select: (AIC,BIC,GIC)}
#'   \item{np}{sequence of complexity measure as the number of effective parameters}
#'   \item{fid}{sequence of fidelity measure as the quasi-likelihood}
#'   \item{nit}{number of iterations}
#'   \item{info}{convergence status}
#'   \item{K}{number of spline basis functions}
#' @import 'stats'
#' @import 'splines'
#' @import 'RhpcBLASctl'
#' @import 'quantreg'
#' @importFrom SparseM as.matrix.csr
#' @export
sqr1.fit <- function(X,y,tau,tau0=tau,spar=1,w=rep(1,length(tau0)-1),mthreads=FALSE,ztol=1e-05,
             solver=c("fnb","sfn"),npar=c(1,2),all.knots=FALSE)
{
	
  create_linear_spline_matrix <- function(tau0,tau,all.knots=FALSE) {
  # create linear B-spline matrix and its derivative for SQR1
  # input:       tau0 = quantiles for fitting
  #               tau = quantiles for interpolation
  #         all.knots = TRUE or FALSE, as in smooth.spline()
  # output: bsMat  = matrix of basis functions
  #         dbsMat = matrix of first derivatives
  #         bsMat2 = maxtrix of basis function for interpolation to tau
  # imports: 'splines'
    knots <- stats::smooth.spline(tau0,tau0,all.knots=all.knots)$fit$knot
	knots <- knots[3:(length(knots)-2)]
    # rescale tau0 and tau to [0,1] for spline matrix to be consistent with smooth.spline
    bsMat <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)),ord=2)
    dbsMat <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)),ord=2,derivs=1)
    bsMat2 <- splines::splineDesign(knots,(tau-min(tau))/(max(tau)-min(tau)),ord=2)
	dbsMat2 <- splines::splineDesign(knots,(tau-min(tau))/(max(tau)-min(tau)),ord=2,derivs=1)
    return(list(bsMat=bsMat,dbsMat=dbsMat,bsMat2=bsMat2,dbsMat2=dbsMat2,knots=knots))
  }

  create_Phi_matrix <- function(bsVec,p) {
  # create spline basis matrix for LP and dual LP
  # input: bsVec = K-vector of bs functions (K=number of basis functions)
  #            p = number of parameters
  # output:  Phi = p*pK matrix of basis functions
    K <- length(bsVec)
    Phi <- matrix(0,nrow=p,ncol=p*K)
    for(i in c(1:p)) Phi[i,((i-1)*K+1):(i*K)] <- bsVec
    Phi
  }
  
  rescale_penalty1 <- function(tau0,X,sdm,w) {
  # rescale penalty parameter for linear SQR (SQR1)
  # input: tau0 = L sequence of quantile levels
  #           X = design matrix (n-by-p)
  #         sdm = object from create_linear_spline_matrix()
  #           w = L-1 sequence of weights for penalty
  # output:  c2 = n * c * w
  # dependencies: create_Phi_matrix()
    L <- length(tau0)
	p <- ncol(X)
	n <- nrow(X)
	r <- 0
	for(ell in c(1:L)) {
      # design matrix for spline coefficients
      X0 <- create_Phi_matrix(sdm$bsMat[ell,],p)
      X0 <- X %*% X0
      r <- r + sum(abs(X0))
    }
	r <- r / sum(abs(w * apply(sdm$dbsMat,2,diff)))
    c2 <- w * r
    c2
  }

  Rho <- function(u, tau) u * (tau - (u < 0))
  
  # turn-off blas multithread to run in parallel
  if(!mthreads) {
    ##sink(tempfile())
    RhpcBLASctl::blas_set_num_threads(1)
    ##sink()
  }
  
  # set default values
  if(is.null(ztol)) {
    ztol <- 1e-05
  } else {
    if(ztol <= 0) ztol <- 1e-05
  }
  if(is.null(solver)) {
    solver <- "fnb"
  } else {
    if( !(solver %in% c("fnb","sfn")) ) solver <- "fnb"
  }
  
  n <- length(y)
  ntau <-length(tau)
  p <- ncol(X)
  L <- length(tau0)
  w <- w[1:(L-1)]
  # create spline design matrix with knots provided by smooth.spline
  sdm <- create_linear_spline_matrix(tau0,tau,all.knots)
  K <- ncol(sdm$bsMat)
  # rescale penalty parameter
  c2 <- rescale_penalty1(tau0,X,sdm,w)
  c0 <- c2 * 1000.0**(spar - 1.0)
  ptm <- proc.time()
  # set up the observation vector in dual LP 
  y2 <- c(rep(y,L),rep(0,p*(L-1)))
  # set up the design matrix and rhs vector in dual LP
  X2 <- matrix(0, nrow=n*L+p*(L-1),ncol=p*K)
  rhs <- rep(0,p*K)
  for(ell in c(1:L)) {
    alpha <- tau0[ell]
    cc <- c0[ell]
    # design matrix for spline coefficients
    X0 <- create_Phi_matrix(sdm$bsMat[ell,],p)
    X0 <- X %*% X0
    sel <- ((ell-1)*n+1):(ell*n)
    X2[sel,] <- X0
	rhs <- rhs + (1-alpha) * apply(X0,2,sum) 
	if(ell < L) {
      # design matrix (already weighted) for penalty function 
      Z0 <- create_Phi_matrix(sdm$dbsMat[ell+1,]-sdm$dbsMat[ell,],p)
      # constraint matrix
      sel <- (((ell-1)*p+1):(ell*p))+n*L
      X2[sel,] <- 2 * cc * Z0
	  rhs <- rhs + cc * apply(Z0,2,sum)
	}
  }
  rm(X0,Z0,sel)
  if(solver[1]=='sfn') X2 <- as.matrix.csr(X2)
  setup.time <- (proc.time()-ptm)[3]
  ptm <- proc.time()
  if(solver[1]=='sfn') {
    fit <- quantreg::rq.fit.sfn(a=X2,y=y2,rhs=rhs)
    theta <- fit$coefficients[c(1:(p*K))]
    nit <- fit$it
	info <- "solved"
	##if(fit$ierr !=0) info <- "failed"
	if(!(fit$ierr %in% c(0,17))) info <- "failed"
  } else {	
    # default initial value of n*L+p*(L-1) dual variables
    zeta0 <- c( rep(1-tau0,each=n), rep(0.5, p*(L-1)) )
	# compute the primal-dual LP solution
    fit <- rq.fit.fnb2(x=X2,y=y2,zeta0=zeta0,rhs=rhs)
	zeta <- fit$dualvars
    theta <- fit$coefficients[c(1:(p*K))]
    nit <- fit$nit
	info <- "solved"
	if(fit$info != 0) info <- "failed"
	rm(zeta0)
  }
  if(info != "solved") theta <- rep(0,p*K)
  # clean up
  rm(X2,y2,rhs,fit) 
  solve.time <- (proc.time()-ptm)[3]
  # map to interpolated regression solution
  theta <- matrix(theta,ncol=1)
  beta <- matrix(0,nrow=p,ncol=ntau)
  dbeta <- matrix(0,nrow=p,ncol=ntau)
  for(ell in c(1:ntau)) {
    Phi <- create_Phi_matrix(sdm$bsMat2[ell,],p)
    beta[,ell] <- Phi %*% theta
	Phi <- create_Phi_matrix(sdm$dbsMat2[ell,],p)
    dbeta[,ell] <- Phi %*% theta
  }
  # np = number of points interpolated by the fitted values
  # bic: aka sic (Koenker 2005, p. 234)
  fid <- rep(0,L)
  np <- rep(0,L)
  beta0 <- matrix(0,nrow=p,ncol=L)  # fitted values
  beta1 <- matrix(0,nrow=p,ncol=L)  # 1st derivatives
  for(ell in c(1:L)) {
    Phi <- create_Phi_matrix(sdm$bsMat[ell,],p)
	beta0[,ell] <- Phi %*% theta    
	resid <- y - X %*% beta0[,ell]
	fid[ell] <- mean(Rho(resid,tau0[ell]))
	np[ell] <- sum( abs(resid) < ztol )
	Phi <- create_Phi_matrix(sdm$dbsMat[ell,],p)
	beta1[,ell] <- Phi %*% theta
  }
  np2 <- rep(0,p)
  del <- mean(abs(diff(tau0)))
  for(j in c(1:p)) {
    np2[j] <- sum( abs(beta1[j,])*del^2 > ztol ) / p
  }
  if(npar[1]==2) np <- np2
  hqc <- 2 * log( mean(fid) ) + log(log(n)) * mean(np) /n
  bic <- 2 * log( mean(fid) ) + log(n) * mean(np) /n
  aic <- 2 * log( mean(fid) ) + 2 * mean(np) /n
  gic <- aic*0.4 + bic*0.6
  crit <- c(aic,bic,gic)
  names(crit) <- c("AIC","BIC","GIC")
  if(!mthreads) {
    ##sink(tempfile())
    RhpcBLASctl::blas_set_num_threads(Sys.getenv('OPENBLAS_NUM_THREADS'))
    ##sink()
  }
  return(list(coefficients=beta,crit=crit,np=np,fid=fid,nit=nit,info=info,
              derivatives=dbeta,npar=npar[1],beta1=beta1,beta0=beta0,tau0=tau0,K=K,
			  setup.time=setup.time,solve.time=solve.time))
}



#' Cubic Spline Quantile Regression with L1-Norm Roughness Penalty (SQR)
#'
#' This function computes spline quantile regression solution with cubic splines and L1-norm roughness penalty (SQR)
#' from the response vector and the design matrix on a given set of quantile levels.
#' It uses the FORTRAN code \code{rqfnb.f} in the "quantreg" package with the kind permission of Dr. R. Koenker
#' or the R function \code{rq.fit.sfn()} in the same package as a sparse-matrix alternative. Both solve
#' the SQR problem as a linear program (LP).
#' @param X design matrix (requirement: \code{nrow(X) = length(y)})
#' @param y response vector
#' @param tau  sequence of quantile levels for evaluation
#' @param tau0 sequence of quantile levels for fitting (\code{min(tau0)} <= \code{tau} <= \code{max(tau0)}; 
#' default = \code{tau})
#' @param spar smoothing parameter (default = 1)
#' @param w weight sequence in penalty (default = \code{rep(1,length(tau0))})
#' @param mthreads if \code{FALSE}, set \code{RhpcBLASctl::blas_set_num_threads(1)} (default = \code{TRUE})
#' @param ztol zero-tolerance parameter to determine the model complexity (default = \code{1e-05})
#' @param solver LP solver: \code{'fnb'} (defaut) or \code{'sfn'}
#' @param all.knots \code{TRUE} or \code{FALSE} (default), same as in \code{smooth.spline()}
#' @return A list with the following elements:
#'   \item{coefficients}{matrix of regression coefficients}
#'   \item{derivatives}{matrix of derivatives of regression coefficients}
#'   \item{crit}{criteria values for spar selection: (AIC,BIC,GIC)}
#'   \item{np}{sequence of complexity measure as the number of effective parameters}
#'   \item{fid}{sequence of fidelity measure as the quasi-likelihood}
#'   \item{info}{convergence status}
#'   \item{nit}{number of iterations}
#'   \item{K}{number of spline basis functions}
#' @import 'stats'
#' @import 'splines'
#' @import 'RhpcBLASctl'
#' @export
sqr.fit <- function(X,y,tau,tau0=tau,spar=1,w=rep(1,length(tau0)),
                    mthreads=TRUE,ztol=1e-05,solver=c("fnb","sfn"),all.knots=FALSE) {  

  create_spline_matrix <- function(tau0,tau,all.knots=FALSE) {
  # create cubic B-spline matrix and its second derivative
  # input:    tau0 = quantiles for fitting
  #            tau = quantiles for interpolation
  # output: bsMat  = matrix of basis functions
  #         dbsMat = matrix of second derivatives
  #         bsMat2 = maxtrix of basis function for interpolation to tau
  # imports: 'splines'
    knots <- stats::smooth.spline(tau0,tau0,all.knots=all.knots)$fit$knot 
    # rescale tau0 and tau to [0,1] for spline matrix to be consistent with smooth.spline
    bsMat <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)))
    dbsMat <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)),derivs=2)
	dbsMatb <- splines::splineDesign(knots,(tau0-min(tau0))/(max(tau0)-min(tau0)),derivs=1)
    bsMat2 <- splines::splineDesign(knots,(tau-min(tau))/(max(tau)-min(tau)))
	dbsMat2 <- splines::splineDesign(knots,(tau-min(tau))/(max(tau)-min(tau)),derivs=1)
    return(list(bsMat=bsMat,dbsMat=dbsMat,dbsMatb=dbsMatb,bsMat2=bsMat2,dbsMat2=dbsMat2,knots=knots))
  }

  create_Phi_matrix <- function(bsVec,p) {
  # create spline basis matrix for LP and dual LP for sqr
  # input: bsVec = K-vector of bs functions (K=number of basis functions)
  #            p = number of parameters
  # output:  Phi = p*pK matrix of basis functions
    K <- length(bsVec)
    Phi <- matrix(0,nrow=p,ncol=p*K)
    for(i in c(1:p)) Phi[i,((i-1)*K+1):(i*K)] <- bsVec
    Phi
  }
  
  rescale_penalty <- function(tau0,X,sdm,w) {
  # rescale penalty par by weights and weighted l1 norm of dbsMat for sqr
  # input: tau0 = L sequence of quantile levels
  #           X = design matrix
  #         sdm = object from create_spline_matrix()
  #           w = weights
  # output:  c2 = n * r * w
  # dependencies: create_Phi_matrix()
    L <- length(tau0)
	p <- ncol(X)
	n <- nrow(X)
	r <- 0
	for(ell in c(1:L)) {
      # design matrix for spline coefficients
      X0 <- create_Phi_matrix(sdm$bsMat[ell,],p)
      X0 <- X %*% X0
      r <- r + sum(abs(X0))
    }
	r <- r / sum(abs(w * sdm$dbsMat))
    c2 <- w * r
    c2
  }

  Rho <- function(u, tau) u * (tau - (u < 0))
  
  # turn-off blas multithread to run in parallel
  if(!mthreads) {
    ##sink(tempfile())
    RhpcBLASctl::blas_set_num_threads(1)
    ##sink()
  }
  
  # set default values
  if(is.null(ztol)) {
    ztol <- 1e-05
  } else {
    if(ztol <= 0) ztol <- 1e-05
  }
  if(is.null(solver)) {
    solver <- "fnb"
  } else {
    if( !(solver %in% c("fnb","sfn")) ) solver <- "fnb"
  }

  n <- length(y)
  ntau <-length(tau)
  p <- ncol(X)
  L <- length(tau0)
  # create spline design matrix with knots provided by smooth.spline
  sdm <- create_spline_matrix(tau0,tau,all.knots)
  K <- ncol(sdm$bsMat)
  # rescale penalty parameter
  c2 <- rescale_penalty(tau0,X,sdm,w)
  c0 <- c2 * 1000.0**(spar - 1.0)
  ptm <- proc.time()
  # set up the observation vector in dual LP
  y2 <- c(rep(y,L),rep(0,p*L))
  # set up the design matrix and rhs vector in dual LP
  X2 <- matrix(0, nrow=n*L+p*L,ncol=p*K)
  rhs <- rep(0,p*K)
  for(ell in c(1:L)) {
    alpha <- tau0[ell]
    cc <- c0[ell]
    # design matrix for spline coefficients
    X0 <- create_Phi_matrix(sdm$bsMat[ell,],p)
    X0 <- X %*% X0
    sel <- ((ell-1)*n+1):(ell*n)
    X2[sel,] <- X0
    # design matrix (already weighted) for penalty function 
    Z0 <- create_Phi_matrix(sdm$dbsMat[ell,],p)
    # constraint matrix
    sel <- (((ell-1)*p+1):(ell*p))+n*L
    X2[sel,] <- 2 * cc * Z0
    rhs <- rhs + (1-alpha) * apply(X0,2,sum) + cc * apply(Z0,2,sum)
  }
  rm(X0,Z0,sel)
  if(solver[1]=='sfn') X2 <- as.matrix.csr(X2)
  setup.time <- (proc.time()-ptm)[3]
  ptm <- proc.time()
  if(solver[1]=='sfn') {
    fit <- quantreg::rq.fit.sfn(a=X2,y=y2,rhs=rhs)
    theta <- fit$coefficients[c(1:(p*K))]
    nit <- fit$it
	info <- "solved"
	if(fit$ierr !=0) info <- "failed"
  } else {	
    # default initial value of n*L+p*L dual variables
    zeta0 <- c( rep(1-tau0,each=n), rep(0.5, p*L) )
    # compute the primal-dual LP solution
    fit <- rq.fit.fnb2(x=X2, y=y2, zeta0=zeta0, rhs = rhs)
    theta <- fit$coefficients[c(1:(p*K))]
    zeta <- fit$dualvars
    nit <- fit$nit
	info <- "solved"
	if(fit$info != 0) info <- "failed"
	rm(zeta0)
  }
  # clean up
  rm(X2,y2,rhs,fit)
  solve.time <- (proc.time()-ptm)[3]
  # map to interpolated regression solution
  theta <- matrix(theta,ncol=1)
  beta <- matrix(0,nrow=p,ncol=ntau)
  dbeta <- matrix(0,nrow=p,ncol=ntau)
  for(ell in c(1:ntau)) {
    Phi <- create_Phi_matrix(sdm$bsMat2[ell,],p)
    beta[,ell] <- Phi %*% theta
	Phi <- create_Phi_matrix(sdm$dbsMat2[ell,],p)
    dbeta[,ell] <- Phi %*% theta
  }
  # np = number of points interpolated by the fitted values
  # bic: aka sic (Koenker 2005, p. 234)
  fid <- rep(0,L)
  np <- rep(0,L)
  beta0 <- matrix(0,nrow=p,ncol=L)  # fitted values
  beta1 <- matrix(0,nrow=p,ncol=L)  # 1st derivatives
  beta2 <- matrix(0,nrow=p,ncol=L)  # 2nd derivatives
  for(ell in c(1:L)) {
    Phi <- create_Phi_matrix(sdm$bsMat[ell,],p)
	beta0[,ell] <- Phi %*% theta    
	resid <- y - X %*% beta0[,ell]
	fid[ell] <- mean(Rho(resid,tau0[ell]))
	np[ell] <- sum( abs(resid) < ztol )
	Phi <- create_Phi_matrix(sdm$dbsMatb[ell,],p)
	beta1[,ell] <- Phi %*% theta
	Phi <- create_Phi_matrix(sdm$dbsMat[ell,],p)
	beta2[,ell] <- Phi %*% theta
  }
  hqc <- 2 * log( mean(fid) ) + log(log(n)) * mean(np) /n
  bic <- 2 * log( mean(fid) ) + log(n) * mean(np) /n
  aic <- 2 * log( mean(fid) ) + 2 * mean(np) /n
  gic <- aic*0.4 + bic*0.6
  crit <- c(aic,bic,gic)
  names(crit) <- c("AIC","BIC","GIC")
  return(list(coefficients=beta,crit=crit,np=np,fid=fid,nit=nit,info=info,derivatives=dbeta,beta0=beta0,
              beta1=beta1,beta2=beta2,tau0=tau0,K=K,setup.time=setup.time,solve.time=solve.time))
}


#' Spline Quantile Regression by Formula
#'
#' This function computes the spline quantile regression solution SQR, SQR1, or SAR3 
#' on a given set of quantile levels from a regression formula with or without user-supplied smoothing parameter.
#' SQR represents the regression coefficients as cubic spline functions of the quantile level and employs the L1-norm of their second derivative as roughness penalty.
#' SQR1 represents the regression coefficients as linear spline functions and employs the total variation of their first derivatives
#' as roughness penalty. SQR3 also represents the regression coefficients as cubic spline functions but employs the L2-norm of their second
#' derivatives as roughness penalty. SQR and SQR1 are solved as linear program (LP) using \code{sqr.fit()} and 
#' \code{sqr1.fit()}, respectively. SQR3 is solved as quadratic program (QP) using \code{sqr3.fit()}.
#' @param formula formula object, with the response on the left of a ~ operator, and the terms, separated by + operators, on the right. 
#' @param tau  sequence of quantile levels for evaluation 
#' @param tau0 sequence of quantile levels for fitting (\code{min(tau0)} <= \code{tau} <= \code{max(tau0)}; 
#' default = \code{tau})
#' @param spar smoothing parameter, selected automatically by \code{criterion} if \code{spar = NULL} (default)
#' @param w weight sequence in penalty (default = \code{rep(1,length(tau0))})
#' @param mthreads if \code{FALSE} (default), set \code{RhpcBLASctl::blas_set_num_threads(1)}
#' @param criterion criterion for smoothing parameter selection (\code{"AIC"}, \code{"BIC"}, or \code{"GIC"})
#' @param method \code{'sqr'} (default), \code{'sqr1'}, or \code{'sqr3'}
#' @param type type of QP formulation for SQR3: \code{'dual'} (default) or \code{'primal'}
#' @param ztol zero-tolerance parameter to determine the model complexity 
#' (default = \code{NULL}: set internally to \code{1e-5} for SQR and SQR1 or \code{1e-4} for SQR3)
#' @param solver \code{'fnb'} or \code{'sfn'} for  SQR and SQR1; \code{'piqp'} or \code{'osqp'} for SQR3
#' (default = \code{NULL}: set internally to \code{'fnb'} for SQR and SQR1 or \code{'piqp'} for SQR3)
#' @param interval interval for \code{spar} optimization (default: \code{c(-1.5,1.5)} for SQR and SQR1 
#' or \code{c(1.0,2.5)} for SQR3)
#' @param npar experimental parameter (default = 1)
#' @param all.knots \code{TRUE} or \code{FALSE} (default), same as in \code{stats::smooth.spline()}
#' @param control list of control parameters for QP solvers \code{'piqp'} and \code{'osqp'}
#' @param data data.frame object containing the observations
#' @param subset an optional vector specifying a subset of observations to be used
#' @param na.action a function to filter missing data (see \code{rq()} in the 'quantreg' package)
#' @param model if \code{TRUE} then the model frame is returned (needed for calling \code{summary} subsequently) 
#' @return object of \code{sqr.fit()}, \code{sqr1.fit()}, or \code{sqr3.fit()}
#' @import 'stats'
#' @import 'splines'
#' @import 'RhpcBLASctl'
#' @import 'quantreg'
#' @import 'osqp'
#' @import 'piqp'
#' @export
#' @examples
#' library(quantreg)
#' data(engel)
#' engel$income <- (engel$income - mean(engel$income))/1000
#' tau <- seq(0.1,0.9,0.05)
#' fit <- rq(foodexp ~ income,tau=tau,data=engel)
#' fit.sqr1 <- sqr(foodexp ~ income,tau=tau,spar=0.5,method="sqr1",data=engel)
#' fit.sqr3 <- sqr(foodexp ~ income,tau=tau,spar=0.5,method="sqr3",data=engel)
#' par(mfrow=c(1,1),pty="m",lab=c(10,10,2),mar=c(4,4,2,1)+0.1,las=1)
#' plot(tau,fit$coef[2,],xlab="Quantile Level",ylab="Coeff1")
#' lines(tau,fit.sqr1$coef[2,])
#' lines(tau,fit.sqr3$coef[2,],col=2)
sqr <- function (formula, tau = seq(0.1,0.9,0.2), tau0 = tau, spar = NULL, w = rep(1,length(tau0)), 
 mthreads = FALSE, criterion = c("AIC","BIC","GIC"), method = c("sqr","sqr1","sqr3"), type = c("dual","primal"), 
 ztol = NULL, solver = NULL, interval=NULL, npar=c(1,2), 
 all.knots=FALSE, control=list(), data, subset, na.action, model = TRUE) 
{	
    contrasts <- NULL 
    call <- match.call()
    mf <- match.call(expand.dots = FALSE)
    m <- match(c("formula", "data", "subset", "na.action"), names(mf), 0)
    mf <- mf[c(1, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval.parent(mf)
    mt <- attr(mf, "terms")
    weights <- as.vector(model.weights(mf))
    tau <- sort(unique(tau))
    eps <- .Machine$double.eps^(2/3)
    if (any(tau == 0)) 
      tau[tau == 0] <- eps
    if (any(tau == 1)) 
      tau[tau == 1] <- 1 - eps
    y <- model.response(mf)
    X <- model.matrix(mt, mf, contrasts)
    vnames <- dimnames(X)[[2]]

    criterion <- criterion[1]
	if(is.null(spar) | length(spar) > 1) {
      if(!(criterion %in% c("AIC","BIC","GIC"))) criterion <- "AIC" 
	}
	
    method <- method[1]
    if(!(method %in% c("sqr","sqr1","sqr3"))) method <- "sqr" 
    
    npar <- npar[1]
    if(!(npar %in% c(1,2))) npar <- 1 
	
	if(method=="sqr3") {
      type <- type[1]
      if(!(type %in% c("dual","primal"))) type <- "dual" 
	}

    if(is.null(interval)) {
      if(method=="sqr3") {
        interval <- c(1.0,2.5)
	  } else {
        interval <- c(-1.5,1.5)
	  }
	}	  
	
    Rho <- function(u, tau) u * (tau - (u < 0))
	
    sqr_obj <- function(spar,X,y,tau,tau0,w,criterion,mthreads,method,type,ztol,solver,npar,all.knots,control) {
		if(method[1]=="sqr3") {
		  crit <- sqr3.fit(X,y,tau,tau0,spar,w=w,mthreads=mthreads,
		                  ztol=ztol,solver=solver,npar=npar,
						  type=type,all.knots=all.knots,control=control)$crit
		} else {
		  if(method[1]=="sqr1") {
            crit <- sqr1.fit(X,y,tau,tau0,spar,w=w,mthreads=mthreads,
		                  ztol=ztol,solver=solver,npar=npar,all.knots=all.knots)$crit
		  } else {
            crit <- sqr.fit(X,y,tau,tau0,spar,w=w,mthreads=mthreads,
		                  ztol=ztol,solver=solver,all.knots=all.knots)$crit		  
		  }
        }
        out <- crit[1]
        if(criterion=="BIC") out <- crit[2]
	    if(criterion=="GIC") out <- crit[3]
		out
    }


    if (length(tau0) > 4) {
        if (any(tau0 <= 0) || any(tau0 >= 1)) 
            stop("invalid tau0:  tau0 should be > 0 and < 1")
	    if( max(tau) > max(tau0) || min(tau) < min(tau0) )
		    stop("invalid tau: tau should be within the range of tau0")

        if(method=="sqr1") {	
          w <- w[1:(length(tau0)-1)]
        } else {
          w <- w[1:length(tau0)]
	    }
	
		if(is.null(spar)) {
          spar0 <- stats::optimize(sqr_obj,interval=interval,
	         X=X,y=y,tau=tau,tau0=tau0,w=w,criterion=criterion,mthreads=mthreads,
			 method=method,type=type,ztol=ztol,solver=solver,npar=npar,all.knots=all.knots,control=control)$minimum
        } else {
          if(length(spar) > 1) {
		    crit <- rep(NA,length(spar))
			for(i in c(1:length(spar))) {
			  crit[i] <- sqr_obj(spar[i],X=X,y=y,tau=tau,tau0=tau0,w=w,criterion=criterion,mthreads=mthreads,
			               method=method,type=type,ztol=ztol,solver=solver,npar=npar,all.knots=all.knots,control=control)
			}
            spar0 <- spar[which.min(crit)]
		  } else {
		    spar0 <- spar
		  }
	    }

        coef <- matrix(0, ncol(X), length(tau))
        rho <- rep(0, length(tau))
        fitted <- resid <- matrix(0, nrow(X), length(tau))
		if(method=="sqr3") {
		  z <- sqr3.fit(X,y,tau=tau,tau0=tau0,spar=spar0,w=w,mthreads=mthreads,ztol=ztol,
		                solver=solver,npar=npar,type=type,all.knots=all.knots,control=control)
		} else {
		  if(method=="sqr1") {
            z <- sqr1.fit(X,y,tau=tau,tau0=tau0,spar=spar0,w=w,mthreads=mthreads,ztol=ztol,solver=solver,npar=npar,all.knots=all.knots)
		  } else {
            z <- sqr.fit(X,y,tau=tau,tau0=tau0,spar=spar0,w=w,mthreads=mthreads,ztol=ztol,solver=solver,all.knots=all.knots)		  
		  }
		}
		
        coef <- z$coefficients
		fitted <- y - X %*% coef
        resid <- y - fitted
        for(i in c(1:length(tau))) rho[i] <- sum(Rho(resid[,i], tau[i]))
        taulabs <- paste("tau=", format(round(tau, 3)))
        dimnames(coef) <- list(vnames, taulabs)
        dimnames(resid)[[2]] <- taulabs
        fit <- z
        fit$coefficients <- coef
        fit$residuals <- resid
		fit$fitted.values <- fitted
		fit$spars <- spar
		fit$spar <- spar0
		fit$criterion <- criterion
        class(fit) <- "rqs"
    }
    else {
	    stop("need at least 5 quantile levels for fitting")
    }
	fit$w <- w
	fit$sqrmethod <- method
	fit$all.knots <- all.knots
    fit$na.action <- attr(mf, "na.action")
    fit$formula <- formula
    fit$terms <- mt
    fit$xlevels <- .getXlevels(mt, mf)
    fit$call <- call
    fit$tau <- tau
	fit$tau0 <- tau0
    fit$weights <- weights
    fit$residuals <- drop(fit$residuals)
    fit$rho <- rho
    fit$method <- "br"
	fit$sqr <- method
	fit$solver <- solver
	fit$ztol <- ztol
	fit$criterion<-criterion
	fit$interval<-interval
	fit$control<-control
    fit$fitted.values <- drop(fit$fitted.values)
    attr(fit, "na.message") <- attr(m, "na.message")
    if (model) 
        fit$model <- mf
    fit
}



#' Trigonometric Spline Quantile Regression (TSQR) of Time Series
#'
#' This function computes trigonometric spline quantile regression (TSQR) for 
#' univariate time series at a single frequency using \code{sqr.fit()}, \code{sqr1.fit()}, or \code{sqr3.fit()}.
#' @param y time series
#' @param f0 frequency in [0,1)
#' @param tau  sequence of quantile levels for evaluation 
#' @param tau0 sequence of quantile levels for fitting (\code{min(tau0)} <= \code{tau} <= \code{max(tau0)}; 
#' default = \code{tau})
#' @param spar smoothing parameter (default = 1)
#' @param w weight sequence in penalty (default = \code{rep(1,length(tau0))})
#' @param mthreads if \code{FALSE} (default), set \code{RhpcBLASctl::blas_set_num_threads(1)}
#' @param prepared if \code{TRUE}, intercept is removed and coef of cosine is doubled when \code{f0 = 0.5}
#' @param method \code{'sqr'} (default), \code{"sqr1"}, or \code{'sqr3'}
#' @param ztol a zero tolerance paramete to determine the model complexity 
#' (default = \code{NULL}: set internally to \code{1e-5} for SQR and SQR1 or \code{1e-4} for SQR3)
#' @param solver \code{'fnb'} or \code{'sfn'} for SQR and SQR1; \code{'piqp'} or \code{'osqp'} for SQR3
#' (default = \code{NULL}: set internally to \code{'fnb'}for SQR and SQR1 or \code{'piqp'} for SQR3)
#' @param all.knots \code{TRUE} or \code{FALSE} (default), as in \code{stats::smooth.spline()}
#' @param control list of control parameters for QP solvers \code{'piqp'} and \code{'osqp'} 
#' @return object of \code{sqr.fit()}, \code{sqr1.fit()}, or \code{sqr3.fit()}
#' @import 'stats'
#' @import 'splines'
#' @import 'RhpcBLASctl'
#' @import 'quantreg'
#' @import 'osqp'
#' @import 'piqp'
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' tau0 <- seq(0.1,0.9,0.2)
#' fit <- tqr.fit(y,f0=0.1,tau=tau)
#' fit.sqr1 <- tsqr.fit(y,f0=0.1,tau=tau,tau0=tau0,spar=0.2,method='sqr1')
#' fit.sqr3 <- tsqr.fit(y,f0=0.1,tau=tau,tau0=tau0,spar=1,method='sqr3')
#' plot(tau,fit$coef[1,],type='p',xlab='QUANTILE LEVEL',ylab='TQR COEF')
#' lines(tau,fit.sqr1$coef[1,])
#' lines(tau,fit.sqr3$coef[1,],col=2)
tsqr.fit <- function(y,f0,tau,tau0=tau,spar=1,w=rep(1,length(tau0)),mthreads=FALSE,prepared=TRUE,
             method=c("sqr","sqr1","sqr3"),ztol=NULL,solver=NULL,all.knots=FALSE,control=list()) {

  tqr_create_design_matrix <- function(f0,n) {
  # create design matrix for trignometric regression of time series
  # input: f0 = frequency in [0,1)
  #         n = length of time series
    tt <- c(1:n)
    if(f0 != 0.5 & f0 != 0) {
      X <- cbind(rep(1,n),cos(2*pi*f0*tt),sin(2*pi*f0*tt))
    }
    if(f0 == 0.5) {
      X <- cbind(rep(1,n),cos(2*pi*0.5*tt))
    }
    if(f0 == 0) {
      X <- matrix(rep(1,n),ncol=1)
    }
    X
  }
  
  tqr_fix_coef <- function(coef) {
  # prepare coef from tqr for qdft
  # input:  coef = p*ntau tqr coefficient matrix from tqr.fit()
  # output: 2*ntau matrix of tqr coefficients
    ntau <- ncol(coef)
    if(nrow(coef)==1) {   
      # for f = 0
      coef <- rbind(rep(0,ntau),rep(0,ntau))
    } else if(nrow(coef)==2) {  
      # for f = 0.5: rescale coef of cosine by 2 so qdft can be defined as usual
      coef <- rbind(2*coef[2,],rep(0,ntau))
    } else {
      # for f in (0,0.5)
      coef <- coef[-1,]
    }
    coef
  }

  if(method=="sqr1") {	
    w <- w[1:(length(tau0)-1)]
  } else {
    w <- w[1:length(tau0)]
  }
  
  n <- length(y)
  X <- tqr_create_design_matrix(f0,n)
  if(method[1]=="sqr3") {
    fit <- sqr3.fit(X,y,tau,tau0=tau0,spar=spar,w=w,mthreads=mthreads,ztol=ztol,solver=solver,all.knots=all.knots,control=control)  
  } else {
    if(method[1]=="sqr1") {
	  fit <- sqr1.fit(X,y,tau,tau0=tau0,spar=spar,w=w,mthreads=mthreads,ztol=ztol,solver=solver,all.knots=all.knots)  
	} else {
      fit <- sqr.fit(X,y,tau,tau0=tau0,spar=spar,w=w,mthreads=mthreads,ztol=ztol,solver=solver,all.knots=all.knots)  
	}
  }
  fit$fitted.values <- matrix(0,ncol=length(tau),nrow=n)
  for(i in c(1:length(tau))) {
    fit$fitted.values[,i] <- X %*% fit$coefficients[,i]
  }
  if(prepared) fit$coefficients <- tqr_fix_coef(fit$coefficients)
  fit
}




#' Spline Quantile Discrete Fourier Transform (SQDFT) of Time Series Given Smoothing Parameter
#'
#' This function computes spline quantile discrete Fourier transform (SQDFT) for univariate or multivariate time series
#' through trigonometric spline quantile regression with user-supplied spar.
#' @param y vector or matrix of time series (if matrix, \code{nrow(y)} = length of time series)
#' @param tau  sequence of quantile levels for evaluation 
#' @param tau0 sequence of quantile levels for fitting (\code{min(tau0)} <= \code{tau} <= \code{max(tau0)}; 
#' default = \code{tau})
#' @param spar smoothing parameter (default = 1)
#' @param w weight sequence in penalty (default = \code{rep(1,length(tau0))})
#' @param method \code{'sqr'} (default), \code{"sqr1"}, or \code{'sqr3'}
#' @param ztol zero-tolerance parameter to determine the model complexity 
#' (default = \code{NULL}: set internally to \code{1e-5} for SQR and SQR1 or \code{1e-4} for SQR3)
#' @param solver \code{'fnb'} or \code{'sfn'} for SQR and SQR1; \code{'piqp'} or \code{'osqp'} for SQR3
#' (default = \code{NULL}: set internally to \code{'fnb'} for SQR and SQR1 or \code{'piqp'} for SQR3)
#' @param all.knots \code{TRUE} or \code{FALSE} (default), as in \code{stats::smooth.spline()}
#' @param control list of control parameters for QP solvers \code{'piqp'} and \code{'osqp'} (default = \code{list()})
#' @param n.cores number of cores for parallel computing (default = 1)
#' @param cl pre-existing cluster for repeated parallel computing (default = \code{NULL})
#' @return A list with the following elements:
#'   \item{coefficients}{matrix of regression coefficients}
#'   \item{qdft}{matrix or array of the spline quantile discrete Fourier transform of \code{y}}
#'   \item{crit}{criteria for smoothing parameter selection: (AIC,BIC,GIC)}
#'   \item{nit}{maximum number of iterations}
#' @import 'stats'
#' @import 'splines'
#' @import 'RhpcBLASctl'
#' @import 'quantreg'
#' @import 'osqp'
#' @import 'piqp'
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' tau0 <- seq(0.1,0.9,0.2)
#' y.sqdft <- sqdft.fit(y,tau,tau0=tau0,spar=0.2,method="sqr1")$qdft
sqdft.fit <- function(y,tau,tau0=tau,spar=1,w=rep(1,length(tau0)),method=c("sqr","sqr1","sqr3"),
             ztol=NULL,solver=NULL,all.knots=FALSE,control=list(),n.cores=1,cl=NULL) {

  z <- function(x) { x <- matrix(x,nrow=2); x[1,]-1i*x[2,] }

  extend_qdft <- function(y,tau,result,sel.f) {
  # define qdft from tqr coefficients
  # input: y = ns*nc-matrix or ns-vector of time series
  #        tau = ntau-vector of quantile levels
  #        result = list of qdft in (0,0.5]
  #        sel.f = index of Fourier frequencies in (0,0.5)
  # output: out = nc*ns*ntau-array or ns*ntau matrix of qdft

    if(!is.matrix(y)) y <- matrix(y,ncol=1)
    nc <- ncol(y)
    ns <- nrow(y)
    ntau <- length(tau)  
    out <- array(NA,dim=c(nc,ns,ntau))
    for(k in c(1:nc)) {
      # define QDFT at freq 0 as ns * quantile
      out[k,1,] <- ns * stats::quantile(y[,k],tau)
      # retrieve qdft for freq in (0,0.5]
      tmp <- matrix(unlist(result[[k]]),ncol=ntau,byrow=TRUE)
      # define remaining values by conjate symmetry (excluding f=0.5) 
      tmp2 <- NULL
      for(j in c(1:ntau)) {
        tmp2 <- cbind(tmp2,rev(Conj(tmp[sel.f,j])))
      }
      # assemble & rescale everything by ns/2 so that periodogram = |dft|^2/ns
      out[k,c(2:ns),] <- rbind(tmp,tmp2) * ns/2
    }
    if(nc == 1) out <- matrix(out[1,,],ncol=ntau)
    out
  }

  mthreads <- FALSE
  if(!is.matrix(y)) y <- matrix(y,ncol=1)
  ns <- nrow(y)
  nc <- ncol(y)
  f2 <- c(0:(ns-1))/ns
  # compute QR at half of the Fourier frequencies, excluding 0
  f <- f2[which(f2 > 0 & f2 <= 0.5)]
  sel.f <- which(f < 0.5)
  nf <- length(f)
  ntau <- length(tau)
  keep.cl <- TRUE
  if(n.cores>1 & is.null(cl)) {
    cl <- parallel::makeCluster(n.cores)
    parallel::clusterExport(cl, c("tsqr.fit","sqr.fit","sqr1.fit","sqr3.fit","rq.fit.fnb2"))
    doParallel::registerDoParallel(cl)
    keep.cl <- FALSE
  }
  `%dopar%` <- foreach::`%dopar%`
  `%do%` <- foreach::`%do%`
  # qdft for f in (0,0.5]
  i <- 0  
  result <- list()
  result2 <- list()
  result3 <- list()
  for(k in c(1:nc)) {
    yy <- y[,k] 
    if(n.cores>1) {
      tmp <- foreach::foreach(i=1:nf, .packages=c('quantreg','osqp','piqp',"Matrix")) %dopar% {
		fit <- tsqr.fit(yy,f[i],tau,tau0=tau0,spar=spar,w=w,
		                mthreads=mthreads,method=method,ztol=ztol,solver=solver,
                        all.knots=all.knots,control=control)
	    coef <- fit$coef
		crit <- fit$crit
		nit <- fit$nit
		list(coef=coef,crit=crit,nit=nit)
      }
    } else {
      tmp <- foreach::foreach(i=1:nf) %do% {
        fit <- tsqr.fit(yy,f[i],tau,tau0=tau0,spar=spar,w=w,
		                mthreads=mthreads,method=method,ztol=ztol,solver=solver,all.knots=all.knots,control=control)
		coef <- fit$coef
		crit <- fit$crit
		nit <- fit$nit
		list(coef=coef,crit=crit,nit=nit)
      }
    }
    # tmp$coef = a list over freq of 2 x ntau coefficiets 
    tmp.coef <- lapply(tmp,FUN=function(x) {apply(x$coef,2,z)})
    result[[k]] <- tmp.coef
	tmp.crit <- lapply(tmp,FUN=function(x) { x$crit } )
	result2[[k]] <- tmp.crit 
	tmp.nit <- lapply(tmp,FUN=function(x) { x$nit } )
	result3[[k]] <- tmp.nit
  }
  if(n.cores>1 & !keep.cl) {
    parallel::stopCluster(cl)
    cl <- NULL
  }
  
  # extend qdft to f=0 and f in (0.5,1) 
  out <- extend_qdft(y,tau,result,sel.f)

  out2 <- rep(0,length(result2[[1]][[1]]))
  for(k in c(1:nc)) {
    for(j in c(1:nf)) out2 <- out2 + result2[[k]][[j]]
  }
  out2 <- out2 / (nc*nf)
  
  out3 <- rep(0,nc)
  for(k in c(1:nc)) {
    out3[k] <- max(unlist(result3[[k]]))
  }
  out3 <- max(out3)
 
  return(list(qdft=out,crit=out2,nit=out3))
}



#' Spline Quantile Discrete Fourier Transform (SQDFT) of Time Series
#'
#' This function computes spline quantile discrete Fourier transform (SQDFT) for univariate or multivariate time series
#' through trigonometric spline quantile regression.
#' @param y vector or matrix of time series (if matrix, \code{nrow(y)} = length of time series)
#' @param tau  sequence of quantile levels for evaluation 
#' @param tau0 sequence of quantile levels for fitting (\code{min(tau0)} <= \code{tau} <= \code{max(tau0)}; 
#' default = \code{tau})
#' @param spar smoothing parameter, selected automatically by \code{criterion} if \code{spar = NULL} or if \code{length(spar) > 1}
#' @param w weight sequence in penalty (default = \code{rep(1,length(tau0))})
#' @param criterion criterion for smoothing parameter selection: \code{"AIC"} (default), \code{"BIC"}, or \code{"GIC"}
#' @param method \code{'sqr'}(default), \code{'sqr1'}, or \code{'sqr3'}
#' @param ztol zero-tolerance parameter to determine the model complexity 
#' (default = \code{NULL}: set internally to \code{1e-5} for SQR and SQR1 or \code{1e-4} for SQR3)
#' @param solver \code{'fnb'} or \code{'sfn'} for SQR and SQR1; \code{'piqp'} or \code{'osqp'} for SQR3
#' (default = \code{NULL}: set internally to \code{'fnb'} for SQR and SQR1 or \code{'piqp'} for SQR3)
#' @param interval interval for \code{spar} optimization (default: \code{c(-1.5,1.5)} for SQR and SQR1 
#' or \code{c(0,2.5)} for SQR3)
#' @param all.knots \code{TRUE} or \code{FALSE} (default), as in \code{stats::smooth.spline()}
#' @param control list of control parameters for QP solvers \code{'piqp'} and \code{'osqp'} (default = \code{list()})
#' @param n.cores number of cores for parallel computing (default = 1)
#' @param cl pre-existing cluster for repeated parallel computing (default = \code{NULL})
#' @return A list with the following elements:
#'   \item{coefficients}{matrix of regression coefficients}
#'   \item{qdft}{matrix or array of the spline quantile discrete Fourier transform of \code{y}}
#'   \item{crit}{criteria for smoothing parameter selection: (AIC,BIC,GIC)}
#'   \item{nit}{maximum number of iterations}
#'   \item{spar}{optimal value of smoothing parameter}
#' @import 'stats'
#' @import 'splines'
#' @import 'RhpcBLASctl'
#' @import 'quantreg'
#' @import 'osqp'
#' @import 'piqp'
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @export
#' @examples
#' y <- stats::arima.sim(list(order=c(1,0,0), ar=0.5), n=64)
#' tau <- seq(0.1,0.9,0.05)
#' y.sqdft <- sqdft(y,tau,spar=0.2,method="sqr1")$qdft
#' plot(y.sqdft[,2])
sqdft <- function(y,tau,tau0=tau,spar=NULL,w=rep(1,length(tau0)),criterion=c("AIC","BIC","GIC"),
             method=c("sqr","sqr1","sqr3"),ztol=NULL,solver=NULL,interval=NULL,
			 all.knots=FALSE,control=list(),n.cores=1,cl=NULL)
{

    criterion <- criterion[1]
	if(is.null(spar) | length(spar) > 1) {
      if(!(criterion %in% c("AIC","BIC","GIC"))) criterion <- "AIC" 
	}
	
    method <- method[1]
    if(!(method %in% c("sqr","sqr1","sqr3"))) method <- "sqr" 
	
    if(is.null(solver)) {
      if(method=="sqr3") {
        solver <- "piqp"
	  } else {
        solver <- "fnb"
	  }
	} else {
      solver <- solver[1]
      if(method=="sqr3") {
	    if( !(solver %in% c("osqp","piqp")) ) solver <- "piqp"
      } else {
	    if( !(solver %in% c("fnb","sfn")) ) solver <- "fnb"
      }
    }

    if(is.null(ztol)) {
      if(method=="sqr3") {
        ztol <- 1e-04
	  } else {
        ztol <- 1e-05
	  }
	} else {
	  ztol <- ztol[1]
	}
	
    if(is.null(interval)) {
      if(method=="sqr3") {
        interval <- c(0.0,2.5)
	  } else {
        interval <- c(-1.5,1.5)
	  }
	}	  
	
    if(method=="sqr1") {	
      w <- w[1:(length(tau0)-1)]
    } else {
      w <- w[1:length(tau0)]
	}
		
    sqdft_obj <- function(spar,y,tau,tau0,w,method,ztol,solver,criterion,all.knots,control,n.cores,cl) {
        crit <- sqdft.fit(y,tau,tau0=tau0,spar=spar,w=w,method=method,
	                  ztol=ztol,solver=solver,all.knots=all.knots,control=control,n.cores=n.cores,cl=cl)$crit
		out <- crit[1]
		if(criterion=="BIC") out <- crit[2]
        if(criterion=="GIC") out <- crit[3]
        out
    }
  
    if(is.null(spar)) {
      spar0 <- stats::optimize(sqdft_obj,interval=interval,y=y,tau=tau,tau0=tau0,w=w,method=method,
			ztol=ztol,solver=solver,criterion=criterion,
            all.knots=all.knots,control=control,n.cores=n.cores,cl=cl)$minimum
    } else {
      if(length(spar)>1) {
  	    crit <- rep(NA,length(spar))
	    for(j in c(1:length(spar))) {
	      crit[j] <- sqdft_obj(spar[j],y,tau,tau0,w,method,ztol,solver,criterion,all.knots,control,n.cores,cl)
		}
		spar0 <- spar[which.min(crit)]
	  } else {
	    spar0 <- spar
      }
    }
    fit <- sqdft.fit(y,tau,tau0=tau0,spar=spar0,w=w,method=method,
	                 ztol=ztol,solver=solver,control=control,all.knots=all.knots,n.cores=n.cores,cl=cl)
    return(list(qdft=fit$qdft,spar=spar0,crit=fit$crit,nit=fit$nit))
}



#' Bootstrap of Spline Quantile Regression (SQR) Coefficients
#'
#' This function generates bootstrap samples for the SQR coefficients and their derivatives
#' @param fit output from \code{sqr()}
#' @param nsim number of bootstrap samples (default = 1000)
#' @param blocklength blocklength for block sampling scheme (default = 1)
#' @param n.cores number of cores for parallel computing (default = 1)
#' @param mthreads if \code{FALSE} (default), set \code{RhpcBLASctl::blas_set_num_threads(1)}
#' @param seed seed of random number generator (default = 1234567)
#' @return A list with the following elements:
#'   \item{coefficients}{array of bootstrap regression coefficients}
#'   \item{derivatives}{array of boostrap derivatives of regression coefficient }
#'   \item{info}{sequence convergence status}
#' @import 'stats'
#' @import 'splines'
#' @import 'RhpcBLASctl'
#' @import 'quantreg'
#' @import 'osqp'
#' @import 'piqp'
#' @import 'boot'
#' @import 'foreach'
#' @import 'parallel'
#' @import 'doParallel'
#' @export
boot.sqr <- function(fit,nsim=1000,blocklength=1,n.cores=1,mthreads=FALSE,seed=1234567) {
# generate boostrap samples for SQR using the xy-pair resampling method
  n <- nrow(fit$model)
  set.seed(seed)
  idx <- list()
  if(blocklength > 1) {
    index <- function(x) { x }
    tmp <-boot::tsboot(c(1:n),statistic=index,l=blocklength,R=nsim,sim="fixed")$t
    for(sim.idx in c(1:nsim)) idx[[sim.idx]] <- tmp[sim.idx,]
  } else {
    for(sim.idx in c(1:nsim)) idx[[sim.idx]] <- sample(c(1:n),n,replace=TRUE)
  }
  `%dopar%` <- foreach::`%dopar%`
  `%do%` <- foreach::`%do%`
  if(n.cores > 1) {
  	if(!mthreads) {
	  ##sink(tempfile())
      RhpcBLASctl::blas_set_num_threads(1)
	  ##sink()
    }
    cl <- parallel::makeCluster(n.cores)
    parallel::clusterExport(cl, c("sqr","sqr.fit","sqr1.fit","sqr3.fit","rq.fit.fnb2"))
    doParallel::registerDoParallel(cl)
    out <- foreach::foreach(sim.idx=c(1:nsim), .packages=c('quantreg','osqp','piqp','Matrix')) %dopar% {
  	  data.sim <- fit$model[idx[[sim.idx]],]
      fit.sim <- sqr(fit$formula,tau=fit$tau,tau0=fit$tau0,spar=fit$spar,method=fit$sqrmethod,
                      npar=fit$npar,ztol=fit$ztol,w=fit$w,all.knots=fit$all.knots,
					  solver=fit$solver,control=fit$control,data=data.sim)
      if(fit.sim$info != "solved") { fit.sim$coefficients <- fit.sim$coefficients * NA }
	  list(beta=fit.sim$coefficients,dbeta=fit.sim$derivatives,info=fit.sim$info)
    }
    parallel::stopCluster(cl)
	if(!mthreads) {
      ##sink(tempfile())
      RhpcBLASctl::blas_set_num_threads(Sys.getenv('OPENBLAS_NUM_THREADS'))
      ##sink()
    }
  } else {
    out <- foreach::foreach(sim.idx=c(1:nsim)) %do% {
      data.sim <- fit$model[idx[[sim.idx]],]
      fit.sim <- sqr(fit$formula,tau=fit$tau,tau0=fit$tau0,spar=fit$spar,method=fit$sqrmethod,
                     npar=fit$npar,ztol=fit$ztol,w=fit$w,solver=fit$solver,control=fit$control,data=data.sim)
	  if(fit.sim$info != "solved") { fit.sim$coefficients <- fit.sim$coefficients * NA }
	  list(beta=fit.sim$coefficients,dbeta=fit.sim$derivatives,info=fit.sim$info)
    }
  }
  beta.sim <- array(0,dim=c(nsim,dim(out[[1]]$beta)))
  dbeta.sim <- beta.sim
  info.sim <- rep(NA,nsim)
  for(sim.idx in c(1:nsim)) {
    beta.sim[sim.idx,,] <- out[[sim.idx]]$beta
	dbeta.sim[sim.idx,,] <- out[[sim.idx]]$dbeta
	info.sim[sim.idx] <- out[[sim.idx]]$info
  }
  return(list(coefficients=beta.sim,derivatives=dbeta.sim,info=info.sim))
}



#' Plot of Spline Quantile Regression Coefficients 
#'
#' This function plots one or all regression coefficients from quantile regression (QR) and spline quantile regression (SQR)
#' on a given sequence of quantiles against the quantile level.
#' @param summary.rq output of \code{summary()} for the QR fit from \code{quantreg::rq()}
#' @param summary.sqr output of \code{summary()} for the primary SQR fit from \code{sqr()} (defalty = \code{NULL})
#' @param summary.sqrb output of \code{summary()} for the secondary SQR fit from \code{sqr()} (default = \code{NULL})
#' @param coef.sim output \code{coefficients} of \code{boot.sqr()} for the primary SQR (default = \code{NULL})
#' @param type as \code{type} parameter in \code{plot()} for plotting SQR (default = \code{"l"})
#' @param lty line types for the primary and secondary SQR (default = \code{c(1,2)})
#' @param lwd line widths for the primary and secondary SQR (default = \code{c(1.5,1)})
#' @param cex as \code{cex} parameter in \code{plot()}
#' @param pch as \code{pch} parameter in \code{plot()} for the QR (default = 1)
#' @param col line colors for the primary and secondary SQR (default = \code{c(2,1)})
#' @param idx index of individual coefficient to be plotted (default = \code{NULL})
#' @param plot.rq \code{TRUE} (default) or \code{FALSE}: if \code{TRUE}, plot QR as points
#' @param plot.rq.line \code{TRUE} or \code{FALSE} (default): if \code{TRUE}, add line plot of QR
#' @param plot.zero \code{TRUE}  or \code{FALSE} (default): if \code{TRUE}, add zero line
#' @param plot.ls \code{TRUE} (default) or \code{FALSE}: if \code{TRUE}, add least-square estimate with 90-percent CI
#' @param plot.ci \code{TRUE} (default) or \code{FALSE}: if \code{TRUE}, add 90-percent bootstrap CI for the primary SQR 
#' when \code{coef.sim} is supplied or approximate CI when \code{coef.sim = NULL}
#' @param var.names user-supplied names of regression coefficients (including the intercept)
#' @param Ylim user-supplied matrix of \code{ylim} for each coefficient (default = \code{NULL})
#' @param xlim user-supplied \code{xlim} (default = \code{NULL})
#' @param set.par \code{TRUE} (default)  or \code{FALSE} if \code{TRUE}, reset \code{par()} 
#' @param mfrow parameter for resetting \code{par()} (default = \code{NULL})
#' @param lab parameter for resetting \code{par()} (default = \code{c(10,7,7)})
#' @param mar parameter for resetting \code{par()} (default = \code{c(2,3,2,1)+0.1})
#' @param las parameter for resetting \code{par()} (default = 1)
#' @return Ylim (invisible)
#' @export
sqr.plot <- function(summary.rq,summary.sqr=NULL,summary.sqrb=NULL,coef.sim=NULL,type="l",
  lty=c(1,2),lwd=c(1.5,1),cex=0.25,pch=1,col=c(2,1),idx=NULL,plot.rq=TRUE,plot.rq.line=FALSE,plot.zero=FALSE,
  plot.ls=TRUE,plot.ci=TRUE,var.names=NULL,Ylim=NULL,xlim=NULL,set.par=TRUE,
  mfrow=NULL,lab=c(10,7,7),mar=c(2,3,2,1)+0.1,las=1) {
  par0 <- par("mfrow","lab","mar","las")
  p <- dim(summary.rq[[1]]$coefficients)[1]
  nq <- length(summary.rq)
  # replace standard error by rq
  tau <- rep(0,nq)
  cf.rq <- matrix(0,p,nq)
  if(!is.null(summary.sqrb)) cf.sqrb<- matrix(0,p,nq)
  for(j in c(1:nq)) {
    tau[j] <- summary.rq[[j]]$tau
	cf.rq[,j] <- summary.rq[[j]]$coefficients[,1]
    if(!is.null(summary.sqrb)) summary.sqrb[[j]]$coefficients[,2] <- summary.rq[[j]]$coefficients[,2]
	if(!is.null(summary.sqrb)) cf.sqrb[,j] <- summary.sqrb[[j]]$coefficients[,1]
  }
  tmp.rq <- plot(summary.rq)
  if(!is.null(summary.sqr)) tmp <- plot(summary.sqr)
  if(is.null(mfrow) & is.null(idx))  { 
    mfrow <- c( ceiling(sqrt(p)), ceiling(sqrt(p)))
  } else {
    if(!is.null(idx) ) mfrow <- c(1,1)
  }
  if(is.null(lab)) lab <- c(10,7,7)
  if(is.null(mar)) mar=c(2,3,2,1)+0.1
  if(is.null(las)) las <- 1
  if(set.par) par(las=las,mar=mar,mfrow=mfrow, lab=lab)
  if(!is.null(idx)) {
    pp <- idx
  } else {
    pp <- c(1:p)
  }
  ylim <- matrix(NA,2,p)
  for(i in pp) {
    if(!is.null(summary.sqr)) {
      ci <- tmp[[1]][i,,]
	} else {
	  ci <- tmp.rq[[1]][i,,]
	}
	if(!is.null(coef.sim)) {
	  ci[2:3,1:(ncol(ci)-1)] <- apply(coef.sim[,i,],2,quantile,probs=c(0.05,0.95),na.rm=T)
    }
	if(is.null(Ylim)) {
	  if(!is.null(summary.sqr)) {
	    rg <- range(tmp$Ylim[,i],cf.rq[i,])
	  } else {
	    rg <- range(tmp.rq$Ylim[,i])
	  }
	} else {
	  rg <- Ylim[,i] 
	}
	ylim[,i] <- rg
    plot(tau,ci[1,1:nq],ylim=rg,type='n',ylab=NA,xlab=NA,xlim=xlim)
    tmp.x <- c(tau,rev(tau))
    tmp.y <- c(ci[3,1:nq],rev(ci[2,1:nq]))
    if(plot.ci) polygon(tmp.x, tmp.y, col = "grey",border=NA)
	if(plot.rq) {
	  points(tau,cf.rq[i,],cex=cex,pch=pch)
	  if(plot.rq.line) lines(tau,cf.rq[i,],lwd=0.5)
	}
    if(plot.ls) abline(h=ci[,nq+1],lty=c(2,3,3))
	if(!is.null(summary.sqrb)) lines(tau,cf.sqrb[i,],type=type,lwd=lwd[2],col=col[2],lty=lty[2])
	if(!is.null(summary.sqr)) lines(tau,ci[1,1:nq],type=type,lwd=lwd[1],col=col[1],lty=lty[1])
	if(is.null(var.names)) {
	  title(rownames(summary.rq[[1]]$coef)[i])
	} else {
	  title(var.names[i])
	}
	if(plot.zero) abline(h=0,lty=4)
  }
  if(set.par) par(mfrow = par0$mfrow, mar = par0$mar, lab=par0$lab, las=par0$las)
  Ylim <- ylim
  invisible(Ylim)
}




#' Plot of Derivative of Spline Quantile Regression Coefficients
#'
#' This function plots the derivative of one or all coefficients from spline quantile regression (SQR)
#' together with the corresponding first difference from quantile regression (QR) on a given sequence of 
#' quantiles against the quantile level.
#' @param fit output of \code{quantreg::rq()} for the QR fit
#' @param fit.sqr output of \code{sqr()} for the primary SQR fit
#' @param fit.sqrb output of \code{sqr()} for the secondary SQR fit (default = \code{NULL})
#' @param deriv.sim output \code{derivatives} of \code{boot.sqr()} for the primary SQR (default = \code{NULL})
#' @param type as \code{type} parameter in \code{plot()} for the primary SQR (default = \code{"l"})
#' @param typeb as \code{type} parameter in \code{plot()} for the secondary SQR (default = \code{"s"})
#' @param var.names user-supplied names of regression coefficients (including the intercept)
#' @param lty line types for the primary and secondary SQR (default = \code{c(1,2)})
#' @param lwd line widths for the primary and secondary SQR (default = \code{c(1.5,1)})
#' @param cex as \code{cex} parameter in \code{plot()}
#' @param pch as \code{pch} parameter in \code{plot()} for QR(default = 1)
#' @param col line colors for the primary and secondary SQR (default = \code{c(2,1)})
#' @param idx index of individual coefficient to be plotted (default = \code{NULL})
#' @param Ylim user-supplied matrix of \code{ylim} for each coefficient (default = \code{NULL})
#' @param xlim user-supplied \code{xlim} (default = \code{NULL})
#' @param plot.zero \code{TRUE}  or \code{FALSE} (default): if \code{TRUE}, add  zero line
#' @param plot.ci \code{TRUE} (default) or \code{FALSE}: if \code{TRUE}, add 90-percent bootstrap CI for the primary SQR
#'  when \code{deriv.sim} is supplied
#' @param set.par \code{TRUE} (default)  or \code{FALSE} if \code{TRUE}, reset \code{par()} 
#' @param mfrow parameter for resetting \code{par()} (default = \code{NULL})
#' @param lab parameter for resetting \code{par()} (default = \code{c(10,7,7)})
#' @param mar parameter for resetting \code{par()} (default = \code{c(2,3,2,1)+0.1})
#' @param las parameter for resetting \code{par()} (default = 1)
#' @return Ylim (invisible)
#' @export
sqr_deriv.plot<-function(fit,fit.sqr,fit.sqrb=NULL,deriv.sim=NULL,type="l",typeb="s",var.names=NULL,
  lty=c(1,2),lwd=c(1.5,1),cex=0.25,pch=1,col=c(2,1),idx=NULL,Ylim=NULL,xlim=NULL,plot.zero=TRUE,plot.ci=TRUE,
  set.par=TRUE,mfrow=NULL,lab=c(10,7,7),mar=c(2,3,2,1)+0.1,las=1) {
  par0 <- par("mfrow","lab","mar","las")
  p <- dim(fit$coefficients)[1]
  if(!is.null(idx)) {
    pp <- idx
  } else {
    pp <- c(1:p)
  }
  if(is.null(mfrow) & is.null(idx))  { 
    mfrow <- c( ceiling(sqrt(p)), ceiling(sqrt(p)))
  } else {
    if(!is.null(idx)) mfrow <- c(1,1)
  }
  if(is.null(lab)) lab <- c(10,7,7)
  if(is.null(mar)) mar=c(2,3,2,1)+0.1
  if(is.null(las)) las <- 1
  if(set.par) par(las=las,mar=mar,mfrow=mfrow, lab=lab)
  if(is.null(xlim)) xlim <- range(fit$tau)
  ylim <- matrix(NA,2,p)
  for(i in pp) {
	tmp <- diff(fit$coef[i,])/diff(fit$tau)
	rg <- range(tmp)
    if(!is.null(deriv.sim)) {
	  if(plot.ci) {
	    ci <- apply(deriv.sim[,i,],2,quantile,probs=c(0.05,0.95),na.rm=T)
		rg <- range(ci,rg)
	  }
	}
    if(!is.null(Ylim)) rg <- Ylim[,i]
    plot(fit$tau[-1],tmp,type='n',ylab=NA,xlab=NA,ylim=rg,pch=pch,cex=cex,xlim=xlim)
    if(!is.null(deriv.sim)) {
	  if(plot.ci) {
        tmp.x <- c(fit$tau,rev(fit$tau))
        tmp.y <- c(ci[2,],rev(ci[1,]))
        polygon(tmp.x, tmp.y, col = "grey",border=NA)
	  }
    }
    points(fit$tau[-1],tmp,pch=pch,cex=cex)
    lines(fit.sqr$tau,fit.sqr$derivatives[i,],lty=lty[1],col=col[1])
    if(!is.null(fit.sqrb)) {
	  tmp <- fit.sqrb$derivatives[i,]
	  tmp[length(tmp)] <- tmp[length(tmp)-1]
	  lines(fit.sqrb$tau,tmp,lty=lty[2],col=col[2],type=typeb)
	}
	if(is.null(var.names)) {
        title(rownames(fit$coef)[i])
    } else {
        title(var.names[i])
    }
	if(plot.zero) abline(h=0,lty=4)
	ylim[,i] <- rg
  }
  if(set.par) par(mfrow = par0$mfrow, mar = par0$mar, lab=par0$lab, las=par0$las)
  Ylim <- ylim
  invisible(Ylim)
}

Try the qfa package in your browser

Any scripts or data that you put into this service are public.

qfa documentation built on March 30, 2026, 5:07 p.m.