This package contains two main modules. The linear module implements convolution-smoothed quantile regression in both low and high dimensions, as well as joint linear quantile and expected shortfall (ES) regression. The nonlinear module contains three classes to implement joint nonlinear/nonparametric quantile and ES regression. This package is also avaiable at https://pypi.org/project/quantes/.
The low_dim class in the linear module applies a convolution smoothing approach to fit linear quantile regression models, known as conquer. It also constructs normal-based and (multiplier) bootstrap confidence intervals for all slope coefficients. The high_dim class fits sparse quantile regression models in high dimensions via L1-penalized and iteratively reweighted L1-penalized (IRW-L1) conquer methods. The IRW method is inspired by the local linear approximation (LLA) algorithm proposed by Zou & Li (2008) for folded concave penalized estimation, exemplified by the SCAD penalty (Fan & Li, 2001) and the minimax concave penalty (MCP) (Zhang, 2010). Computationally, each weighted l1-penalized conquer estimator is solved using the local adaptive majorize-minimization algorithm (LAMM). For comparison, the proximal ADMM algorithm (pADMM) is implemented in the admm module. The joint class fits joint linear quantile and ES regression models (Dimitriadis & Bayer, 2019; Patton, Ziegel & Chen, 2019) using either FZ loss minimization (Fissler & Ziegel, 2016) or two-step procedures (Barendse, 2020; Peng & Wang, 2023; He, Tan & Zhou, 2023). For the second step of ES estimation, setting robust=TRUE uses the Huber loss with an adaptively chosen robustification parameter to gain robustness against heavy-tailed error/response; see He, Tan & Zhou (2023) for more details. Moreover, a combination of the iteratively reweighted least squares (IRLS) algorithm and quadratic programming is used to compute non-crossing ES estimates. This ensures that the fitted ES does not exceed the fitted quantile at each observation.
Estimation and inference for sparse ES regression are currently under construction (Zhang et al., 2023).
The nonlinear modules contains three classes, KRR, LocPoly and FNN, which implement three nonparametric joint quantile and expected shortfall estimators using kernel ridge regression (Yu et al., 2024a), local polynomial regression (Olma, 2021), and feedforward neural network regression (Yu et al. 2024b), respectively. For fitting nonparametric QR through the qt() method in both KRR and FNN, there is a smooth option available. When set to TRUE, it uses the Gaussian kernel convoluted check loss. For fitting nonparametric ES regression using (nonparametrically) generated surrogate response variables, the es() method in FNN provides two options: squared loss (robust=FALSE) and the Huber loss (robust=TRUE); the res() method in KRR fits robust ES kernel ridge regression using the Huber loss. Currently, the es() method in the LocPoly module uses only the squared loss to compute the two-step local polynomial ES regression estimator.
python >= 3.9, numpy, matplotlib, scipy, scikit-learn, cvxopt, qpsolvers, torch
optional: osqp
import numpy as np
import numpy.random as rgt
from scipy.stats import t
from quantes.linear import low_dim, high_dim
from quantes.admm import proximal
Generate data from a linear model with random covariates. The dimension of the feature/covariate space is p, and the sample size is n. The itercept is 4, and all the p regression coefficients are set as 1 in magnitude. The errors are generated from the t2-distribution (t-distribution with 2 degrees of freedom), centered by subtracting the population τ-quantile of t2.
When p < n, the module low_dim constains functions for fitting linear quantile regression models with uncertainty quantification. If the bandwidth h is unspecified, the default value max{0.01, {τ(1- τ)}^0.5 {(p+log(n))/n}^0.4} is used. The default kernel function is Laplacian. Other choices are Gaussian, Logistic, Uniform and Epanechnikov.
n, p = 8000, 400
itcp, beta = 4, np.ones(p)
tau, t_df = 0.75, 2
X = rgt.normal(0, 1.5, size=(n,p))
Y = itcp + X.dot(beta) + rgt.standard_t(t_df, n) - t.ppf(tau, t_df)
qr = low_dim(X, Y, intercept=True).fit(tau=tau)
qr.keys()
# qr['beta'] : conquer estimate (intercept & slope coefficients)
# qr['res'] : n-vector of fitted residuals
# qr['bw'] : bandwidth
# qr['message']: convergence message
At each quantile level τ, the norm_ci and boot_ci methods provide four 100* (1-alpha)% confidence intervals (CIs) for regression coefficients: (i) normal distribution calibrated CI using estimated covariance matrix, (ii) percentile bootstrap CI, (iii) pivotal bootstrap CI, and (iv) normal-based CI using bootstrap variance estimates. For multiplier/weighted bootstrap implementation with boot_ci, the default weight distribution is Exponential. Other choices are Rademacher, Multinomial (Efron's nonparametric bootstrap), Gaussian, Uniform and Folded-normal. The latter two require a variance adjustment; see Remark 4.7 in Paper.
n, p = 500, 20
itcp, beta = 4, np.ones(p)
tau, t_df = 0.75, 2
X = rgt.normal(0, 1.5, size=(n,p))
Y = itcp + X.dot(beta) + rgt.standard_t(t_df, n) - t.ppf(tau, t_df)
qr = low_dim(X, Y, intercept=True)
ci1 = qr.norm_ci(tau=tau)
ci2 = qr.mb_ci(tau=tau)
# ci1['normal'] : p+1 by 2 numpy array of normal CIs based on estimated asymptotic covariance matrix
# ci2['percentile'] : p+1 by 2 numpy array of bootstrap percentile CIs
# ci2['pivotal'] : p+1 by 2 numpy array of bootstrap pivotal CIs
# ci2['normal'] : p+1 by 2 numpy array of normal CIs based on bootstrap variance estimates
The module high_dim contains functions that fit high-dimensional sparse quantile regression models through the LAMM algorithm. The default bandwidth value is max{0.05, {τ(1- τ)}^0.5 { log(p)/n}^0.25}. To choose the penalty level, the self_tuning method implements the simulation-based approach proposed by Belloni & Chernozhukov (2011).
The l1 and irw methods compute L1- and IRW-L1-penalized conquer estimators, respectively. For the latter, the default concave penality is SCAD with constant a=3.7 (Fan & Li, 2001). Given a sequence of penalty levels, the solution paths can be computed by l1_path and irw_path.
p, n = 1028, 256
tau = 0.8
itcp, beta = 4, np.zeros(p)
beta[:15] = [1.8, 0, 1.6, 0, 1.4, 0, 1.2, 0, 1, 0, -1, 0, -1.2, 0, -1.6]
X = rgt.normal(0, 1, size=(n,p))
Y = itcp + X@beta + rgt.standard_t(2,size=n) - t.ppf(tau,df=2)
sqr = high_dim(X, Y, intercept=True)
lambda_max = np.max(sqr.tuning(tau))
lambda_seq = np.linspace(0.25*lambda_max, lambda_max, num=20)
## l1-penalized conquer
l1_model = sqr.l1(tau=tau, Lambda=0.75*lambda_max)
## iteratively reweighted l1-penalized conquer (default is SCAD penality)
irw_model = sqr.irw(tau=tau, Lambda=0.75*lambda_max)
## solution path of l1-penalized conquer
l1_path = sqr.l1_path(tau=tau, lambda_seq=lambda_seq)
## solution path of irw-l1-penalized conquer
irw_path = sqr.irw_path(tau=tau, lambda_seq=lambda_seq)
## model selection via bootstrap
boot_model = sqr.boot_select(tau=tau, Lambda=0.75*lambda_max, weight="Multinomial")
print('selected model via bootstrap:', boot_model['majority_vote'])
print('true model:', np.where(beta!=0)[0])
The module admm contains two classes, proximal and ncvx. The former employs the proximal ADMM algorithm to solve weighted L1-penalized quantile regression (without smoothing).
beta_true = np.insert(beta, 0, itcp)
lambda_max = np.max(high_dim(X, Y, intercept=True).tuning(tau))
lambda_seq = np.linspace(0.25*lambda_max, lambda_max, num=20)
admm = proximal(X, Y, intercept=True)
## l1-penalized QR
l1_admm = admm.l1(tau=tau, Lambda=0.5*lambda_max)
print(np.mean((l1_admm['beta'] - beta_true)**2).round(5))
print()
## iteratively reweighted l1-penalized QR (default is SCAD penality)
irw_admm = admm.irw(tau=tau, Lambda=0.75*lambda_max)
print(np.mean((irw_admm['beta'] - beta_true)**2).round(5))
print()
## solution path of l1-penalized QR
l1_admm_path = admm.l1_path(tau=tau, lambda_seq=lambda_seq)
print(np.mean((l1_admm_path['beta_seq'] - beta_true.reshape(-1,1))**2, axis=0).round(5))
print()
## solution path of irw-l1-penalized QR
irw_admm_path = admm.irw_path(tau=tau, lambda_seq=lambda_seq)
print(np.mean((irw_admm_path['beta_seq'] - beta_true.reshape(-1,1))**2, axis=0).round(5))
print()
The joint class in quantes.linear contains functions that fit joint (linear) quantile and expected shortfall models. The fz method computes joint quantile and ES regression estimates based on FZ loss minimization (Fissler & Ziegel, 2016). The twostep method implements two-stage procesures to compute quantile and ES regression estimates, with the ES part depending on a user-specified loss. Options are L2, TrunL2, FZ and Huber. The nc method computes non-crossing counterparts of the ES estimates when loss = L2 or Huber.
import numpy as np
import pandas as pd
import numpy.random as rgt
from quantes.linear import joint
p, n = 10, 5000
tau = 0.1
beta = np.ones(p)
gamma = np.r_[0.5*np.ones(2), np.zeros(p-2)]
X = rgt.uniform(0, 2, size=(n,p))
Y = 2 + X @ beta + (X @ gamma) * rgt.normal(0, 1, n)
lm = joint(X, Y)
## two-step least squares
m1 = lm.twostep(tau=tau, loss='L2')
## two-step truncated least squares
m2 = lm.twostep(tau=tau, loss='TrunL2')
## two-step FZ loss minimization
m3 = lm.twostep(tau=tau, loss='FZ', G2_type=1)
## two-step adaptive Huber
m4 = lm.twostep(tau=tau, loss='Huber')
## non-crossing two-step least squares
m5 = lm.nc(tau=tau, loss='L2')
## non-crossing two-step adaHuber
m6 = lm.nc(tau=tau, loss='Huber')
## joint quantes regression via FZ loss minimization (G1=0)
m7 = lm.fz(tau=tau, G1=False, G2_type=1, refit=False)
## joint quantes regression via FZ loss minimization (G1(x)=x)
m8 = lm.fz(tau=tau, G1=True, G2_type=1, refit=False)
out = pd.DataFrame(np.c_[(m1['coef_e'], m2['coef_e'], m3['coef_e'], m4['coef_e'],
m5['coef_e'], m6['coef_e'], m7['coef_e'], m8['coef_e'])],
columns=['L2', 'TLS', 'FZ', 'AH', 'NC-L2', 'NC-AH', 'Joint0', 'Joint1'])
out
The KRR class in quantes.nonlinear contains functions to compute quantile and expected shortfall kernel ridge regression (KRR) estimators, respectively. Recast as a quadratic program (Takeuchi et al., 2006), the qt() method computes quantile KRR using QP solvers from qpsolvers (Caron et al., 2024). Setting smooth=TRUE, a convolution-smoothed quantile KRR is computed by the BFGS algorithm using scipy.optimize.minimize. By plugging in estimated quantiles, the two-step ES KRR estimator and its robust counterpart (using the Huber loss) are computed by the es() method with robust=FALSE and robust=TRUE, respectively. The methods ls(), qt() and es() from the FNN class in quantes.nonlinear compute nonparametric mean, quantile and expected shortfall regression estimators using feedforward neural networks, respectively.
import numpy as np
import numpy.random as rgt
from scipy.stats import norm
from quantes.nonlinear import KRR, LocPoly
m_fn = lambda x: np.cos(2*np.pi*(x[:,0])) \
+ (1 + np.exp(-x[:,1]-x[:,2]))**(-1) + (1 + x[:,3] \
+ x[:,4])**(-3) + (x[:,5] + np.exp(x[:,6]*x[:,7]))**(-1)
s_fn = lambda x: np.sin(np.pi*(x[:,0] + x[:,1])*0.5) \
+ np.log(1 + (x[:,2]*x[:,3]*x[:,4])**2) \
+ x[:,7]*(1 + np.exp(-x[:,5]-x[:,6]))**(-1)
n, p = 2048, 8
X = rgt.uniform(0, 1, (n, p))
Y = m_fn(X) + s_fn(X)*rgt.normal(0, 1, n)
tau = 0.2
qt = norm.ppf(tau)
es = norm.expect(lambda x : (x if x < qt else 0))/tau
X_test = rgt.uniform(0, 1, (1024, p))
qt_test = m_fn(X_test) + s_fn(X_test)*qt
es_test = m_fn(X_test) + s_fn(X_test)*es
# kernel ridge regression (KRR)
kr = KRR(X, Y, kernel='polynomial',
kernel_params={'degree': 3, 'gamma': 1, 'coef0': 1})
kr.qt(tau=tau, alpha=.5, solver='cvxopt')
qt_pred = kr.qt_predict(X_test)
print('Mean squared prediction error of quantile-KRR:', np.mean((qt_test - qt_pred)**2))
kr.ES(tau=tau, alpha=2, x=X_test)
es_pred = kr.pred_e
print('Mean squared prediction error of ES-KRR:', np.mean((es_test - es_pred)**2))
print()
# local linear regression (LLR)
lp = LocPoly(X, Y) # default kernel is Gaussian
qt_fit = lp.qt_predict(x0=X, tau=tau, bw=0.5, degree=1)
qt_pred = lp.qt_predict(x0=X_test, tau=tau, bw=0.5, degree=1)
print('Mean squared prediction error of quantile-LLR:', np.mean((qt_test - qt_pred)**2))
es_pred = lp.es(x0=X_test, tau=tau, bw=1, degree=1, fit_q=qt_fit)
print('Mean squared prediction error of ES-LLR:', np.mean((es_test - es_pred)**2))
Fernandes, M., Guerre, E. and Horta, E. (2021). Smoothing quantile regressions. J. Bus. Econ. Statist. 39(1) 338–357. Paper
He, X., Pan, X., Tan, K. M. and Zhou, W.-X. (2023). Smoothed quantile regression with large-scale inference. J. Econom. 232(2) 367-388. Paper
He, X., Tan, K. M. and Zhou, W.-X. (2023). Robust estimation and inference for expected shortfall regression with many regressors. J. R. Stat. Soc. B. 85(4) 1223-1246. Paper
Koenker, R. (2005). Quantile Regression. Cambridge University Press, Cambridge. Book
Pan, X., Sun, Q. and Zhou, W.-X. (2021). Iteratively reweighted l1-penalized robust regression. Electron. J. Stat. 15(1) 3287-3348. Paper
Tan, K. M., Wang, L. and Zhou, W.-X. (2022). High-dimensional quantile regression: convolution smoothing and concave regularization. J. R. Stat. Soc. B. 84(1) 205-233. Paper
Yu, M., Wang, Y., Xie, S., Tan, K. M. and Zhou, W.-X. (2024a). Estimation and inference for nonparametric expected shortfall regression over RKHS. Paper
Yu, M., Tan, K. M., Wang, H. J. and Zhou, W.-X. (2024b). Deep neural expected shortfall regression with tail-robustness.
Zhang, S., He, X., Tan, K. M. and Zhou, W.-X. (2023). High-dimensional expected shortfall regression. Paper
This package is released under the GPL-3.0 license.