The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
64,192 characters
Factor Model of Mixtures
\maketitle
\begin{abstract}
This paper proposes a new approach to
estimating the distribution of a response variable conditioned on observing some factors. The proposed approach possesses desirable properties of flexibility, interpretability, tractability and extendability.
The conditional quantile function is modeled by a mixture (weighted sum) of basis quantile functions, with the weights depending on factors.
The calibration problem is formulated as a convex optimization problem. It can be viewed as conducting quantile regressions for all confidence levels simultaneously while avoiding quantile crossing by definition. The calibration problem is equivalent to minimizing the continuous ranked probability score (CRPS).
Based on the canonical polyadic (CP) decomposition of tensors, we propose a dimensionality reduction method that reduces the rank of the parameter tensor and propose an alternating algorithm for estimation.
Additionally, based on Risk Quadrangle framework, we generalize the approach to conditional distributions defined by Conditional Value-at-Risk (CVaR), expectile and other functions of uncertainty measures. Although this paper focuses on using splines as the weight functions, it can be extended to neural networks.
Numerical experiments demonstrate the effectiveness of our approach.
\end{abstract}
\section{Introduction}
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 \citep{Boor1974nonnegativeB} are used as a primary example for the weight function.
The model is analogous to classic mixture regression model \citep{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 \citep{KoenkerBassett1978quantile}, while calibrated by convex optimization, only estimate conditional quantile at some selected confidence levels and suffers from quantile crossing \citep{Koenker2001quantilereg}.
Our approach has the following features:
\begin{itemize}
\item \textit{Flexibility} --
The conditional quantile function is flexible in shape to accommodate fat tails and multimodality. The weight function can capture nonlinear relations between factors and quantiles of the response variable (Section \ref{sec_model}). The model can approximate any bounded conditional quantile model in the limit (Section \ref{sec_approximation}).
\item \textit{Interpretability} -- The model can be viewed as a factor model such that the impact of factors on the shape of conditional quantile function can be traced analytically. Both the conditional quantile function and the quantile (hyper)surface have closed-form expressions (Section \ref{sec_model}). Furthermore, the quantile (hyper)surfaces do not cross by definition.
\item \textit{Tractability} --
The model calibration is formulated as a linear regression problem similar to quantile regression, which can be efficiently solved by convex and linear programming (Section \ref{sec_calibration}). Inspired by Canonical Polyadic (CP) decomposition of tensors \citep{Hitchcock1927cp,
Harshman1970cp,Carroll1970cp}, we also propose a dimensionality reduction method by reduced rank tensor and an alternating algorithm for estimation of large models (Section \ref{sec_dim_reduction}).
\item \textit{Extendability} --
The weight function can be modeled by functions beyond splines such as neural networks. Furthermore, we generalize the approach to modeling other functions of uncertainty measures beyond quantile, such as Conditional Value-at-Risk (CVaR, \cite{cvar1,cvar2}) and expectile function \cite{Newey1987expectile} based on the Risk Quadrangle framework \citep{quadrangle} (Section \ref{sec_quadrangle}).
\end{itemize}
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 \citep{Boor1974nonnegativeB}.
Based on \citep{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, \cite{Hersbach2000CRPS}), a popular measure of quality of distributional prediction.
The remainder of this paper is organized as follows. Section \ref{sec_review} reviews related work and discusses their relations to our proposed approach. Section \ref{sec_model} formulates the model and gives several illustrative examples. Section \ref{sec_calibration} formulates a convex optimization problem of model calibration
and discusses its relation to joint quantile regression and CRPS minimization.
Section \ref{sec_approximation} proves an approximation theorem for the model. Section \ref{sec_asym}
proves the asymptotic normality of the estimator.
Section \ref{sec_dim_reduction} introduces a dimensionality reduction methods by reduced rank tensor and an alternating algorithm for estimation.
Section \ref{sec_quadrangle} extends the method to estimate other uncertainty measure based on Risk Quadrangle framework. Section \ref{sec_numerical} presents numerical experiments with real-world data.
\section{Discussion on Related Work}\label{sec_review}
This section reviews related work and discusses the relation with our approach.
\textit{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 \citep{Sillitto1969}, a mixture of normal and Cauchy distributions with linear and quadratic terms \citep{Karvanen2006mixtures}, a modified logistic distribution \citep{Keelin2016metalog}, and the quantile functions of the Generalized Beta Distribution of the Second Kind \citep{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{sec_model}). Note that mixture quantiles and mixture densities are two different families of distributions.
\textit{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 \citep{Quandt1972switching,
desarbo1988regression,
ghahramani1996mixture,
VILLANI2009155mixture,
Yuksel2012reviewmixture}.
Similarly, in our model formulation, the weights of quantile functions are functions of the factors (Section \ref{sec_model}). 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.
\textit{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 \citep{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 \citep{bondell2010noncrossing} or rearrangement \citep{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{sec_model}). The objective function in our approach is a linear combination of pinball losses at different confidence levels (Section \ref{sec_calibration}), which has been adopted in various model calibration methods \citep{KOENKER1984wqr,
Zou2008additive,
Sottile2022wqr}.
\textit{Quantile regression process} refers to the regression coefficient as a function of the confidence level of quantile regression. \cite{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 \citep{Reich2011Bqrprocess}, P-spline basis \citep{LIAN2015qrprocess}, common parameterized functions \citep{Frumento2016WQR}, and monotone B-spline \citep{Yuan2017Splineqrprocess}.
By rearranging the model formulation in Section \ref{sec_model},
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 \cite{Frumento2016WQR} and \cite{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. \cite{Reich2011Bqrprocess} introduces prior latent unconstrained variables;
\cite{Frumento2016WQR} checks the nonnegativity of derivative;
\cite{Yuan2017Splineqrprocess} uses linear constraints when the feasible set is a bounded convex polytope; \cite{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.
\textit{Continuous Ranked Probablity Score} (CRPS, \cite{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.
\cite{Laio2006VerificationTF} proposes an equivalent definition that is useful to formulate convex optimization. \cite{Hothorn2014Transformation,
Gouttes2021NN,BERRISCH2021CRPS} use CRPS as the objective function of optimization in learning tasks.
\cite{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.
\textit{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. \cite{BERRISCH2021CRPS} uses B-splines as weight function across confidence levels, while \cite{
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, \cite{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.
\textit{The linear and polynomial model} in \cite{chernozhukov2001conditional} can be regarded as a special case of our model.
\textit{Conditional transformation model} \citep{Hothorn2014Transformation} finds the optimal conditional transformation of a fixed quantile function such that it best fits the data.
\textit{Spline quantile function RNN} \citep{Gasthaus2019quantile} models the parameter of piecewise linear splines with neural networks and calibrate the model by CRPS minimization. \cite{Bremnes2020polyNN} has a similar model formulation with Bernstein polynomials. Our model is connected to the \textit{generalized additive models for location, scale and shape} (GAMLSS, \cite{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 \citep{Pan2016cad}. For a comprehensive review of tensor decomposition, readers are referred to \cite{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.
\cite{Lian2019reducedrank} proposes reduced matrix rank regression for homoscedastic multiple quantile modeling. \cite{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 \cite{Chen2021QFM} (Section \ref{sec_dim_reduction}).
The Fundamental Risk Quadrangle \citep{quadrangle} provides an axiomatic framework to study uncertainty measures and regressions.
Besides quantile regression, regressions for other uncertainty measures, such as CVaR (superquantile) regression \citep{Rockafellar2013cvarreg, Rockafellar2014cvarreg, Uryasev2019cvarreg} and expectile regression \citep{Newey1987expectile}, have been developed and are studied within the Risk Quadrangle framework \citep{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{sec_quadrangle}).
\section{Model Description}\label{sec_model}
This section describes the Factor Model of Mixture Quantiles and discusses noncrossing quantile and various perspectives to view the model.
\noindent
\begin{itemize}
\item $p = $ confidence level of a quantile function
\item $\bm{x} = $ vector of factors
\item $\bm{a}=$ vector of parameters (unknown coefficients to be estimated)
\item $Q(p)=$ quantile function
\item $G(p,\bm{x},\bm{a})=$ model with parameters $\bm{a}$ that outputs the $p$-quantile of the response variable conditioned on observing factor $\bm{x}$
\item $\emph{I} = $ number of basis functions in the mixture
\item $\emph{J} =$ number of basis functions
\end{itemize}
\subsection{Factor Model of Mixture Quantiles}
The Factor Model of Mixture Quantiles is defined as
\begin{equation} \label{model}
G(p,\bm{x},\bm{a})
=
\sum_{i=0}^\emph{I} f_i (\bm{x},\bm{a}_i) Q_i(p) \;,
\end{equation}
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 \eqref{model} 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
\begin{align}\label{model_2}
G(p,\bm{x},\bm{a})
=
\sum_{i=0}^\emph{I} \sum_{j=0}^\emph{J}
a_{i j }
B_{j}(\bm{x})
Q_i(p)
\;,
\end{align}
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
\begin{align}
G(p,\bm{x},\bm{a})
=
a_{00}
+
\sum_{j=1}^\emph{J}
a_{0 j }
B_{j}(\bm{x})
+
\sum_{i=1}^\emph{I}
a_{i 0 }
Q_i(p)
+
\sum_{i=1}^\emph{I} \sum_{j=1}^\emph{J}
a_{i j }
B_{j}(\bm{x})
Q_i(p) \;,
\end{align}
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 \citep{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.
\subsection{Noncrossing Quantile Model}\label{sec_noncrossing}
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,
\begin{align}\label{noncrossing}
\forall \bm{x} \in \mathcal{X},p_1,p_2,\; 0< p_1 < p_2<1,\; G(p_1,\bm{x},\bm{a}) \leq G(p_2,\bm{x},\bm{a}) \;,
\end{align}
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 \eqref{noncrossing}, we propose using a sufficient condition. We show in Section \ref{sec_approximation} that such formulation is still highly flexible. Note that the model \eqref{model_2} satisfies noncrossing condition \eqref{noncrossing} if
\begin{itemize}
\item $\{Q_i(p)\}_{i=1}^\emph{I}$ are nondecreasing;
\item $\{ B_j (\bm{x}) \}_{j=1}^\emph{J}$ are nonnegative;
\item $\{a_{ij}\}_{i\neq 0}$ are nonnegative.
\end{itemize}
Nonnegativity constraint is relatively easy to impose in optimization. Although a sufficient condition may seem too restrictive, Theorem \ref{theorem_approx} in Section \ref{sec_approximation} shows that the model retains good approximation ability.
For $\{Q_i(p)\}_{i=1}^\emph{I}$, we propose two types of nondecreasing functions:
\begin{itemize}
\item quantile functions of common distributions;
\item monotone basis spline functions such as I-spline \citep{Ramsay1988Ispline}.
\end{itemize}
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:
\begin{itemize}
\item nonnegative basis spline functions such as B-spline with nonnegative coefficients \citep{Boor1974nonnegativeB,
Papp2011polynomial, Papp2014spline};
\item arbitrary basis spline functions with constraints that characterize nonnegative polynomials in the optimization \citep{
Papp2011polynomial, Papp2014spline}.
\end{itemize}
While the latter is more flexible, it is more computationally expensive to optimize.
Besides B-spline, other choices include Bernstein polynomial and M-spline \citep{Curry1988Mspline}. B-spline and M-spline differ by a constant. \cite{Papp2011polynomial}, \cite{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.
\subsection{Examples}
This subsection provides several examples of the model described in Section \ref{sec_model}. 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{sec_review}.
\begin{example}\label{eg_locationscale}
Location-scale model of normal distribution
\begin{align}
G(p,\bm{x},\bm{a})
=
a_{00} + \sqrt{2}{\mathop{{\rm erf}}}^{-1}(2p-1) a_{10} \;,
\end{align}
where $\sqrt{2}{\mathop{{\rm erf}}}^{-1}(2p-1)$ is the quantile function of standard normal distribution.
\end{example}
\begin{example} \label{eg_mixture}
Logistic-normal mixture of quantiles without factor dependence \citep{Karvanen2006mixtures,
Keelin2016metalog,
Cheng2022mixquantile}
\begin{equation}
G(p,\bm{x},\bm{a})
=
a_{00} + \log(\frac{p}{1-p}) a_{10}
+ \sqrt{2}{\mathop{{\rm erf}}}^{-1}(2p-1) a_{20} \;,
\end{equation}
where $\log(\frac{p}{1-p})$ is the quantile function of standard logistic distribution.
\end{example}
\begin{example}\label{eg_linearnormal}
Linear model with normal noise
\begin{equation}
G(p,\bm{x},\bm{a}) = a_{00} + x a_{01} + \sqrt{2}{\mathop{{\rm erf}}}^{-1}(2p-1) \;.
\end{equation}
\end{example}
\begin{example}\label{eg_qrproess}
Quantile regression process model \citep{Reich2011Bqrprocess,
LIAN2015qrprocess,
Frumento2016WQR,
Yuan2017Splineqrprocess}
\begin{align}
G(p,\bm{x},\bm{a})
=
\sum_{k_1=1}^{K} B_{k_1}(p) x_1 a_{1 k_1} +
\sum_{k_2=1}^{K} B_{k_2}(p) x_2 a_{2 k_2} \;,
\end{align}
where $B_{k}(p)$ can be Bernstein basis polynomials or B-spline basis.
\end{example}
\begin{example}\label{eg_arch}
Autoregressive conditional heteroskedasticity-1 (ARCH(1)) model with normal noise \citep{Engle1982arch}
\begin{equation}
G_t(p,\bm{x},\bm{a})
=
\sqrt{2}{\mathop{{\rm erf}}}^{-1}(2p-1) \cdot \sqrt{ a_{10} + x_{t-1}^2 a_{11} } \;,
\end{equation}
where $G_t$ is the parameterized quantile function of $x_t$ at time $t$, $x_{t-1}$ is the response variable at at time $t-1$. At every time $t$, the quantile function $G_t(p,\bm{x},\bm{a})$ depends on the previous response variable $x_{t-1}$.
\end{example}
\begin{example} \label{eg_dynamic}
Dynamic quantile model \citep{
Gourieroux2008dynamic}
\begin{equation}
G_t(p,\bm{x},\bm{a})
=
a_{00} + x_{t-1} a_{01} + \log(\frac{p}{1-p}) \cdot ( a_{10} + x_{t-1}^2 a_{11}) \;.
\end{equation}
where $G_t$ is the parameterized quantile function of $x_t$ at time $t$, $x_{t-1}$ is the response variable at at time $t-1$. At every time $t$, the quantile function $G_t(p,\bm{x},\bm{a})$ depends on past response variable $x_{t-1}$.
The response variable at time $t$ is generated by plugging a sample $p$ from uniform distribution $U(0,1)$ and the past response variable $x_{t-1}$ to $G_t(p,\bm{x},\bm{a})$.
\end{example}
Example \ref{eg_locationscale}, \ref{eg_linearnormal} and \ref{eg_arch} 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{sec_calibration}. Dynamic models such as Example \ref{eg_arch} and \ref{eg_dynamic} 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{eg_arch} 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.
\cite{
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.
\section{Model Calibration by Convex Optimization}
\label{sec_calibration}
This section first presents the basics of quantile regression \citep{KoenkerBassett1978quantile} and then introduces the convex formulation of model calibration.
\subsection{Quantile Regression}\label{sec_quantile_reg}
Consider a model of $p$-quantile of response variable $y|\bm{x}$ conditioned on factor $\bm{x}$ with parameter $\bm{a}$
\begin{align}
q_p(y|\bm{x}) = g(p,\bm{x},\bm{a}) \;.
\end{align}
Quantile regression estimates the parameter $\bm{a}$ by minimizing the scaled Koenker-Bassett error of the residuals
\begin{align}
\min_{\bm{a}} \bm{\mathcal{E}}_p(y-G(p,\bm{x},\bm{a}))
\end{align}
where the error is the expected pinball loss
\begin{align}
\bm{\mathcal{E}}_p(Z) = E[\rho_p(Z)] \;, \quad
\rho_p(Z) = pZ \mathbb{1}_{\{ Z > 0 \}} - (1-p) Z \mathbb{1}_{\{ Z \leq 0 \}} \;, \label{def_pinball}
\end{align}
$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
\begin{align}
\min_{\bm{a}} \bm{\mathcal{E}}_p (y_i-g(p,\bm{x}_i,\bm{a})) \;.
\end{align}
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{sec_quadrangle}.
\subsection{Calibration: Optimization Problem Statement}
\begin{itemize}
\item $\emph{M}=$ number of grid points of discretized confidence level
\item $\emph{N}=$ sample size
\item $\bm{Y}^\emph{N} = (y_1,\cdots,y_\emph{N})=$ response variables of dependent variable
\item $\bm{x}_i=$ vector of factors corresponding to response variable $y_n$, $n=1,\cdots,\emph{N}$;
\item $\bm{x}^\emph{N}=(\bm{x}_1,\cdots,\bm{x}_\emph{N})= $ set of factor vectors
\item $\bm{a}=$ vector of parameters
\item $\mathcal{A}=$ feasible set of $\bm{a}$ determined by constraints
\item $L(\bm{a};\bm{Y}^\emph{N},\bm{x}^\emph{N},p)=$ discrete residual random variable taking with equal probabilities the following values
$$y_n -
\sum_{i=0}^\emph{I} \sum_{j=0}^\emph{J}
a_{i j }
B_{j}(\bm{x}_n)
Q_i(p)
\;,
\quad n=1,\ldots,\emph{N}
$$
\end{itemize}
\bigskip
The optimization problem statement for finding optimal parameters $\bm{a}$ is formulated as follows
\begin{problem}\label{problem_1}
\begin{equation}
\label{continuous}
\begin{aligned}
\min_{\bm{a}} \; & \int_0^1 w(p) \bm{\mathcal{E}}_{p} \left( L(\bm{a};\bm{Y}^\emph{N},\bm{x}^\emph{N},p) \right) \mathrm{d}p \\
& \text{subject to} \; \bm{a}\in\mathcal{A}
\;,
\end{aligned}
\end{equation}
where $w(p)$ is a nonnegative weight function satisfying $\int_0^1 w(p) \mathrm{d}p= 1$.
\end{problem}
$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 \eqref{model}.
Next, we discretize the problem \eqref{continuous} by using a grid on $p$. The resultant optimization problem is still a convex programming problem.
\begin{problem}
\label{first_method}
Discrete variant:
\begin{equation}\label{problem_statement}
\begin{aligned}
\min_{\bm{a}} \; & \sum_{m=1}^\emph{M} w_m \bm{\mathcal{E}}_{p_m} \left( L(\bm{a};\bm{Y}^\emph{N},\bm{x}^\emph{N},p_m) \right) \\
& \text{subject to} \; \bm{a}\in\mathcal{A}
\;,
\end{aligned}
\end{equation}
where $w_m$ are nonnegative weights satisfying $\sum_{m=1}^\emph{M} w_m = 1$.
\end{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.
\subsection{Equivalence to Constrained Joint Quantile Regression}\label{two_step_procedure}
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 \eqref{problem_statement} equivalently as
\begin{problem} \label{equivalent_formulation_0}
\begin{align}
\min_{\bm{a},\, \bm{\lambda} }\; &
\sum_{m=1}^\emph{M} w_m \bm{\mathcal{E}}_{p_m} \left( \widehat{L}_m(\bm{a};\bm{Y}^\emph{N},\bm{x}^\emph{N})\right) \label{objective_0}
\\
\text{subject to} \; & \bm{a}\in\mathcal{A} \label{constraint_set_0}
\\
& \sum_{i=0}^\emph{I}
Q_i(p_m) a_{i}
= \lambda_{jm} \;, \label{sys_equations_0}
\\
& m = 1,\ldots,\emph{M},\; j=0,\cdots,\emph{J}
\;. \nonumber
\end{align}
\end{problem}
This problem formulation makes it clear that the calibration in Problem \ref{first_method} 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 \eqref{objective_0} 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 \eqref{sys_equations_0}, where the coefficients $\bm{a}$ are constrained by feasible set $\mathcal{A}$.
As $\emph{J}$ systems of equations, \eqref{sys_equations_0} can be written in matrix format
\begin{align}\label{matrix_equations_0}
\bm{Q} \bm{A} = \bm{\Lambda} \;.
\end{align}
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 \eqref{objective_0} is in the feasible set defined by \eqref{constraint_set_0}\eqref{sys_equations_0}, we can solve Problem \ref{equivalent_formulation_0} with the following two steps.
\begin{enumerate}
\item Solve $\emph{M}$ spline quantile regressions separately
\begin{align}\label{problem_statement_2}
\min_{ \bm{\lambda} } \; & \bm{\mathcal{E}}_{p_m} \left(\widehat{L}_m(\bm{a};\bm{Y}^\emph{N},\bm{x}^\emph{N}) \right), \quad m=1,\cdots,M \;.
\end{align}
\item Solve the systems of equations with constraints
\begin{equation}\label{second_step}
\begin{aligned}
\bm{a}& \in\mathcal{A} \\
\bm{Q} \bm{A} &= \bm{\Lambda} \;.
\end{aligned}
\end{equation}
\end{enumerate}
In both steps, the solution can be nonunique. The uniqueness of solution to quantile regression is discussed in \cite{koenker2005book}.
\cite{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$.
\subsection{Equivalence to CRPS Minimization}\label{sec_crps}
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{problem_1} is equivalent to CRPS minimization.
For a CDF $F$ and an observation $y_i$ of the response variable $y$, CRPS is defined by
\begin{align}
CRPS(F,y_i) = \int_\mathcal{R} \left( F(x) - 1_{\{y_i \leq x \} } \right)^2 \mathrm{d}x \;.
\end{align}
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
\begin{align}
CRPS(Q,y_i) = 2\int_0^1 \rho_p \left( y_i-Q(p) \right)\mathrm{d}p \;.
\end{align}
Consider calibrating the model $G(p,\bm{\theta},\bm{a})$ by minimizing the sum of CRPS of the response variables
\begin{align}\label{crps_min}
\min_{\bm{a}}\; \sum_{n=1}^\emph{N} \int_0^1 \rho_p \left( y_n-G(p,\bm{x}_n,\bm{a}) \right)\mathrm{d}p \;.
\end{align}
We see that \eqref{crps_min} is equal to the objective function in Problem \ref{problem_1} with uniform weight $w(p)$ by exchanging the integral and sum
\begin{align}
\int_0^1 \bm{\mathcal{E}}_p \left( y_n-G(p,\bm{x}_n,\bm{a}) \right)\mathrm{d}p
=
\int_0^1 \sum_{n=1}^\emph{N} \rho_p \left( y_n-G(p,\bm{x}_n,\bm{a}) \right) \mathrm{d}p \;.
\end{align}
When there is no basis quantile function, the CRPS minimization reduces the model to the least absolute deviation regression splines.
\section{Approximation Theorem}\label{sec_approximation}
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.
\begin{itemize}
\item $\mathcal{X}=$ hyperrectangle of domain of factors
\item $T= (t_{ij})_{i=1,\cdots,\emph{I+1}, j=1,\cdots,\emph{J} }=$ matrix representation of subdivision of $[0,1] \times \mathcal{X} $
\item $||T||=\max_{i,j}{|t_{i+1,j}-t_{i,j}|}$ mesh size of subdivision $Z$
\item $\mathrm{cone}(U)=$ cone of nonnegative linear combination of functions in a set of basis functions $U$
\item $\mathrm{int}=$ interior of a set
\item $\mathcal{P}(U,T)=$ set of piecewise functions where each piece (in
the scaled representation) defined on a subdivision defined by $Z$ is in $\mathrm{cone}(T)$
\item $G_{[0,1]\times \mathcal{X}}=$ cone of all continuous and nonnegative conditional quantile models defined on $[0,1]\times \mathcal{X}$ that is nondecreasing in the first variable
\end{itemize}
\begin{theorem}\label{theorem_approx}
Consider I-spline basis $\{Q_i(p)\}_{i=0}^\emph{I}$ defined on $[0,1]$ and B-spline basis $\{B_j(\bm{x})\}_{j=0}^\emph{J}$ defined on $\mathcal{X}$.
Furthermore, let $T_i$ be an asymptotically
nested sequence of subdivisions with mesh sizes approaching zero. Then the set $\cup_i \mathcal{P}(Q_i(p) \otimes B_j(\bm{x}), T_i)$ is a
dense subcone of $G_{[0,1]\times \mathcal{X}}$.
\end{theorem}
Theorem \ref{theorem_approx} is an application of \cite{Papp2011polynomial,
Papp2014spline}.
\section{Asymptotic Property}\label{sec_asym}
This section contains the asymptotic properties of the estimator.
\begin{itemize}
\item $\widehat{\mathcal{E}}_\emph{N}(\bm{a}) = \sum_{m=1}^\emph{M} w_m \bm{\mathcal{E}}_{p_m} \left( L(\bm{a};\bm{Y}^\emph{N},\bm{x}^\emph{N},p_m) \right)$
\item $\widehat{\bm{a}}_\emph{N} = \mathop{{\rm argmin}}_{\bm{a} \in \mathcal{A}} \widehat{\mathcal{E}}_\emph{N}(\bm{a})$
\item ${\bm{a}}_0=$ true parameter
\end{itemize}
Theorem \ref{theorem_consistency}, \ref{theorem_normalilty} in the following from \cite{Frumento2016WQR} are applications of \cite{Newey2}.
\begin{theorem}\label{theorem_consistency}
Assume that $\mathcal{A}$ is a compact set. If there is a function $\widehat{\mathcal{E}}_0(\bm{a})$ such that (i) $\widehat{\mathcal{E}}_N(\bm{a})$ converges uniformly in probability to $\widehat{\mathcal{E}}_0(\bm{a})$; (ii) $\widehat{\mathcal{E}}_0(\bm{a})$ is uniquely minimized by $\widehat{\bm{a}}$; (iii) $\widehat{\mathcal{E}}_0(\bm{a})$ is continuous at $\widehat{\bm{a}}$. Then $\widehat{\bm{a}}_\emph{N} \overset{d}{\rightarrow} {\bm{a}}_0$.
\end{theorem}
\begin{theorem}\label{theorem_normalilty}
Suppose that the conditions Theorem \ref{theorem_consistency} are satisfied, and (i) $\bm{a}_0$ is an interior point of $\mathcal{A}$; (ii) $\widehat{\mathcal{E}}_\emph{N}(\bm{a})$ is twice continuously differentiable in a neighborhood of $\bm{a}_0$; (iii) $\sqrt{n} \nabla_{\bm{a}} \widehat{\mathcal{E}}_\emph{N} (\bm{a}_0) \overset{d}{\rightarrow} \mathcal{N}(\bm{0},\Omega)$; (iv) there is $H(\bm{a})$ that is continuous at $\bm{a}_0$ and $\sup_{\bm{\theta}_0 \in \mathcal{A}} ||\nabla_{\bm{a}\bm{a}} \widehat{\mathcal{E}}_\emph{N} (\bm{a}_0) - H(\bm{a}) || \overset{p}{\rightarrow} \bm{0}$; (v) $H = H(\bm{a}_0)$ is nonsingular. Then
\begin{align}
\sqrt{N}(\widehat{\bm{a}}_\emph{N} -\bm{a}_0) \overset{d}{\rightarrow} \mathcal{N}(\bm{0},H^{-1} \Omega H^{-1} ) \;.
\end{align}
\end{theorem}
\section{Dimensionality Reduction by Reduced Rank Tensor}\label{sec_dim_reduction}
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.
\subsection{Tensor and CP Decomposition}
\begin{itemize}
\item $\bm{A}=$ parameter tensor
\item $K=$ number of scalar factors
\item $R=$ rank of parameter tensor
\item $\otimes =$ tensor product of two vectors
\item $\bigotimes_{k =1}^\emph{K}=$ tensor product of $\emph{K}$ vectors
\item $\cdot=$ sum of elements of the elementwise product of two tensors of the same size
\end{itemize}
\bigskip
A tensor is a multidimensional array. CP decomposition for a rank-$\emph{R}$ tensor $\bm{A}$ is defined as
\begin{align}\label{cp_decomp}
\bm{A} = \sum_{r=1}^\emph{R}
\bm{u}_r^{(0)} \otimes \cdots \otimes \bm{u}_r^{(\emph{K})} \;,
\end{align}
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
\begin{align}\label{tensor_dot_product}
\bigotimes_{k =0}^\emph{K} \bm{u}_r^{(k)}
\cdot
\bigotimes_{k =0}^\emph{K} \bm{v}_r^{(k)}
=
\prod_{k=0}^\emph{K}
\bm{u}_r^{(\emph{k})} \cdot \bm{v}_r^{(\emph{k})}\;,
\end{align}
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)
$.
\subsection{Reduced Rank Method for Parameter Tensor}
\begin{itemize}
\item $\emph{K}=$ number of scalar factors
\item $\emph{L} =$ number of basis spline functions for each scalar factor
\item $\emph{J} = \emph{K} \times \emph{L}$
\item $\bm{A}=$ parameter tensor with the same elements as parameter vector $\bm{a}$
\item $\bm{B}_k(x) = (B_{k1}(x),\cdots,B_{k\emph{L}}(x))=$ vector of univariate spline basis of scalar factor $x_k$
\item $\bigotimes_{k =1}^\emph{K} \bm{B}_k(x_k) =$ tensor product of $\emph{K}$ vectors of univariate spline basis of scalar factors
\item $\bm{Q}(p) = (Q_0(p),\cdots,Q_\emph{I}(p))=$ vector of basis quantile functions
\end{itemize}
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 \eqref{model_2} in tensor format
\begin{align}
G(p,\bm{x},\bm{a})
=
\bm{A} \cdot
\left(
\bm{Q}(p) \otimes
\bigotimes_{k =1}^\emph{K} \bm{B}_k(x_k)
\right)
\;.
\end{align}
Note that in our model formulation, $\bm{A}\geq 0$. Consider a nonnegative version of CP decomposition \eqref{cp_decomp} on the parameter tensor $\bm{A}$ where $\bm{u}_r^{(k)} \geq 0$. With \eqref{cp_decomp}\eqref{tensor_dot_product}, we have
\begin{align}\label{model_rank_decomp}
G(p,\bm{x},\bm{a})
=
\sum_{r=1}^\emph{R}
\left(
\bm{u}_r^{(0)} \cdot \bm{Q} \right)
\left(
\bm{u}_r^{(1)} \cdot \widehat{\bm{B}}_1
\right)
\cdots
\left(
\bm{u}_r^{(\emph{K})} \cdot \bm{B}_\emph{K}
\right)
\;.
\end{align}
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{first_method} is often sparse with many small nonzero values. Thus the low rank decomposition is expected to produce good approximation, although the decomposition \eqref{cp_decomp} 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.
\subsection{Alternating Algorithm}
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} \;.
$$
\begin{algorithm}[H]
\caption{Alternating Algorithm for Estimating Reduced Rank Parameter Tensor}
\begin{algorithmic}[1]
\REQUIRE Data $\{(y_i,\bm{x}_i)\}_{i=1}^\emph{N}$; confidence levels $\{p_m\}_{m=1}^\emph{M}$; model $G(p,\bm{\theta},{\bm{a}})$; error function $\bm{\mathcal{E}}_p$; threshold $\epsilon$
\STATE Initialize $\{\bm{U}^{(k)}\}_{k=0}^\emph{K}$
\REPEAT
\FOR{$k=0$ to $\emph{K}$}
\STATE Update $\bm{U}_k^{(s)}$ to $\bm{U}_k^{(s+1)}$ with fixed $\{\bm{U}_{k_t}^{(s)}\}_{k_t \neq k}$
\begin{equation}\label{problem_statement_iterative_2}
\begin{aligned}
\bm{U}_k^{(s+1)} = \mathop{{\rm argmin}}_{\bm{U}_k \geq 0} \; & \sum_{m=1}^\emph{M} w_m \bm{\mathcal{E}}_{p_m} \left( \widetilde{L}( \{\bm{U}_k^{(s)}\}_{k=0}^\emph{K} ; \bm{Y}^\emph{N},\bm{x}^\emph{N},p_m ) \right)
\end{aligned}
\end{equation}
\ENDFOR
\STATE $s\,+=1$
\UNTIL{
$
||\bm{A}^{(s+1)} - \bm{A}^{(s)}||_F \leq \epsilon
$
}
\ENSURE Parameter tensor $\bm{A}$
\end{algorithmic}
\end{algorithm}
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.
\section{Generalization Based on Risk Quadrangle Framework}\label{sec_quadrangle}
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.
\subsection{Risk Quadrangle}
\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$
\begin{align}
\mathrm{CVaR}_p(Z) = \Bar{q}_p(Z) = E[Z|Z>q_p(X)] \;.
\end{align}
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
\begin{align}
S(Z) = \arg \min_{C\in \mathcal{R} } \bm{\mathcal{E}}(Z-C) \;.
\end{align}
Two prominent examples are the quantile-based quadrangle and the superquantile-based Quadrangle.
Quantile-based quadrangle (at any confidence level $p\in (0,1)$):
\begin{align}\label{quantile_quad}
\mathcal{S}(Z)=\mathrm{VaR}_p(Z)=q_p(Z) = \text{quantile}, \quad
\bm{\mathcal{E}}_p(Z) = E[\frac{p}{1-p}Z_+ +Z_{-}]
=\text{normalized Koenker-Bassett error}.
\end{align}
Superquantile-based quadrangle (at any confidence level $p\in (0,1)$):
\begin{align}\label{cvar_quad}
\mathcal{S}(Z)=\mathrm{CVaR}_p(Z)=\Bar q_p(Z) = \text{superquantile}, \quad
\bm{\mathcal{E}}_p(Z) = \frac{1}{1-p}\int_0^1 \max \{0,\Bar q_\beta(Z)\}\mathrm{d}\beta - EZ.
\end{align}
\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$
\begin{align}\label{quadrangle_reg}
\min_{g\in\mathcal{C}} \bm{\mathcal{E}}(Y-g(\bm{x})) \;,
\end{align}
for given random variables $\bm{x}, Y$, and some given class $\mathcal{C}$ of functions $g$.
With $g$ obtained in \eqref{quadrangle_reg}, we can estimate the statistic of $Y$ conditioned on observing $\bm{x}$ by plugging in values of $\bm{x}$ to $g$.
\subsection{Generalization}
By replacing the error $\bm{\mathcal{E}}$ in Problem \ref{problem_1} 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.
\begin{example} \label{eg_factor}
Two-factor model of normal-logistic CVaR mixtures
\begin{equation}
\Bar{q}(p,\bm{x},\bm{a})
=
a_{00} + x_1 a_{01} + \Bar{Q}_\mathtt{N}(p) x_1 a_{11} + \frac{ H(p)}{1-p} x_2 a_{21} \;,
\end{equation}
where $\Bar{q}$ is a parameterized CVaR function,
$
H(p) = - p log(p) - (1-p) \log(1-p))
$, $ \frac{ H(p)}{1-p}$ is the CVaR function of standard logistic distribution, $\Bar{Q}_\mathtt{N}(p) = f_\mathtt{N}(\sqrt{2}{\mathop{{\rm erf}}}^{-1}(2p-1)) / (1-p) $ is the CVaR function of standard normal distribution, $f_\mathtt{N}$ is the density function of standard normal distribution
\citep{norton2021calculating}.
\end{example}
\section{Numerical Experiments}\label{sec_numerical}
This section presents two numerical experiments with real-world data.
\subsection{Factor Model of Mixture Quantiles}\label{sec_case_1}
The first experiment focuses on the Factor Model of Mixture Quantiles (Section \ref{sec_model}) and its calibration by convex optimization (Section \ref{sec_calibration}). 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 \cite{mccracken2016data} and downloaded from \cite{Fredcpi}. The data are standardized to have zero median and unit interquartile range. The data is randomly shuffled before experiments.
\paragraph{Model}
Let
\begin{itemize}
\item $Q_\mathtt{N}(p) = \sqrt{2}{\mathop{{\rm erf}}}^{-1}(2p-1) = $ quantile function of standard normal distribution
\item $Q_{\mathtt{E}_1}(p) = - \ln(4-4p)=$ quantile function of right-side transformed exponential distribution that is nonzero in $(0.75,1)$
\item $Q_{\mathtt{E}_1}(p) = \ln(4p)=$ quantile function of left-side transformed exponential distribution that is nonzero in $(0,0.25)$
\end{itemize}
Consider the Two-Factor Model of Mixture Quantiles
\begin{align}\label{model_experiment}
G(p,\bm{x},\bm{a})
=
f_0(x_1,x_2) +
Q_\mathtt{N}(p) f_1(x_1,x_2) +
Q_{\mathtt{E}_1}(p) f_2(x_1,x_2) +
Q_{\mathtt{E}_2}(p) f_3(x_1,x_2) \;,
\end{align}
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 \citep{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 \citep{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 \citep{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 \citep{Benaglia2009mixtools} and GAMLSS \citep{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{table_coverage} and Table \ref{table_crps}. 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.
\begin{table}[!htbp] \centering
\caption{}
\label{}
\begin{tabular}{@{\extracolsep{5pt}} cccccccc}
\\[-1.8ex]\hline
\hline \\[-1.8ex]
& 0.98 & 0.9 & 0.7 & 0.5 & 0.3 & 0.1 & Ave. Diff. \\
\hline \\[-1.8ex]
Factor model of mixture quantiles & $0.977$ & $0.895$ & $0.741$ & $0.508$ & $0.297$ & $0.089$ & $0.003$ \\
Gaussian mixture regression & $0.969$ & $0.907$ & $0.789$ & $0.603$ & $0.372$ & $0.113$ & $0.044$ \\
Generalized additive model & $0.975$ & $0.903$ & $0.728$ & $0.527$ & $0.315$ & $0.112$ & $0.011$ \\
\hline \\[-1.8ex]
\end{tabular}
\caption{This table presents the average coverage rates of the prediction intervals in 10-fold cross-validation for three models: Factor model of mixtures, Gaussian mixture regression, and generalized additive model. The average differences between the realized coverage rate and the target coverage rate are calculated. The target coverage rates are obtained by confidence intervals $(0.01,0.99)$, $(0.05,0.95)$, $(0.15,0.85)$,$(0.25,0.75)$,$(0.35,0.65)$,$(0.45,0.55)$. }
\label{table_coverage}
\end{table}
\begin{table}[!htbp]
\centering
\begin{minipage}{1\linewidth}
\centering
\begin{tabular}{ccccccc}
\\[-1.8ex]\hline
\hline \\[-1.8ex]
Model & 1 & 2 & 3 & 4 & 5 & 6 \\ \hline
Factor model of mixture quantiles & 0.199 & 0.241 & 0.236 & 0.220 & 0.204 & 0.219 \\
Gaussian mixture regression & 0.209 & 0.262 & 0.247 & 0.238 & 0.215 & 0.241 \\
Generalized additive model & 0.186 & 0.213 & 0.210 & 0.204 & 0.193* & 0.208* \\ \hline
\end{tabular}
\end{minipage}
\\
\begin{minipage}{1\linewidth}
\centering
\begin{tabular}{ccccccc}
\\[-1.8ex]\hline
\hline \\[-1.8ex]
Model & 7 & 8 & 9 & 10 & Mean & Std \\ \hline
Factor model of mixture quantiles & 0.197 & 0.255 & 0.193 & 0.186 & 0.215 & 0.023 \\
Gaussian mixture regression & 0.205 & 0.268 & 0.214 & 0.188 & 0.229 & 0.027 \\
Generalized additive model & 0.187* & 0.226 & 0.176 & 0.253 & 0.206 & 0.023 \\ \hline
\end{tabular}
\end{minipage}
\caption{This table presents the average out-of-sample Continuous Ranked Probability Score (CRPS) for three models: Factor model of mixtures, Gaussian mixture regression, and generalized additive model, obtained through a 10-fold cross-validation procedure. The mean and standard deviation of the CRPS values are reported for each model. The asterisk denotes instances in which a model produces predictions of nonexistent second moment and therefore nonexistent CRPS values for some data points. The CRPS is calculated only for predictions that have finite CRPS values. }\label{table_crps}
\end{table}
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{figure} shows three selected quantile surfaces for visualization, demonstrating the nonlinear relations incorporated in the model. The surfaces do not to cross each other.
\begin{figure}[!htbp]
\centering
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[scale=0.4]{ss1.png}
\caption{Quantile surface, confidence level $p=0.05$}
\label{fig_quantile_surfaces_1}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[scale=0.4]{ss2.png}
\caption{Quantile surface, confidence level $p=0.5$}
\label{fig_quantile_surfaces_2}
\end{subfigure} \\
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[scale=0.4]{ss3.png}
\caption{Quantile surface, confidence level $p=0.95$}
\label{fig_quantile_surfaces_3}
\end{subfigure}
\hfill
\caption{Quantile surfaces of confidence levels $0.05,0.5,0.95$ in the experiments in Section \ref{sec_case_1}. The black dots are the data points. }\label{figure}
\end{figure}
\subsection{Dimensionality Reduction by Reduced Rank Tensor for Factor Model of Mixture CVaRs}\label{sec_case_2}
The second numerical experiment examines the reduced rank method (Section \ref{sec_dim_reduction}) for the Factor Model of Mixture CVaRs (Section \ref{sec_quadrangle}). 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{sec_case_1} is used in the calibration.
\paragraph{Model}
Consider the Factor Model of Mixture CVaRs in tensor format
\begin{align}
G(p,\bm{x},\bm{a})
=
\bm{A} \cdot
\left(
\Bar{\bm{Q}}(p) \otimes
\bigotimes_{k =1}^\emph{K} \bm{B}_k(x_k)
\right)
\;,
\end{align}
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 \eqref{model_rank_decomp}, 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{problem_1} with the error defined in \eqref{cvar_quad}.
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{sec_case_1}, 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{fig_obj}. 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.
\begin{figure}[!htbp]
\centering
\includegraphics[scale=0.8]{Rplot_error_curve.pdf}
\caption{Stepwise objective value of alternating algorithm of rank $1$, $3$ and $10$ of experiments in Section \ref{sec_case_2}. }
\label{fig_obj}
\end{figure}
\clearpage
\bibliographystyle{unsrtnat}
\bibliography{FM_reference.bib}{}