---
title: "Dimension Reduction Techniques for Conditional Quantiles in R: A Vignette"
author: 
- name: Eliana Christou
  affiliation: University of North Carolina at Charlotte
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Dimension Reduction Techniques for Conditional Quantiles in R: A Vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
abstract: "While dimension reduction techniques have been well developed, they focus on the entire conditional distribution of the data.  However, there are a lot of research areas where extremes are important and require the study of the conditional quantiles.  `quantdr` is an `R` package that performs dimension reduction to conditional quantiles by determining the directions that span the central quantile subspace (CQS)."  
keywords: "central quantile subspace, dimension reduction, quantile regression"
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

-------------

# Introduction

Data visualization and non-parametric fitting can become very challenging when data is *high-dimensional*.  This has led to the development of many dimension reduction techniques that focus on the entire conditional distribution of the data.  However, when specific aspects of the conditional distribution are of interest, these methods can provide more directions than necessary.  

Quantile regression (QR) has received a lot of attention since its inception from Koenker and Bassett (1978) and it has been an attractive tool for research areas where non-central parts of the data are important.  Dimension reduction techniques for conditional quantiles have recently made an appearance and are still in development.  For some references see Wu et al. (2010), Kong and Xia (2012, 2014), Luo et al. (2014), and Christou and Akritas (2016).

More recently, Christou (2020) proposed a new technique for finding the fewest linear combinations of the predictor variable $\mathbf{X}$ that contain all the information about the conditional quantile.  To give further intuition, for a univariate response $Y$ and a $p$-dimensional predictor $\mathbf{X}$, the method focuses on finding a $p \times d_{\tau}$ matrix $\mathbf{B}_{\tau}$, for $\tau \in (0, 1)$ and $d_{\tau} \leq p$, such that $Y$ and $Q_{\tau}(Y|\mathbf{X})$ are conditionally independent given $\mathbf{B}_{\tau}^{\top} \mathbf{X}$, where $Q_{\tau}(Y|\mathbf{X})$ denotes the $\tau$th conditional quantile of $Y$ given $\mathbf{X}$.  The space spanned by $\mathbf{B}_{\tau}$ is called the *$\tau$th quantile dimension reduction subspace* for the regression of $Y$ on $\mathbf{X}$.  The intersection of all $\tau$th dimension reduction subspaces is called the *$\tau$th central quantile subspace* ($\tau$-CQS) and is denoted by $\mathcal{S}_{Q_{\tau}(Y|\mathbf{X})}$.  

The purpose of this vignette is to demonstrate the implementation of the functions appearing in the `quantdr` package, along with some examples.  The main function of the package is `cqs` which returns the estimated directions of the $\tau$-CQS.  

# Getting Started

To get started you need to first install the package using the command 
```{r, eval=FALSE}
install.packages("quantdr")
```
and then use it during any `R` session by issuing the command
```{r setup}
library(quantdr)
```
For an overview of the available help files of the package use
```{r, eval=FALSE}
help(package = "quantdr")
```
The package includes the following functions:

* `llqr` Local linear quantile regression
* `llqrcv` Cross-Validation for bandwidth selection of local linear quantile regression
* `cqs` Central quantile subspace
* `ValAR` Value-at-Risk estimation using the central quantile subspace
<!-- * `bic_d` Modified-BIC type criterion for structural dimension selection -->

To get help on a specific function type `help(function.name)` or `?function.name` for a convenient shorthand.  Try
```{r, eval=FALSE}
help(cqs)
?cqs
```
For further examples, refer to the coding examples listed in the documentation for each `quantdr` function.  

# Central quantile subspace
## Overview

The overall goal of the `cqs` function is to identify the coefficients of the linear combination $\mathbf{B}_{\tau}^{\top}\mathbf{X}$ in order to replace the $p \times 1$ predictor vector $\mathbf{X}$ with the low-dimensional $d_{\tau} \times 1$ predictor vector $\mathbf{B}_{\tau}^{\top}\mathbf{X}$.  To do that, Christou (2020) proved two important results:

