Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
64,192 characters · 23 sections · 78 citation commands
Factor Model of Mixtures
This paper proposes a new approach to estimating the distribution of the response variable $y$ conditioned on observing some factors $\bm{x}$ given data pairs $\{(y_i,\bm{x}_i)\}_{i=1}^N$. The quantile function of random variable $y$ is defined by $ Q_y(p) = \inf_q\{ p \leq P(y\leq q) \} $. The conditional quantile function $Q_{y|\bm{x}}(p)$ is modeled by a mixture (weighted sum) of basis quantile functions, e.g. quantile functions of normal distribution and exponential distribution. The weight of each basis quantile function is a function of the factors. B-Splines with nonnegative coefficients Boor1974nonnegativeB are used as a primary example for the weight function.
The model is analogous to classic mixture regression model Quandt1972switching where the conditional density function is a mixture of basis density functions, while the parameters of each density function are functions of the factors. However, estimation of mixture regression model relies on computationally expensive nonconvex optimization. Quantile regression KoenkerBassett1978quantile, while calibrated by convex optimization, only estimate conditional quantile at some selected confidence levels and suffers from quantile crossing Koenker2001quantilereg.
Our approach has the following features:
While the quantile functions and weight functions are arbitrary in principle, we primarily focus on quantile functions of common distributions and B-splines with nonnegative coefficients Boor1974nonnegativeB. Based on Papp2011polynomial, Papp2014spline, we prove that the proposed model can approximate any bounded continuous quantile model. The calibration problem is equivalent to minimization of Continuous Ranked Probability Score (CRPS, Hersbach2000CRPS), a popular measure of quality of distributional prediction.
The remainder of this paper is organized as follows. Section (ref) reviews related work and discusses their relations to our proposed approach. Section (ref) formulates the model and gives several illustrative examples. Section (ref) formulates a convex optimization problem of model calibration and discusses its relation to joint quantile regression and CRPS minimization. Section (ref) proves an approximation theorem for the model. Section (ref) proves the asymptotic normality of the estimator. Section (ref) introduces a dimensionality reduction methods by reduced rank tensor and an alternating algorithm for estimation. Section (ref) extends the method to estimate other uncertainty measure based on Risk Quadrangle framework. Section (ref) presents numerical experiments with real-world data.
This section reviews related work and discusses the relation with our approach.
Mixture quantiles use a combination of some basis functions functions to model the quantile function. Various basis functions have been used, such as orthogonal polynomials Sillitto1969, a mixture of normal and Cauchy distributions with linear and quadratic terms Karvanen2006mixtures, a modified logistic distribution Keelin2016metalog, and the quantile functions of the Generalized Beta Distribution of the Second Kind Cheng2022mixquantile. However, these methods lack factor dependence and some do not guarantee that the resulting quantile function is nondecreasing. Our model uses mixture quantiles as the conditional quantile function (Section (ref)). Note that mixture quantiles and mixture densities are two different families of distributions.
Mixture regression shares a similar idea to our approach, where the conditional distribution is modeled as a linear combination of basis functions. The conditional density function of the widely applied Gaussian mixture regression is modeled by a mixture (weighed sum) of Gaussian densities, where the weight, mean and variance of each component can be modeled as functions of factors. Gaussian mixture regression has been studied under different names, such as switching regression model, mixture of experts, mixture of factor analyzers Quandt1972switching, desarbo1988regression, ghahramani1996mixture, VILLANI2009155mixture, Yuksel2012reviewmixture. Similarly, in our model formulation, the weights of quantile functions are functions of the factors (Section (ref)). Mixture regression is generally computationally expensive to calibrate, while our model can be calibrated by convex optimization. The weight function in mixture regression must sum up to one, which is typically achieved by the softmax transformation, while our model only requires nonnegative weight functions.
Quantile regression estimates the conditional quantile of a given confidence level given some factors. The application of spline functions to conduct nonparametric quantile regression has been presented in literature koenker1994spline,He1998bispline,Koenker2011additive. However, conducting multiple quantile regressions separately can result in quantile crossing, particularly in nonlinear models, which impedes the interpretation of the results. To mitigate this problem, various methods have been proposed, such as imposing extra constraints bondell2010noncrossing or rearrangement Chernozhukov2010rearrange. The constraints are often imposed at observed data points or at some confidence levels. We present a novel approach to quantile regression that guarantees noncrossing quantiles without any post-processing. Our method estimates the entire conditional distribution of the response variable, not just the quantiles at several confidence levels, as a result of conducting joint quantile regression. The proper choice of basis functions ensures the validity of the noncrossing property at all points (Section (ref)). The objective function in our approach is a linear combination of pinball losses at different confidence levels (Section (ref)), which has been adopted in various model calibration methods KOENKER1984wqr, Zou2008additive, Sottile2022wqr.
Quantile regression process refers to the regression coefficient as a function of the confidence level of quantile regression. Angrist2006qrpasymptotics shows that the rescaled quantile regression process converges to a zero-mean Gaussian process. A number of independent but closely related studies have explored the use of various basis functions for modeling the quantile regression process, such as Bernstein basis polynomials Reich2011Bqrprocess, P-spline basis LIAN2015qrprocess, common parameterized functions Frumento2016WQR, and monotone B-spline Yuan2017Splineqrprocess. By rearranging the model formulation in Section (ref), it can be observed that our model coincides with the quantile regression process models if we view the individual polynomial terms as factors.
There are several notable differences between our approach and existing research. First, the model calibration in Frumento2016WQR and Reich2011Bqrprocess involves numerical integration and Bayesian method, respectively, while our approach uses convex optimization, which is more efficient. Second, it it can be challenging to ensure noncrossing quantiles when directly modeling the linear quantile regression process. Reich2011Bqrprocess introduces prior latent unconstrained variables; Frumento2016WQR checks the nonnegativity of derivative; Yuan2017Splineqrprocess uses linear constraints when the feasible set is a bounded convex polytope; LIAN2015qrprocess does not guarantee noncrossing quantiles. In contrast, our model guarantees noncrossing conditional quantiles at any two distinct confidence levels and any given factor value. Second, using a spline to model quantile regression process results in bounded conditional quantile, which may lead to severe underestimation of tail risk. Thus it is important to supplement splines with quantile functions of common distributions. Furthermore, all aforementioned quantile regression process models study linear models and can be regarded as special cases of our approach. For linear quantile models to be noncrossing everywhere, the coefficients for each factor must be equal across all confidence levels.
Continuous Ranked Probablity Score (CRPS, Matheson1976score) is a proper scoring rule that is frequently used to measure the quality of a distributional prediction when the response variable is a scalar. Laio2006VerificationTF proposes an equivalent definition that is useful to formulate convex optimization. Hothorn2014Transformation, Gouttes2021NN,BERRISCH2021CRPS use CRPS as the objective function of optimization in learning tasks. Zhang2022CDFmixture proposes a related model where the CDF is a linear combination of basis CDFs. Although not mentioned in the paper, their estimation procedure minimizes CRPS.
Quantile model aggregation aims to aggregate multiple models of conditional quantile function to improve overall performance. The aggregation can be performed across different factor values and different quantile levels. BERRISCH2021CRPS uses B-splines as weight function across confidence levels, while Fakoor2021Aggregation uses neural networks as weight function across different factor values and different quantile levels. Our model can be adapted for quantile model aggregation. On the other hand, Fakoor2021Aggregation considers a regression model called the deep quantile regression by using constant functions as individual models. The conditional quantile function is a neural network, which requires extra efforts to ensure monotonicity.
The linear and polynomial model in chernozhukov2001conditional can be regarded as a special case of our model. Conditional transformation model Hothorn2014Transformation finds the optimal conditional transformation of a fixed quantile function such that it best fits the data. Spline quantile function RNN Gasthaus2019quantile models the parameter of piecewise linear splines with neural networks and calibrate the model by CRPS minimization. Bremnes2020polyNN has a similar model formulation with Bernstein polynomials. Our model is connected to the generalized additive models for location, scale and shape (GAMLSS, Rigby2005gamlss). When there is only one basis quantile function, e.g., quantile function of normal distribution, the mean and variance are modeled as splines of the factor, which is equivalent to generalized additive models for location and scale.
Multivariate B-spline is obtained by tensor product of univariate B-splines. The CP decomposition for spline tensors has been applied in Computer Aided Design (CAD) to reduce model complexity Pan2016cad. For a comprehensive review of tensor decomposition, readers are referred to Kolda2009tensor_intro, Rabanser2017tensor_intro. This study introduces its application to statistical estimation. In high dimensions, the tensor product B-spline has a gigantic number of basis functions. CP decomposition has the appealing interpretation that it reduces the number of basis by assembling them to new ones. Our approach is different from existing work on dimensionality reduction for quantile regression. Lian2019reducedrank proposes reduced matrix rank regression for homoscedastic multiple quantile modeling. Chen2021QFM proposes a different approach to dimensionality reduction for quantile regression determining different principle components for different confidence levels. Our approach results in different factors for different confidence levels, connecting our work with Chen2021QFM (Section (ref)).
The Fundamental Risk Quadrangle quadrangle provides an axiomatic framework to study uncertainty measures and regressions. Besides quantile regression, regressions for other uncertainty measures, such as CVaR (superquantile) regression Rockafellar2013cvarreg, Rockafellar2014cvarreg, Uryasev2019cvarreg and expectile regression Newey1987expectile, have been developed and are studied within the Risk Quadrangle framework rockafellar2018secondorder, kuzmenko2020expectile. Similar to quantile regression, these regressions provide estimation on the uncertainty measure at a given confidence level. By simply replacing the error function in optimization and basis functions in model formulation, our approach can be used to model other uncertainty measures (Section (ref)).
This section describes the Factor Model of Mixture Quantiles and discusses noncrossing quantile and various perspectives to view the model.
The Factor Model of Mixture Quantiles is defined as
where $Q_{0}(p) = 1$. The basis functions $\{Q_i(p)\}_{i=0}^\emph{I}$ are defined on the unit interval, and are linearly independent, i.e., any $Q_i(p)$ is not equal to a linear combination of other $Q_j(p)$, $j\neq i$. The weight function $f_i (\bm{x},\bm{a}_i)$ of each basis function $Q_i(p)$ is a function of the factors $\bm{x}$.
The model formulation (ref) is quite general since $\{Q_i(p)\}_{i=1}^\emph{I}$ and $\{ f_i (\bm{x},\bm{a}_i) \}_{i=0}^\emph{I}$ can be arbitrary functions, providing great flexibility in modeling heteroskedastic data and nonlinear relations. The weight functions determine how the factors impact the scale of basis quantile functions. It has the appealing interpretation that different factors may impact the different part of the distribution. The conditional quantile function has a closed-form expression if all basis quantile functions do. With our model, Monte Carlo simulation can be easily conducted with inverse transform sampling.
We focus on the case where $\{ f_i (\bm{x},\bm{a}_i) \}_{i=0}^\emph{I}$ are splines with basis spline functions $\{B_{ij}(\bm{x})\}_{i,j=0}^{\emph{I},\emph{J}}$. Furthermore, we use the same basis spline functions for all basis quantile functions $\{Q_i(p)\}_{i=1}^\emph{I}$, i.e., $\forall i,j, \;B_{ij}(\bm{x}) = B_{j}(\bm{x})$. Then the model is defined as
where $B_0(x) = 1$, $Q_0(p)=1$. Spline is adaptive to the data and can be optimized with the model in one shot. The linearity with respect to the coefficients not only results in a convenient formulation of convex optimization, but also retains an interpretable factor model structure such that we can analyze the impact of each factor on the shape of the conditional quantile function.
To distinguish $\{Q_i(p)\}_{i=0}^\emph{I}$ and $\{B_j(\bm{x})\}_{j=0}^\emph{J}$, we hereafter refer to $\{Q_i(p)\}_{i=0}^\emph{I}$ as basis quantile functions and $\{B_j(\bm{x})\}_{j=0}^\emph{J}$ as basis spline functions. The meaning of the terms may be clearer if we expand the summation
where $a_{00}$ is the constant location parameter, $ \sum_{j=1}^\emph{J} a_{0 j } B_{j}(\bm{x})$ determines the conditional location, $\sum_{i=1}^\emph{I} a_{i 0 } Q_i(p)$ is the base quantile function that does not vary with factors, $\sum_{i=1}^\emph{I} \sum_{j=1}^\emph{J} a_{i j } B_{j}(\bm{x}) Q_i(p)$ determines the conditional quantile function.
The positions of knots and maximal degree of polynomials of the basis spline functions $\{Q_i(p)\}_{i=1}^\emph{I}$ are chosen or tuned. Splines are known to have poor performance on two ends where there is not much data. To mitigate this problem, we can constrain the function to be linear on both ends. The knot selection problem is often nonconvex. Alternatively, it can be addressed by P-spline Eilers1996Pspline, Eilers2003multiPspline which uses a large number of knots and penalizes the absolute value of the second derivative. The complexity of the model is determined by the number of basis quantile functions, the number of factors and the form of weight functions. The number of parameters in the model is $(\emph{I}+1) \times (\emph{J}+1)$. The degrees of freedom are usually smaller due to constraints in optimization.
We refer to a model as noncrossing if it satisfies that the conditional quantile functions of any two different confidence levels do not cross conditioned on any value of factors. That is,
where $\mathcal{X}$ is the set of all possible values of factors $\bm{x}$. In the literature, the condition is often relaxed to that $ G(p_1,\bm{x},\bm{a}) \leq G(p_2,\bm{x},\bm{a})$ for a certain $\bm{x} \in \mathcal{X}$ or for a certain selected confidence levels $\{p_m\}_{m=1}^\emph{M}$. Even with such simplification, it could require a large number of constraints to guarantee noncrossing.
Instead of relaxing the condition (ref), we propose using a sufficient condition. We show in Section (ref) that such formulation is still highly flexible. Note that the model (ref) satisfies noncrossing condition (ref) if
Nonnegativity constraint is relatively easy to impose in optimization. Although a sufficient condition may seem too restrictive, Theorem (ref) in Section (ref) shows that the model retains good approximation ability.
For $\{Q_i(p)\}_{i=1}^\emph{I}$, we propose two types of nondecreasing functions:
Quantile functions of common distributions are prefered when the shape of the distribution is known to be close to common distributions. When the distribution is multimodal, splines offer greater flexibility. Since splines are bounded by definition, it works well for bounded variable, but needs to be combined with common quantile functions when the distribution has fat tails.
For $\{ B_j (\bm{x}) \}_{j=1}^\emph{J}$, we propose two ways to guarantee nonnegativity:
While the latter is more flexible, it is more computationally expensive to optimize. Besides B-spline, other choices include Bernstein polynomial and M-spline Curry1988Mspline. B-spline and M-spline differ by a constant. Papp2011polynomial, Papp2014spline show that piecewise Bernstein polynomial includes B-spline as a subset. The monotone I-spline used for $\{Q_i(p)\}_{i=1}^\emph{I}$ is obtained by integrating M-spline. The aforementioned choices guarantee nonnegativity in the domain of the spline. When noncrossing property is desired beyond the domain, we can require the spline to have nonnegative partial derivative on the boundary and extrapolate with linear function. The obtained function is nonnegative everywhere. Moreover, one can use neural networks with nonnegative output. This variant will be be considered in the extensions of this paper.
The number of basis functions determines the ability of approximation. Intuitively, suppose we want to approximate only two smooth quantile functions with the proposed model, one quantile function is sufficient.
This subsection provides several examples of the model described in Section (ref). We refer to models from other research to demonstrate the wide applicability of our approach. The connection between the our approach and related work is discussed in detail in Section (ref).
Example (ref), (ref) and (ref) shows that our model incorporates some most common models as special cases. Maximum Likelihood estimation for these models are well studied. We propose a different estimation method in Section (ref). Dynamic models such as Example (ref) and (ref) provide useful complements to classic time series models. Although we focus on linear function in $\bm{a}$ in subsequent sections on model calibration, the function $f(\bm{x},\bm{a})$ in Example (ref) is nonlinear in both $\bm{x}$ and $\bm{a}$. The model can still be calibrated by optimizing our proposed objective function, but the optimization is not convex. Gourieroux2008dynamic also uses $|x|$ and $x^+ = \max\{0,x\}$, $x^- = \max\{0,-x\}$ to replace $x$, so that the factors remain nonnegative. Examples involving spline function can be formulated by replacing factors with splines.
This section first presents the basics of quantile regression KoenkerBassett1978quantile and then introduces the convex formulation of model calibration.
Consider a model of $p$-quantile of response variable $y|\bm{x}$ conditioned on factor $\bm{x}$ with parameter $\bm{a}$
Quantile regression estimates the parameter $\bm{a}$ by minimizing the scaled Koenker-Bassett error of the residuals
where the error is the expected pinball loss
$Z$ is a random variable, $\mathbb{1}_{\{\cdot\}}$ is the indicator function that equals $1$ when the equation in the bracket is true and $0$ if otherwise. Suppose we have $\emph{N}$ pairs of response variables and factors $\{y_i,\bm{x}_i\}_{i=1}^\emph{N}$. The calibration problem is
We write error $\bm{\mathcal{E}}_p$ instead of average of pinball loss $\frac{1}{\emph{N}} \sum_{i=1}^\emph{N} \rho_p(y_i-g(p,\bm{x}_i,\bm{a}))$, since some errors cannot be written as an average in the generalization in Section (ref).
The optimization problem statement for finding optimal parameters $\bm{a}$ is formulated as follows
$w(p)$ can be chosen to focus on the distribution tail or body. Inherent from the pinball loss, the calibration is robust to outliers. While we consider the case where $f_{i}(\bm{x},\bm{a})$ are spline functions, the calibration method works for the general case (ref).
Next, we discretize the problem (ref) by using a grid on $p$. The resultant optimization problem is still a convex programming problem.
Similar to quantile regression, the problem can be reduced to linear programming. Note that although we select only a finite number of confidence levels in the discrete variant, the procedure still estimates the whole conditional quantile function, and the quantiles are noncrossing. The calibration problem also allows one to focus on the tail of the distribution by assigning higher weights on errors with tail confidence levels.
Define $\emph{M}$ discrete random variables $\widehat{L}_m(\bm{a};\bm{Y}^\emph{N},\bm{x}^\emph{N})$, each taking with equal probabilities the values $$ y_n - \sum_{j=0}^\emph{J} \lambda_{ jm } B_{j}(\bm{x}_n) , \quad n=1,\ldots,\emph{N} \;. $$ We can write (ref) equivalently as
This problem formulation makes it clear that the calibration in Problem (ref) is equivalent to conducting several quantile regressions in one shot with constraints on the parameters. Since $\sum_{j=0}^\emph{J} \lambda_{ jm } B_{j}(\bm{x}_n) $ is a spline of scalar factor $x_j$, the objective function in (ref) is the sum of objective functions of spline quantile regressions with confidence levels $p_m$, $m=1,\cdots,\emph{M}$. $\{\lambda_{jm}\}_{j=0,m=1}^{\emph{J},\emph{M}}$ are associated by systems of equations (ref), where the coefficients $\bm{a}$ are constrained by feasible set $\mathcal{A}$.
As $\emph{J}$ systems of equations, (ref) can be written in matrix format
where $\bm{Q} = [Q_i(p_m)]_{\emph{M} \times (\emph{I}+1)}', \bm{A} = [a_{ij}]_{(I+1) \times (\emph{J}+1)}, \bm{\Lambda} = [\lambda_{jm}]_{\emph{M} \times (\emph{J}+1)}'$.
If the solution to the unconstrained problem (ref) is in the feasible set defined by (ref)(ref), we can solve Problem (ref) with the following two steps.
In both steps, the solution can be nonunique. The uniqueness of solution to quantile regression is discussed in koenker2005book. Portnoy1991qrbreakpoints shows that the number of distinct solutions to quantile regression when the confidence level $p$ varies in $(0,1)$ is $\mathcal{O}(N\log N)$ in probability, which grows slower than the upper bound $\binom{N}{M \times (\emph{J}+1)}$. The uniqueness of solution to the systems of equations depends on $\bm{Q}$ when $\mathcal{A} = \mathrm{R}^{(I+1) \times (\emph{J}+1)}$. For proper choice of independent basis quantile functions $\{Q_i(p)\}_{i=1}^\emph{I}$, $\bm{Q}$ has full column rank when $\emph{M}\geq \emph{I}+1$. To check feasibility when $\mathcal{A} = \mathrm{R}_{\geq 0}^{(I+1) \times (\emph{J}+1)}$, we can solve the linear programming problem $\min_{\bm{a} \geq 0} || \bm{Q} \bm{A} - \bm{\Lambda} ||_1$.
Continuous Ranking Probability Score (CRPS) is frequently used to evaluate the quality of a distributional forecast. One may want to directly optimize the measure by which the model is evaluated in model calibration. We show that Problem (ref) is equivalent to CRPS minimization.
For a CDF $F$ and an observation $y_i$ of the response variable $y$, CRPS is defined by
This is the squared distance between $F$ and the CDF of a single observation $y_i$ of the response variable $y$. For the corresponding quantile function $Q$, i.e., the generalized inverse function of $F$, and the response variable $y_i$, CRPS has the following equivalent definition
Consider calibrating the model $G(p,\bm{\theta},\bm{a})$ by minimizing the sum of CRPS of the response variables
We see that (ref) is equal to the objective function in Problem (ref) with uniform weight $w(p)$ by exchanging the integral and sum
When there is no basis quantile function, the CRPS minimization reduces the model to the least absolute deviation regression splines.
Our model can approximate any bounded conditional quantile model to arbitrary precision as the number of knots tends to infinity. Since the considered model has nonnegative parameters, classic approximation theorem of polynomials cannot be directly applied.
Theorem (ref) is an application of Papp2011polynomial, Papp2014spline.
This section contains the asymptotic properties of the estimator.
Theorem (ref), (ref) in the following from Frumento2016WQR are applications of Newey2.
This section introduces the basics of CP decomposition of tensors, and an alternating algorithm for dimensionality reduction. Spline-based methods are often deemed inadequate for high-dimensional problems due to the exponential growth of parameters in multivariate splines, which poses challenges in both statistical inference and optimization. We propose an dimensionality reduction method, with linear growth rate. Additionally, our proposed alternating algorithm solves a smaller convex optimization problem in each iteration. It differs from the classic Alternating Least-Squares Algorithm through the minimization of a distinct objective function and the addition of a nonnegative constraint in each step, given that our model incorporates nonnegative parameters.
A tensor is a multidimensional array. CP decomposition for a rank-$\emph{R}$ tensor $\bm{A}$ is defined as
where $\{\bm{u}_r^k\}_{r=1,k=0}^{\emph{R},\emph{K}}$ are vectors.
Let $\bm{U}_k = (\bm{u}_1^{(k)},\cdots,\bm{u}_\emph{R}^{(k)})$, $k=1,\cdots,\emph{K}$. CP decomposition is often estimated by Alternating Least-squares Algorithm, which estimates $\bm{U}_k$ one by one with fixed $\{\bm{U}_{k'}\}_{k'\neq k}$. Each step is a convex optimization problem.
The following equation is useful for subsequent interpretation of reduced rank method
i.e., $ \bm{u}_r^{(0)} \otimes \cdots \otimes \bm{u}_r^{(\emph{K})} \cdot \bm{v}_r^{(0)} \otimes \cdots \otimes \bm{v}_r^{(\emph{K})} = \left( \bm{u}_r^{(0)} \cdot \bm{v}_r^{(0)} \right) \cdots \left( \bm{u}_r^{(\emph{K})} \cdot \bm{v}_r^{(\emph{K})} \right) $.
For B-spline, $\emph{L}=$ degree of univariate polynomial$+$number of knots$+1$. Multivariate B-spline basis is obtained by tensor product of univariate B-spline basis.
With the above notations, we can write the model (ref) in tensor format
Note that in our model formulation, $\bm{A}\geq 0$. Consider a nonnegative version of CP decomposition (ref) on the parameter tensor $\bm{A}$ where $\bm{u}_r^{(k)} \geq 0$. With (ref)(ref), we have
Thus reduced rank method has a straightforward interpretation. The model is reduced to a linear combination of $\emph{R}$ functions from $(\emph{I}+1)(\emph{J}+1)$. In the extreme case where $\emph{R} = 1$, the model reduced to a heteroscedastic model where the conditional quantile function is $\bm{u}_r^{(0)} \cdot \bm{Q}$, whose scale is controlled by a scalar function $\prod_{k=1}^\emph{K} \bm{u}_r^{(k)} \cdot \bm{B}_k$. It can be regarded as obtaining $\emph{R}$ new basis quantile functions $\{\bm{u}_r^{(0)} \cdot \bm{Q}\}_{r=1}^\emph{R}$, likewise for basis spline functions. The new basis, having a smaller number, can still have undesirable smoothness condition. We expect better performance when it is used along with penalties in P-spline. The reduced basis functions are a linear combination of all original basis functions, while sparse optimization leads to a linear combination of a subset by forcing zeros among the parameters. We find that the solution to Problem (ref) is often sparse with many small nonzero values. Thus the low rank decomposition is expected to produce good approximation, although the decomposition (ref) is not always valid for any $\bm{A}$ for a low rank $\emph{R}$.
Nonnegativity constraint is imposed in each step of the algorithm. However, in the model formulation, there is no sign constriant on the coefficients of the splines that determine the conditional location. Thus we add a constant to all the observations so that the conditional location also satisfies the constraint.
We propose an alternating algorithm to find $\bm{A}$. Define $\widetilde{L}\left( \{\bm{U}^{(k)}\}_{k=0}^\emph{K};\bm{Y}^\emph{N},\bm{x}^\emph{N},p \right) =$ discrete residual random variable taking with equal probabilities the following values $$y_n - \sum_{r=1}^\emph{R} \left( \bm{u}_r^{(0)} \cdot \bm{Q}(p) \right) \left( \bm{u}_r^{(1)} \cdot \bm{B}_1(\bm{x}_n) \right) \cdots \left( \bm{u}_r^{(\emph{K})} \cdot \bm{B}_\emph{K}(\bm{x}_n) \right) \;, \quad n=1,\ldots,\emph{N} \;. $$
Other stopping rules can be used as well, such as the number of iterations, the difference between consecutively updated objective values. Each step in the algorithm is a convex optimization problem. The objective function is nonincreasing in each step. The algorithm will find a local minimum at the end. Different initialization is needed for finding a smaller local minimum.
This section first introduces CVaR and the Risk Quadrangle framework. Then the framework is used to extend our approach to estimation of conditional functions of other uncertainty measures such as CVaR.
\paragraph{Conditional Value-at-Risk (CVaR)}
For a continuous random variable $Z$, the $p$-CVaR is defined as the average value that exceeds the $p$-quantile of $Z$
As an uncertainty measure, it takes into account not only the probability, but also the scale of extreme losses.
\paragraph{Error and Statistic} The Risk Quadrangle is a general framework describing the relation between risk, deviation, regret, error and statistics. This section presents the basics of the regression theory in Risk Quadrangle framework necessary for understanding the subsequent generalization of the model.
A functional of random variable $Z$ is called a regular measure of error if it satisfies the following conditions: it has values in $[0,\infty)$ ; it is closed convex in $Z$ with $ \bm{\mathcal{E}}(0)=0,\; \bm{\mathcal{E}}(Z) > 0 \; \text{when} \; Z \neq 0 \;; $ for sequences of random variables $\{Z_k\}_{k=1}^\infty$, if $ \lim_{k \rightarrow \infty} \bm{\mathcal{E}} (Z_k)=0 $, then $\lim_{k\rightarrow 0}E[Z_k] = 0 \;. $
The statistic associated with $Z$ by $\bm{\mathcal{E}}$ is defined by
Two prominent examples are the quantile-based quadrangle and the superquantile-based Quadrangle.
Quantile-based quadrangle (at any confidence level $p\in (0,1)$):
Superquantile-based quadrangle (at any confidence level $p\in (0,1)$):
\paragraph{Regression in Risk Quadrangle Framework} We can conduct the following regression to estimate the statistic $S$ associated with $Y$ by $\bm{\mathcal{E}}_p$
for given random variables $\bm{x}, Y$, and some given class $\mathcal{C}$ of functions $g$. With $g$ obtained in (ref), we can estimate the statistic of $Y$ conditioned on observing $\bm{x}$ by plugging in values of $\bm{x}$ to $g$.
By replacing the error $\bm{\mathcal{E}}$ in Problem (ref) with other errors pamameterized by a confidence level, our approach can be generalized to model other functions of uncertainty measures. Examples of Risk Quadrangle parameterized by a parameter in $(0,1)$ include quantile-based, CVaR-based, and expectile-based quadrangle.
An example of Factor Model of Mixture CVaRs is given as follows.
This section presents two numerical experiments with real-world data.
The first experiment focuses on the Factor Model of Mixture Quantiles (Section (ref)) and its calibration by convex optimization (Section (ref)). The model is used to estimate the conditional distribution of Durable Goods as a function of Nondurable Goods and Services.
\paragraph{Data} The response variable is the percentage change from previous year of the Consumer Price Index for All Urban Consumers: Durables in US City Average, while the factors are the percentage change from previous year of the Consumer Price Index for All Urban Consumers: Nondurables in US City Average and the percentage change from previous year of the Consumer Price Index for All Urban Consumers: Services in US City Average. For convenience, we denote them by Durable Goods, Nondurable Goods and Services in subsequent sections. The monthly data, covering the period from January 1957 to January 2023, is described in mccracken2016data and downloaded from Fredcpi. The data are standardized to have zero median and unit interquartile range. The data is randomly shuffled before experiments.
\paragraph{Model} Let
Consider the Two-Factor Model of Mixture Quantiles
For each bivariate B-spline $\{f_i\}_{i=0,1,2,3}$, we use basis of degree three with six knots. The interpretation is that $f_1(x_1,x_2)$ determines the body of the distribution, while $f_2(x_1,x_2)$ and $f_3(x_1,x_2)$ determine the left and right tails, respectively.
\paragraph{Calibration} The confidence levels used in the optimization are $0.01,0.05,0.15,0.25,\cdots,0.95,0.99$. The weights for the confidence levels are $20,10,1,,:\cdots,1,10,20$ and normalized such that they sum up to $1$. We choose $10$ equidistant knots for each factor. We add a penalty term on the squared difference between adjacent coefficients Eilers1996Pspline, Eilers2003multiPspline for each spline function $\{f_i\}_{i=0,1,2,3}$ to reduce overfit. The penalty coefficients are fixed at $0.01$ and not optimized. The model calibration is conducted with R package Portfolio Safeguard zabarankin2016book.
\paragraph{Benchmarks} We choose the Gaussian mixture regression and the generalized additive model (GAM) as the benchmarks. For the Gaussian mixture regression, we use BIC as the model selection criteria, which is one of the built-in methods in the package Mixtools. The weights of the components are constants in the model. For the generalized additive model, the mean, standard deviation, shape and skewness of the Skew $t$ type 2 distribution Azzalini2003student are functions of the factors. The additive formula of the functions are sum of univariate cubic splines of each factor with an additional linear interaction term. The univariate cubic splines are penalized on the second derivative of the functions. The benchmark models are implemented by R packages Mixtools Benaglia2009mixtools and GAMLSS Rigby2005gamlss, respectively.
\paragraph{Performance Measure} We conduct 10-fold cross-validation and use the out-of-sample CRPS and coverage rate as performance measures. The CRPS is calculated by sum of CRPS of distributional prediction on response variable. The coverage rate is calculated by the percentage of the response variable that falls within the predicted conditional interval of two specified quantiles. A small CRPS and a coverage rate close to the target indicate a high-quality distributional prediction.
\paragraph{Results} Results are presented in Table (ref) and Table (ref). In general, our approach has comparable performance with the benchmarks. Our approach exhibits improved average out-of-sample coverage rate. The average out-of-sample CRPS of our approach is lower than the Gaussian mixture regression but higher than the generalized additive model. However, the generalized additive model often yield poor extrapolation results when the test data lies outside the domain of training data, resulting in nonexistent second moments and therefore nonexistent CRPS values.
We can obtain the $p$-quantile surface of the response variable by varying the values of $(x_1,x_2)$ with a fixed $p$. We use the model calibrated in the $10$th fold of the cross-validation as an example. Figure (ref) shows three selected quantile surfaces for visualization, demonstrating the nonlinear relations incorporated in the model. The surfaces do not to cross each other.
The second numerical experiment examines the reduced rank method (Section (ref)) for the Factor Model of Mixture CVaRs (Section (ref)). We illustrate that the low-rank representation of parameter tensor is efficient to compute by showing the in-sample error curves. The full data set described in Section (ref) is used in the calibration.
\paragraph{Model} Consider the Factor Model of Mixture CVaRs in tensor format
where $\Bar{\bm{Q}}(p)$ is the vector of basis CVaR functions. We use $I$-spline basis of degree three with five knots and CVaR function of exponential distribution as basis CVaR functions. for splines $\bm{B}_k(x_k)$, we still use B-spline of degree three with six knots.
\paragraph{Calibration} The experiment considers cases where the rank, $\emph{R}$ in (ref), equals 1, 3, and 10. For each case, 21 optimization steps are conducted using the alternating algorithm. The error function for the regression is defined by replacing $\bm{\mathcal{E}}_p$ in Problem (ref) with the error defined in (ref). Since CVaR is often used to measure tail risk, we choose more tail confidence levels in the optimization. The confidence levels are $0.25,0.55,0.60,0.65,\cdots,0.95$. No penalty is applied to the smoothness of splines as in Section (ref), since the aim is to compare the models with varying ranks. Note that the low-rank representation alone does not guarantee improvements in out-of-sample test, since the smoothness condition is unconstrained and may lead to overfit. The combination of the low-rank representation with P-spline is left for future study.
\paragraph{Results} The objective values at each step are displayed in Figure (ref). A higher rank corresponds to a smaller final objective value and a greater number of steps to converge. Nonetheless, the optimization converges in a few steps for all three cases, demonstrating the efficiency of the algorithm for estimation.