* $\boldsymbol{\beta}_{\tau} \in \mathcal{S}_{Q_{\tau}(Y|\mathbf{X})}$, where $\boldsymbol{\beta}_{\tau}$ is the slope vector from regressing the conditional quantile on $\mathbf{X}$, i.e., 
\begin{eqnarray*}
(\alpha_{\tau}, \boldsymbol{\beta}_{\tau}) = \arg \min_{(a_{\tau}, \mathbf{b}_{\tau})} E \{ Q_{\tau}(Y| \mathbf{A}^{\top} \mathbf{X}) - a_{\tau} - \mathbf{b}_{\tau}^{\top}\mathbf{X}\}^2,
\end{eqnarray*}
and $\mathbf{A}$ spans the central subspace (Li 1991).  This implies an initial dimension reduction using $\mathbf{A}$, making the nonparametric estimation of $Q_{\tau}(Y|\mathbf{A}^{\top}\mathbf{X})$ tractable.   

* $E\{Q_{\tau} (Y|U_{\tau})\mathbf{X}\} \in \mathcal{S}_{Q_{\tau}(Y|\mathbf{X})}$, where $U_{\tau}$ is a measurable function of $\mathbf{B}_{\tau}^{\top}\mathbf{X}$, provided that $Q_{\tau}(Y|U_{\tau})\mathbf{X}$ is integrable.  

The above two results imply the following:

1. If the dimension of the $\tau$-CQS is known to be one, then we can fit a linear regression model of $Q_{\tau}(Y|\mathbf{A}^{\top}\mathbf{X})$ on $\mathbf{X}$ and report the slope vector as the basis vector.

1. If the dimension of the $\tau$-CQS is known to be greater than one, then we can set $\boldsymbol{\beta}_{\tau, 0}=\boldsymbol{\beta}_{\tau}$ as the initial vector and then create more vectors using $\boldsymbol{\beta}_{\tau,j}=E\{Q_{\tau}(Y | \boldsymbol{\beta}_{\tau, j-1}^{\top}\mathbf{X})\mathbf{X}\}$, for $j=1, \dots, p-1$.  In order to obtain linearly independent vectors, we can perform an eigenvalue decomposition on $\mathbf{V}_{\tau}\mathbf{V}_{\tau}^{\top}$, where $\mathbf{V}_{\tau}=(\boldsymbol{\beta}_{\tau, 0}, \dots, \boldsymbol{\beta}_{\tau, p-1})$ and choose the eigenvectors corresponding to the $d_{\tau}$ nonzero eigenvalues.  

1. If the dimension of the $\tau$-CQS is unknown, we can apply the procedure from 2. Then, we can estimate the dimension $d_{\tau}$ using the eigenvalues resulting from the eigenvalue decomposition.  Existing techniques such as the Cross-Validation (CV) criterion or the modified-BIC type criterion of Zhu et al. (2010) can be used to determine the structural dimension. 

All of the above are incorporated in the `cqs` function.

## The algorithm at a glance

The first step of the algorithm requires to fit a linear regression model of $Q_{\tau}(Y|\mathbf{A}^{\top}\mathbf{X})$ on $\mathbf{X}$, which implies the non-parametric estimation of the conditional quantile.  The second step of the algorithm relies on estimating an expectation and performing an eigenvalue decomposition.  Specifically,

1.  We use a dimension reduction technique to estimate the basis matrix $\mathbf{A}$ of the central subspace, denoted by $\widehat{\mathbf{A}}$, and form the new sufficient predictors $\widehat{\mathbf{A}}^{\top} \mathbf{X}_{i}$, $i=1, \dots, n$.  For this we use the function `dr` from the package `dr`.

1.  For each $i=1, \dots, n$, we use the local linear conditional quantile estimation method of Guerre and Sabbah (2012) to estimate $Q_{\tau}(Y|\widehat{\mathbf{A}}^{\top} \mathbf{X}_{i})$.  For this we use the `llqr` of the presented package.

1.  We set $\widehat{\boldsymbol{\beta}}_{\tau}$ to be 
\begin{eqnarray*}
(\widehat{\alpha}_{\tau}, \widehat{\boldsymbol{\beta}}_{\tau}) = \arg \min_{(a_{\tau}, \mathbf{b}_{\tau})} \sum_{i=1}^{n} \{\widehat{Q}_{\tau}(Y|\widehat{\mathbf{A}}^{\top}\mathbf{X}_{i}) - a_{\tau} - \mathbf{b}_{\tau}^{\top} \mathbf{X}_{i}\}^2.
\end{eqnarray*}
For this we use the `lm` function of the `stats` package.

1. If $d_{\tau}=1$ we stop and report $\widehat{\boldsymbol{\beta}}_{\tau}$ as the estimated basis vector for $\mathcal{S}_{Q_{\tau}(Y|\mathbf{X})}$.  Otherwise, we move to Step 5.  

1. We set $\widehat{\boldsymbol{\beta}}_{\tau, 0}=\widehat{\boldsymbol{\beta}}_{\tau}$, and, for $j=1, \dots, p-1$, we
    A.  form the predictors $\widehat{\boldsymbol{\beta}}_{\tau, j-1}^{\top} \mathbf{X}_{i}$, $i=1, \dots, n$ and estimate $Q_{\tau}(Y| \widehat{\boldsymbol{\beta}}_{\tau, j-1}^{\top} \mathbf{X}_{i})$ using the local linear conditional quantile estimation method.  For this we use the `llqr` function of the presented paper.
    A.  form $\widehat{\boldsymbol{\beta}}_{\tau, j} = n^{-1} \sum_{i=1}^{n} \widehat{Q}_{\tau}(Y|\widehat{\boldsymbol{\beta}}_{\tau, j-1}^{\top} \mathbf{X}_{i}) \mathbf{X}_{i}$.
    
1.  We form the $p \times p$ matrix $\widehat{\mathbf{V}}_{\tau}=(\widehat{\boldsymbol{\beta}}_{\tau, 0}, \dots, \widehat{\boldsymbol{\beta}}_{\tau, p-1})$ and choose the $d_{\tau}$ eigenvectors corresponding to the $d_{\tau}$ largest eigenvalues of $\widehat{\mathbf{V}}_{\tau} \widehat{\mathbf{V}}_{\tau}^{\top}$.  For this we use the `eigen` function of the `base` package.  

If $d_{\tau}$ is unknown, a suggested structural dimension can be estimated using existing techniques such as the CV criterion or the modified-BIC type criterion of Zhu et al. (2010).  In terms of the `cqs` function, the modified-BIC type criterion is utilized and returns a suggested dimension.  Note that the user can extract the eigenvalues and apply them to any other existing technique.  

## Output of `cqs`

For the purpose of this section, we simulate data from the homoscedastic single-index model
\begin{eqnarray*}
Y=3X_{1}+X_{2}+\epsilon,
\end{eqnarray*}
where $\mathbf{X}=(X_{1}, \dots, X_{10})^{\top}$ and the error $\epsilon$ are generated according to a standard normal distribution.  The $\tau$-CQS is spanned by $(3, 1, 0, \dots, 0)^{\top}$ for all $\tau$.  We focus on the estimation of the 0.5-CQS.    
``` {r}
set.seed(1234)
n <- 100
p <- 10
tau <- 0.5
x <- matrix(rnorm(n * p), n, p)
error <- rnorm(n)
y <- 3 * x[, 1] + x[, 2] + error
```
The `cqs` function proceeds by specifying an $n \times p$ design matrix $\mathbf{X}$, a vector of the response variable $Y$, and a quantile level $\tau$.  The dimension of the $\tau$-CQS, i.e., $d_{\tau}$, is optional. Note that for this example $d_{\tau}=1$ for every quantile level $\tau$.  

We now discuss the output of the `cqs` function. 
```{r}
out1 <- cqs(x, y, tau = tau, dtau = 1)
out1
```
When the dimension of the $\tau$-CQS is known to be one, the algorithm reports the ordinary least-squares slope vector as the basis vector of the subspace.  This is denoted by `qvectors`.  Moreover, the function `cqs` outputs `dtau`, which in this case is specified by the user.

When the dimension of the $\tau$-CQS is unknown (or known to be greater than one), the algorithm continues, creates more vectors, and performs an eigenvalue decomposition.  Therefore, the output includes one more element, the eigenvalues provided by the eigenvalue decomposition and denoted by `qvalues`.
```{r}
out2 <- cqs(x, y, tau = tau)
out2
```
Note that the dimension `dtau` is correctly estimated by the modified-BIC type criterion.  Therefore, this output suggests to take the first eigenvector as the basis vector for the $\tau$-CQS.  We can extract the direction using
``` {r}
out2$qvectors[, 1:out2$dtau]
```
If we want to measure the estimation error we can use the angle between the true and the estimated subspaces.  This is performed using the `subspace` function of the `pracma` package.  Note that the angle is measured in radians, and so we divide by $\pi / 2$.
``` {r}
library(pracma)
beta_true <- c(3, 1, rep(0, p - 2))
beta_hat1 <- out1$qvectors
beta_hat2 <- out2$qvectors[, 1:out2$dtau]
subspace(beta_true, beta_hat1) / (pi / 2)
subspace(beta_true, beta_hat2) / (pi / 2)
```
Note that the estimation error is slightly higher when the dimension of the subspace is unknown and needs to be estimated. 

We can continue further and estimate the conditional quantile function.  To do this, we first form the new sufficient predictor $\widehat{\mathbf{B}}_{\tau}^{\top}\mathbf{X}$ using
``` {r}
newx <- x %*% beta_hat1
```
We then estimate the conditional quantile using the `llqr` function of the presented package.  The main arguments of the function are: the design matrix, which in this case is $\widehat{\mathbf{B}}_{\tau}^{\top}\mathbf{X}$, the vector of the response variable $Y$, and the quantile level $\tau$.  The rest of the arguments, i.e., the bandwidth $h$, the method for the estimation of $h$, and the single observation $x_0$ for which to perform the estimation, are optional.  If $h$ is not specified, then it will be determined using either the rule-of-thumb bandwidth of Yu and Jones (1998) or the CV criterion; the default choice is the rule-of-thumb.  However, the user needs to be careful about the bandwidth selection. When the dimension of the predictor variable is large compared to the sample size, local linear fitting meets the 'curse of dimensionality' problem. In situations like that, the bandwidth selected by the rule-of-thumb or the CV criterion might be small and cause the function to fail. For these cases, we advice the user to pre-specify the bandwidth in the function.  Finally, if $x_0$ is not specified, the estimation will be performed on the entire design matrix.  For this example we use
``` {r}
qhat1 <- llqr(newx, y, tau)
qhat1
```
The output consists of the estimation of the conditional quantile function at each point of the design matrix, i.e., $\widehat{Q}_{\tau}(Y|\widehat{\mathbf{B}}_{\tau}^{\top} \mathbf{X}_{i})$, $i=1, \dots, n$, and the estimation of the bandwidth using the rule-of-thumb.  For illustration purposes, we repeat the above using the CV criterion for the bandwidth selection.
``` {r}
qhat2 <- llqr(newx, y, tau, method = "CV")
qhat2
```
If $h$ is pre-specified, then the output reports the value given by the user, i.e.,
``` {r}
qhat3 <- llqr(newx, y, tau, h = 1)
qhat3
```
To illustrate the median regression, compared with the original data, we can use
``` {r fig1, fig.height = 4.5, fig.width = 4.5, fig.align = "center"}
true_dir <- x %*% beta_true

plot(true_dir, y, xlab = "sufficient direction", ylab = "y", pch = 16)
points(true_dir, qhat1$ll_est, pch = 16, col = 'red')
```

The same procedure can be performed for multiple quantile levels.  For example, we estimate the $\tau$-CQS for $\tau= 0.1, 0.25, 0.5, 0.75, 0.9$.  
``` {r}
taus <- c(0.1, 0.25, 0.5, 0.75, 0.9)

out3 <- matrix(0, p, length(taus))
for (i in 1:length(taus)) {
  out3[, i] <- cqs(x, y, tau = taus[i], dtau = 1)$qvectors
}
out3
```
Similarly, we can plot the data along with the estimated conditional quantile functions for the various quantile levels.  
``` {r fig2, fig.height = 5.5, fig.width = 7, fig.align = "center"}
newx <- x %*% out3

oldpar <- par(no.readonly = TRUE)
par(mfrow=c(2,3))
qhat_tau <- as.null()
for (i in 1:length(taus)) {
  plot(true_dir, y, xlab = "sufficient direction", ylab = "y", main = taus[i], pch = 16)
  qhat_tau <- llqr(newx[, i], y, tau = taus[i])$ll_est
  points(true_dir, qhat_tau, pch = 16, col = "red")
}
par(oldpar)
```

<!--  ## A Real Data Application -->

# Value-at-Risk Estimation
## Overview

Value-at-risk (VaR) is an important financial quantity directly related to the QR framework.  Specifically, for $\{r_{t}\}_{t=1}^{n}$ a time series of returns, the $\tau$th VaR at time $t$, denoted by $VaR_{\tau}(t)$, is the smallest number for which $P\{r_{t} < -VaR_{\tau}(t) | \mathcal{F}_{t-1}\} = \tau$, where $\mathcal{F}_{t-1}$ denotes the information set at time $t-1$.  Therefore, $-VaR_{\tau}(t)$ is the $\tau$th conditional quantile of $r_{t}$, and we can write $Q_{\tau}(r_{t} | \mathcal{F}_{t-1}) = -VaR_{\tau}(t)$.   

Christou and Grabchak (2019) used the methodology proposed by Christou (2020) to estimate VaR and compare it with existing commonly-used methods.  The estimation of $VaR_{\tau}(t)$ at time $t$ reduces to the estimation of $Q_{\tau}(r_{t} | \mathcal{F}_{t-1})$, which can be performed using the `cqs` and `llqr` functions of the presented package.  Specifically, let $r_{t}$ denote a daily return (usually log-return) and $\mathbf{r}_{t}$ denote a $p$-dimensional vector of covariates consisting of the $p$ past returns, i.e., $\mathbf{r}_{t,p} = (r_{t-1},\dots, r_{t-p})^{\top}$.  We want to estimate $Q_{\tau}(r_{t} | \mathbf{r}_{t,p})$ at time $t$ by assuming that there exists a $p \times d_{\tau}$ matrix $\mathbf{B}_{\tau}$, $d_{\tau} \leq p$, such that $Q_{\tau}(r_{t} | \mathbf{r}_{t,p}) = Q_{\tau}(r_{t} | \mathbf{B}_{\tau}^{\top}\mathbf{r}_{t,p})$.  This implies that $\mathbf{B}_{\tau}^{\top}\mathbf{r}_{t,p}$ contains all the information about $r_{t}$ that is available from $Q_{\tau}(r_{t} | \mathbf{r}_{t,p})$.  Once $\mathbf{B}_{\tau}$ is estimated using the `cqs` function, the sufficient predictors $\widehat{\mathbf{B}}_{\tau}^{\top}\mathbf{r}_{t,p}$ are formed and $Q_{\tau}(r_{t} | \widehat{\mathbf{B}}_{\tau}^{\top}\mathbf{r}_{t,p})$ is estimated using the `llqr` function.    

## Output of `ValAR`

For an illustration of VaR estimation, we use the `edhec` data set available in the `PerformanceAnalytics` package in R.  The data set includes the EDHEC composite hedge fund style index returns.  
```{r, message = FALSE}
library(PerformanceAnalytics)
data(edhec, package = "PerformanceAnalytics")
head(edhec)
```
To estimate the one-step ahead VaR we need a vector of returns in a standard chronological order.  For this, we will select the 'CTA Global Distressed' as a vector of returns.  The `ValAR` function proceeds by specifying the vector of returns $y$, the number $p$ of past observations to be used as the predictor variables, the quantile level $\tau$, a moving window for the number of observations to be used for fitting the model, and a logical operator to indicate whether the returns are given in a standard chronological order (oldest to newest) or not.  If a moving window is not specified, then the default value of min(250, n - p) will be used to estimate the conditional quantile.  Note that, typical values for a moving window correspond to one or two years of returns.  A moving window will be useful if the number of observations is large and old returns will not be of much value to the prediction.  Moreover, if the number of observations is large, a moving window will also speed up the algorithm, since the number of observations to be used will be smaller.  Finally, if the returns are in reversed chronological order, then the user should either reverse the vector of returns or indicate that using the `chronological` option of the function. 

For this example, we will use $p = 5$, $\tau = 0.05$, and will not specify any moving window so that the default value of min(250, 275 - 5) = 250 will be used.  Also, the returns are given in standard chronological order and, since this is the default setting, we do not specify anything for `chronological`.  Note that, the last day of returns is 2019-11-30 and therefore, we are predicting the 5\% VaR for 2019-12-31.     

```{r}
y <- as.vector(edhec[, 2])
n <- length(y)
p <- 5
tau <- 0.05
ValAR(y, p = p, tau = tau)
```

If we want to compare this number with the 5\% VaR estimated by the historical method, we can use
```{r}
VaR(y, 0.95, method = 'historical')
```

For illustration purposes, we will use the last 100 observations as a testing data set for one-step ahead VaR estimation and the rest of the observations as historical data.  Specifically, for each of the time points, $t$, of the observations not included in the historical data, we estimate the model based on all of the data up to time $t-1$ and report the 5\% VaR.     
```{r, eval=FALSE}
size <- 100

VaRest <- as.null(size)
for (i in 1:size){
  VaRest[i] <- ValAR(y[1:(n - size + i - 1)], p, tau)
}
```
To visualize the results, we can plot the time series of returns with the one-step ahead 5\% VaR forecasts.  
```{r fig3, fig.height = 4.8, fig.width = 6.5, fig.align = "center", eval=FALSE}
plot.ts(y[(n - size + 1):n], ylim = range(y[(n - size + 1):n], VaRest), ylab = 'returns')
lines(VaRest, col = 'red')
```

```{r fig.height = 4.8, fig.width = 6.5, include=TRUE, fig.align = "center", echo = FALSE}
img <- png::readPNG("returns1.png")
grid::grid.raster(img)
```

Moreover, we can calculate the proportion of exceptions, i.e., the number of times $r_{t} < -\widehat{VaR}_{\tau}(t)$ for each time point $t$ of the last 100 observations.  This number should be a reasonable estimate for $\tau$.  In this case we can see that the estimated proportion is actually 0.05.   
```{r, eval=FALSE}
sum(y[(n - size + 1):n] < VaRest) / size 
#> [1] 0.05
```

We can repeat the above procedure for multiple quantile levels, i.e., for $\tau= 0.01, 0.025, 0.05$.
```{r fig4, fig.height = 4.8, fig.width = 6.5, fig.align = "center", eval=FALSE}
taus <- c(0.01, 0.025, 0.05)
VaRest <- matrix(0, size, length(taus))
for (i in 1:size) {
  for (j in 1:length(taus)) {
  VaRest[i, j] <- ValAR(y[1:(n - size + i - 1)], p, taus[j])
  }
}  

# plots
plot.ts(y[(n - size + 1):n], ylim = range(y[(n - size + 1):n], VaRest), ylab = 'returns')
lines(VaRest[, 1], col = 'red')
lines(VaRest[, 2], col = 'blue')
lines(VaRest[, 3], col = 'green')
legend('top', legend = c("1%", "2.5%", "5%"), col = c("red", "blue", "green"),
       lty=1, cex=1, horiz = T, bty = "n")
```

```{r fig.height = 4.8, fig.width = 6.5, include=TRUE, fig.align = "center", echo = FALSE}
img <- png::readPNG("returns2.png")
grid::grid.raster(img)
```

```{r, eval=FALSE}
# proportion of exceptions
sum(y[(n - size + 1):n] < VaRest[, 1]) / size 
#> [1] 0.03
sum(y[(n - size + 1):n] < VaRest[, 2]) / size 
#> [1] 0.03
sum(y[(n - size + 1):n] < VaRest[, 3]) / size 
#> [1] 0.05
```
We can see that the estimated proportion of exceptions are very close to the true values of 0.01, 0.025, and 0.05.  


<!-- ## A Real Data Application -->


<!-- !!!!!!!! If you uncomment the following in a future version of the package, remember to include BatchGetSymbols in Suggests of Description!!!!!!!!!!!!! -->

<!-- ## A Real Data Application -->

<!-- In this section, we demonstrate an example presented in Christou (2020) concerning the Bitcoin data set (Nakamoto 2008), a cryptocurrency characterized by great volatility.  The data can be downloaded in `R` from [https://finance.yahoo.com/](https://finance.yahoo.com/) using the `BatchGetSymbols` package.  The data consists of the daily returns from January 2016 to December 2018. -->
<!-- ``` {r, message=FALSE, warning = FALSE} -->
<!-- library(BatchGetSymbols) -->
<!-- BTC <- BatchGetSymbols("BTC-USD", first.date = "2016-01-01", last.date = "2018-12-31") -->
<!-- BTC <- BTC$df.tickers -->

<!-- names(BTC) -->
<!-- dim(BTC) -->
<!-- ``` -->
<!-- We consider $Y$ as the log return and $\mathbf{X}=(X_{1}, \dots, X_{10})^{\top}$, where $X_{1}, \dots, X_{5}$ are the returns in the past five days, and $X_{6}, \dots, X_{10}$ are the absolute values of those returns. -->
<!-- ``` {r} -->
<!-- y <- BTC$price.adjusted -->
<!-- r <- diff(log(y)) -->
<!-- lag <- 5 -->
<!-- newy <- r[ - c(1:lag)] -->
<!-- n <- length(newy) -->
<!-- p <- 10 -->

<!-- xtemp <- matrix(0, n, lag) -->
<!-- for (i in (lag + 1):(n + lag)){ -->
<!-- xtemp[i - lag, ] <- r[(i - lag):(i - 1)] -->
<!-- } -->
<!-- xtemp2 <- cbind(xtemp[, 5], xtemp[, 4], xtemp[, 3], xtemp[, 2], xtemp[, 1]) -->

<!-- x <- cbind(xtemp2, abs(xtemp2)) -->
<!-- dim(x) -->
<!-- ``` -->
<!-- We want to estimate $Q_{\tau}(Y|\mathbf{X})$ for $\tau=0.01, 0.025, 0.05$.  The reason for choosing these values is because $Q_{\tau}(Y|\mathbf{X})$ is an important quantity in a financial risk setup, called the value-at-risk (VaR). -->

<!-- To do this, we first want to reduce the dimension of the predictor vector $\mathbf{X}$, and therefore, we estimate the $\tau$-CQS.  Since the dimension of each $\tau$-CQS is unknown, we do not specify `dtau` and instead, allow the function to estimate it using the modified-BIC type criterion.  This is done by -->
<!-- ``` {r} -->
<!-- taus <- c(0.01, 0.025, 0.05) -->
<!-- out1 <- cqs(x, newy, tau = taus[1]) -->
<!-- out2 <- cqs(x, newy, tau = taus[2]) -->
<!-- out3 <- cqs(x, newy, tau = taus[3]) -->

<!-- out1$dtau; out2$dtau; out3$dtau -->
<!-- ``` -->
<!-- Since the dimension of the $\tau$-CQS for all quantile levels is estimated to be one, we take the first eigenvector as its estimated basic vector. -->
<!-- ``` {r} -->
<!-- beta_hat1 <- out1$qvectors[, 1] -->
<!-- beta_hat2 <- out2$qvectors[, 1] -->
<!-- beta_hat3 <- out3$qvectors[, 1] -->
<!-- beta_hat <- cbind(beta_hat1, beta_hat2, beta_hat3) -->
<!-- ``` -->
<!-- To estimate the conditional quantile function, we form the new sufficient predictors $\widehat{\mathbf{B}}_{\tau}^{\top}\mathbf{X}$ and use the `llqr` function. -->
<!-- ``` {r fig3, fig.height = 5.5, fig.width = 7, fig.align = "center"} -->
<!-- newx <- x %*% beta_hat -->
<!-- qhat <- matrix(0, n, length(taus)) -->
<!-- for (i in 1:length(taus)) { -->
<!--   qhat[, i] <- llqr(newx[, i], newy, tau = taus[i])$ll_est -->
<!-- } -->

<!-- plot.ts(newy, ylab = "returns") -->
<!-- lines(qhat[, 1], col = "red") -->
<!-- lines(qhat[, 2], col = "blue") -->
<!-- lines(qhat[, 3], col = "green") -->
<!-- legend(840, 0.22, c("tau=0.01","tau=0.025","tau=0.05"), col=c("red", "blue", "green"), lty = 1) -->
<!-- ``` -->

<!-- To further illustrate the estimation of the conditional quantile function, we estimate the proportion of times the return is less than $\widehat{Q}_{\tau}(Y|\widehat{\mathbf{B}}_{\tau}^{\top}\mathbf{X})$, i.e., the proportion of times $Y < \widehat{Q}_{\tau}(Y|\widehat{\mathbf{B}}_{\tau}^{\top}\mathbf{X})$. -->
<!-- ``` {r} -->
<!-- colMeans(matrix(newy, dim(qhat)[1], dim(qhat)[2]) < qhat) -->
<!-- taus -->
<!-- ``` -->
<!-- Note that the estimated proportions are very close to the true values of $\tau$. -->

# Conclusion

This vignette provides a brief introduction to the R package `quantdr` and presents a tutorial on how to implement the basic function `cqs`.  Performing dimension reduction techniques to conditional quantiles is an active research topic and therefore updates and new functions will be incorporated into forthcoming versions of the package.  For this reason, this document will be updated accordingly and will be made available through the package.  Suggestions and recommendations from the users will be highly appreciated and welcomed.  This package is also available through [GitHub](https://github.com/elianachristou/quantdr).  Feel free to contact the author at echris15@uncc.edu.  

# References

1. Christou, E. (2020) Central quantile subspace. *Statistics and Computing*, 30, 677–695

1. Christou, E. and Akritas, M. (2016) Single index quantile regression for heteroscedastic data. *Journal of Multivariate Analysis*, 150, 169-182

1.  Christou, E. and Grabchak, M. (2019) Estimation of value-at-risk using single index quantile regression.  *Journal of Applied Statistics*, 46(13), 2418–2433

1. Guerre, E., and Sabbah, C. (2012) Uniform bias study and Bahadur representation for local polynomial estimators of the conditional quantile function. *Econometric Theory*, 28, 87-129

1. Koenker, R., and Bassett, G. (1978) Regression quantiles. *Econometrica*, 46, 33-50

1. Kong, E., and Xia, Y. (2012) A single-index quantile regression model and its estimation. *Econometric Theory*, 28, 730-768

1. Kong, E., and Xia, Y. (2014) An adaptive composite quantile approach to dimension reduction. *Annals of Statistics*, 42, 1657-1688

1. Li, K.-C. (1991) Sliced inverse regression for dimension reduction. *Journal of the American Statistical Association*, 86, 316-327

1. Luo, W., Li, B., and Yin, X. (2014) On efficient dimension reduction with respect to a statistical functional of interest. *Annals of Statistics*, 42, 382-412 
<!-- 1. Nakamoto, S. (2008) Bitcoin: a peer-to-peer electronic cash system. Available online. [https://bitcoin.org/bitcoin.pdf](https://bitcoin.org/bitcoin.pdf)  -->
1. Wu, T. Z., Yu, K., and Yu, Y. (2010) Single index quantile regression. *Journal of Multivariate Analysis*, 101, 1607-1621

1. Yu, K., and Jones, M. C. (1998) Local linear quantile regression. *Journal of the American Statistical Association*, 93, 228-238

1. Zhu, L.-P., Zhu, L.-X., and Feng, Z.-H. (2010) Dimension reduction in regression through cumulative slicing estimation. *Journal of the American Statistical Association*, 105, 1455-1466
