Quantile Fourier Regression

• QFR is a modeling framework that combines quantile regression with fixed trigonometric bases to capture periodicity, oscillations, and nonlinear distributional structures.
• It is applied across various domains including functional regression, cyclostationary uncertainty modeling, extrapolative probabilistic forecasting via QFNN, and quantile spectral analysis.
• The methodology balances interpretability and flexibility, while addressing challenges like tail behavior, quantile crossing, and optimal basis selection.

Searching arXiv for recent and foundational uses of “Quantile Fourier Regression” and closely related quantile-Fourier formulations. Quantile Fourier Regression (QFR) denotes a class of quantile-based modeling strategies in which Fourier or trigonometric representations are used to describe distributional structure, periodic conditional quantiles, or frequency-domain dependence. On arXiv, the term appears in several closely related senses: as a Fourier-basis instantiation of quantile functional regression over probability levels; as a low-order Fourier model for cyclostationary uncertainty and time-inhomogeneous Markov decision processes; as the Quantile Fourier Neural Network for extrapolative probabilistic forecasting; and as trigonometric quantile regression underlying the quantile discrete Fourier transform, the quantile series, and estimators of quantile spectra (Yang et al., 2017, Khojaste et al., 2024, Hatalis et al., 2017, Li, 2022, Li, 2024). This suggests that QFR is best understood as a methodological family whose common feature is the use of quantile regression with Fourier structure, rather than as a single canonical estimator.

1. Scope and principal usages

The main uses of QFR differ by the object being modeled. In functional-distributional regression, the response is a subject-specific quantile function $Q_i(p)$ indexed by probability level $p$, and the Fourier basis is one option for expanding regression coefficient functions over $p$ (Yang et al., 2017). In cyclostationary modeling, the response is a time-varying conditional quantile $q_p(t)$ over a known period $T$, represented by low-order trigonometric polynomials and coupled to periodically varying Markov transitions (Khojaste et al., 2024). In probabilistic forecasting, the conditional quantile is treated as a deterministic function of time, learned by a single-hidden-layer Fourier neural network with time as the only input (Hatalis et al., 2017). In quantile spectral analysis, QFR means trigonometric quantile regression at a given frequency $\omega$ and quantile level $\tau$, producing the quantile discrete Fourier transform (QDFT) and the quantile series (QSER) (Li, 2022, Li, 2024).

Usage of QFR	Modeled object	Typical application
Fourier-basis quantile functional regression	$Q_i(p)$ over $p\in\mathcal P$	Distribution-on-covariates regression
Cyclostationary QFR	$q_p(t)$ over a known period $p$0	Periodic uncertainty and MDPs
QFNN	$p$1 with time-only input	Multi-step probabilistic forecasting
Spectral QFR	Trigonometric QR at $p$2	Quantile spectra and cross-spectra

These formulations share two structural commitments. First, they target quantiles or quantile-derived objects rather than conditional means. Second, they use Fourier bases to impose smooth periodic or oscillatory structure while avoiding a fully parametric specification of the outcome distribution.

2. Fourier bases within quantile functional regression

In the quantile functional regression framework, each subject $p$3 has a real-valued outcome distribution with CDF $p$4 and subject-specific quantile function

$p$5

The functional regression model treats the entire quantile function as the response: $p$6 with $p$7 and $p$8, where $p$9 is the sample size used to construct the empirical quantile for subject $p$0. The residual process $p$1 is assumed to be a mean-zero Gaussian process with covariance surface $p$2, independent across subjects (Yang et al., 2017).

Within this framework, QFR is the same quantile functional regression model but with a Fourier basis $p$3 replacing quantlets. The coefficient functions are expanded as

$p$4

On $p$5, a standard choice rescales $p$6 to $p$7 and uses

$p$8

so that $p$9 (Yang et al., 2017).

The Bayesian implementation proceeds by projecting each empirical quantile function onto the basis to obtain coefficient vectors $q_p(t)$0, collecting them into $q_p(t)$1, and fitting

$q_p(t)$2

with $q_p(t)$3 having i.i.d. rows $q_p(t)$4, often with $q_p(t)$5. Variance priors are $q_p(t)$6, and the basis coefficients $q_p(t)$7 receive spike-and-slab priors with cluster-specific hyperparameters. The MCMC algorithm updates $q_p(t)$8 and $q_p(t)$9 column-wise, and back-transformation yields $T$0 (Yang et al., 2017).

Inference is carried out at the level of whole functions. Simultaneous credible bands are constructed from posterior draws, and multiplicity-adjusted local inference is summarized through Simultaneous Band Scores (SimBaS), while the global Bayesian p-value is $T$1. The framework also supports posterior probability scores for differences in mean, variance, skewness, kurtosis-like fourth moments, and a “Gaussianness” diagnostic when the basis contains an explicit Gaussian subspace (Yang et al., 2017).

The original paper introduces quantlets rather than Fourier bases as the preferred basis system. Quantlets are described as sparse, regularized, near-lossless, and empirically defined, with the first two basis elements spanning the Gaussian quantile subspace and higher bases capturing departures such as kurtosis, skewness, and tails. The Fourier version preserves the same modeling core but changes the basis. Compared to quantlets, Fourier functions are fixed, global, and periodic; they may require larger $T$2 for tail behavior near $T$3, do not isolate a Gaussian subspace, and may produce endpoint artifacts unless trimming, cosine-only variants, or tapering are used (Yang et al., 2017).

3. Cyclostationary QFR for uncertainty modeling and Markov decision processes

A distinct formulation of QFR models periodic conditional quantiles of a stochastic process in time. For $T$4, the $T$5-quantile of $T$6 is represented as

$T$7

where the basis $T$8 spans periodic functions of period $T$9. The paper specializes to a Fourier basis with odd dimension $\omega$0,

$\omega$1

so that small $\omega$2 yields smooth, low-order trigonometric polynomials for $\omega$3 (Khojaste et al., 2024).

The setting is cyclostationary rather than stationary. A discrete-time process $\omega$4 is cyclostationary if the joint distribution of $\omega$5 is periodic in $\omega$6 with integer period $\omega$7, written in the paper as

$\omega$8

For a time-inhomogeneous Markov process, periodicity is encoded by a transition kernel satisfying

$\omega$9

This framework is motivated by seasonal or diurnal phenomena in energy systems, where stationary models or independent distributions fitted separately at each phase of the cycle can induce discontinuities such as month-end effects (Khojaste et al., 2024).

The time-varying quantiles are converted into a finite-state process by quantile partitioning. Given quantile levels $\tau$0, define

$\tau$1

By construction, the cyclostationary per-cycle probability of state $\tau$2 is $\tau$3 at every $\tau$4. Serial dependence is then modeled by periodically varying transition probabilities

$\tau$5

with row-stochasticity and optional column-sum constraints expressed directly on the Fourier coefficients. In the two-state example, validity for all $\tau$6 is enforced by

$\tau$7

together with linear constraints ensuring row sums equal one (Khojaste et al., 2024).

Estimation separates the quantile and transition components. For each quantile level $\tau$8, the coefficients $\tau$9 are obtained by minimizing Koenker’s check loss,

$Q_i(p)$0

which yields a linear programming problem. Transition coefficients $Q_i(p)$1 are estimated by maximizing the likelihood of observed state transitions, equivalently minimizing

$Q_i(p)$2

subject to probability and stationary-mass constraints for all $Q_i(p)$3. The resulting uncertainty model is embedded in a cyclostationary MDP with time-dependent costs $Q_i(p)$4, time-dependent transition probabilities $Q_i(p)$5, and an occupancy-measure LP for long-run average per-cycle cost (Khojaste et al., 2024).

Two numerical examples illustrate the framework. In the hydropower scheduling problem for the Waitaki River, weekly inflows from 1948–2021 are modeled with annual period $Q_i(p)$6, quantiles $Q_i(p)$7, QFR order $Q_i(p)$8, and a four-state Markov chain whose 16 transition probabilities are fitted as sinusoids with $Q_i(p)$9. The reservoir is discretized into energy blocks of $p\in\mathcal P$0 GWh, capacity is $p\in\mathcal P$1 GWh with 51 storage states, hydropower capacity is $p\in\mathcal P$2 MW, thermal capacity is $p\in\mathcal P$3 MW, costs are $p\in\mathcal P$4 for unserved demand, and the resulting MDP has $p\in\mathcal P$5 LP variables. The reported optimal policy exploits summer inflows to build storage and intentionally keeps storage high by early May to cover winter low inflows (Khojaste et al., 2024).

In the offshore wind integration example, hourly ISO-NE demand and buoy-based wind power are combined into a net-demand series with daily and annual commensurate periods. Using the multi-period QFR expansions, $p\in\mathcal P$6 clearly improves pseudo-$p\in\mathcal P$7, with $p\in\mathcal P$8 for $p\in\mathcal P$9 versus $q_p(t)$0–$q_p(t)$1 for $q_p(t)$2, and better captures diurnal winter double peaks. In 13 years of historical simulation, the time-dependent action plan yields total unmet demand $q_p(t)$3 and total cost $q_p(t)$4, while a fixed all-year plan yields total unmet demand $q_p(t)$5 and total cost $q_p(t)$6 (Khojaste et al., 2024).

4. QFR as extrapolative probabilistic forecasting: the Quantile Fourier Neural Network

In probabilistic forecasting of nonstationary univariate time series, QFR is operationalized by the Quantile Fourier Neural Network (QFNN), a single-hidden-layer Fourier neural network with sinusoidal activations and time as the only input. For quantile level $q_p(t)$7, the model is written as

$q_p(t)$8

where $q_p(t)$9 is the number of sinusoidal hidden units, $p$00 are learnable frequencies shared across quantiles, $p$01 are learnable phase shifts encoded as hidden-layer biases, $p$02 are quantile-specific amplitudes, $p$03 models nonperiodic components, and $p$04 is a quantile-specific bias (Hatalis et al., 2017).

The construction is called “composite quantile” because the network predicts multiple quantile levels simultaneously with shared hidden representation and linear output heads. It is explicitly extrapolative: because time is the only regressor, multi-step forecasts are obtained by direct evaluation at future time indices, without lagged values or covariates. The paper argues that this avoids error propagation across horizons and makes the same functional form naturally applicable out of sample (Hatalis et al., 2017).

Training uses the pinball loss

$p$05

and, because $p$06 is nondifferentiable at zero, a smooth approximation due to Zheng,

$p$07

For quantile set $p$08, the objective averaged over $p$09 samples and $p$10 quantiles is

$p$11

The implemented model does not impose explicit non-crossing constraints (Hatalis et al., 2017).

The reported training protocol is specific. Training times are normalized to $p$12; future times satisfy $p$13. If the maximum training value exceeds 10, values are scaled to $p$14. If multiplicative seasonality is present, a log transform is applied before training and inverted after prediction. Frequencies are initialized to multiples of $p$15, biases near 0, output weights near 1, and the augmentation unit is linear in all experiments. Optimization uses batch gradient descent with backpropagation in Keras/TensorFlow, smooth pinball loss with $p$16, learning rate $p$17 in all experiments except $p$18 for the Air Passengers median run, and a maximum of $p$19 training iterations. Dropout is applied to the hidden layer with rate selected by grid search in $p$20, and no $p$21 regularization is used on QFNN (Hatalis et al., 2017).

The paper evaluates QFNN on eight univariate datasets: Air Passengers, Sunspots, ISO New England real-time load, Internet traffic, Apple closing price, Solar power, Wind power, and a synthetic ocean wave elevation series. Nine benchmarks are used: Uniform method, Persistence method, ARIMA, SARIMA, ETS, linear quantile regression, polynomial quantile regression, composite SVQR, and composite QRNN. With 50% train / 50% test splits and long multi-step horizons, QFNN is reported to attain the lowest standardized quantile score on most datasets in the multiple-quantiles setting, specifically Air, Sunspots, Load, Internet, Stock, and Waves; Solar and Wind favor ETS and/or Persistence. QFNN generally produces the narrowest intervals except on Air Passengers, where Persistence is narrower, and the narrow intervals contribute to undercoverage, with average ACE around the low 20% range across datasets (Hatalis et al., 2017).

The paper presents QFNN as particularly effective when a series exhibits periodicity or quasi-periodicity, multiple seasonality, and smooth trend. Its limitations are equally explicit: because the input is time only, the model cannot predict unseen structural breaks or weather-driven intermittency from time alone; regime shifts and irregular bursts violate the extrapolation assumption; quantile crossing can occur because no monotonicity constraints are enforced; and the harmonic budget $p$22 must balance underfitting against overfitting (Hatalis et al., 2017).

5. QFR in quantile spectral analysis: QDFT, QSER, lag-window estimation, and SAR

In the quantile-spectrum literature, QFR means trigonometric quantile regression at a fixed frequency $p$23 and quantile level $p$24. For a series $p$25, the regression solves

$p$26

with

$p$27

This is the quantile analogue of least-squares Fourier regression and uses the standard quantile check loss $p$28 (Li, 2022, Li, 2024).

The fitted coefficients are then mapped to the quantile discrete Fourier transform. On the Fourier grid $p$29, the QDFT is

$p$30

The quantile periodogram is

$p$31

The inverse transform gives the quantile series,

$p$32

which is real-valued by conjugate symmetry and has sample mean equal to the sample $p$33-quantile (Li, 2022, Li, 2024).

The quantile series is important because it approximates the quantile-crossing process

$p$34

with $p$35. Writing $p$36, the quantile spectrum is defined as

$p$37

where $p$38 is the correlation of the level-crossing indicators at lag $p$39. Equivalently, $p$40 is the ordinary spectrum in $p$41 of the $p$42-indexed stationary process $p$43 (Li, 2024).

One route to estimation is nonparametric lag-window smoothing of the QSER autocovariance: $p$44 with Tukey–Hanning window

$p$45

The paper then smooths the resulting estimates across $p$46 using cubic smoothing splines or generalized additive mixed models. In a nonlinear mixture experiment with $p$47, the reported KLD is $p$48 at $p$49 without smoothing, $p$50 with smooth.spline and GCV, $p$51 with smooth.spline and $p$52, and $p$53 with gamm and correlated residuals. In a bivariate ARMA(2,1) experiment, the best mean KLD is $p$54 without smoothing at $p$55, $p$56 with GCV-based smooth.spline, $p$57 with smooth.spline and $p$58, and $p$59 with gamm (Li, 2022).

A more structured route is the spline autoregression (SAR) method, which fits a functional AR model to the centered QSER across a grid of quantile levels: $p$60 The coefficient matrices $p$61 are represented as smoothing splines of $p$62, estimated by penalized least squares,

$p$63

The resulting spectral estimator is

$p$64

The paper recommends selecting $p$65 by average AIC across quantiles and $p$66 by GCV (Li, 2024).

The SAR simulations are favorable relative to both unsmoothed AR and lag-window competitors. For the nonlinear mixture design with $p$67, SAR with GCV yields $p$68 versus $p$69 for AR without quantile smoothing and $p$70 for lag-window. For the ARMA design, the reported values are $p$71 for SAR, $p$72 for AR, and $p$73 for lag-window. The same paper also reports SAR-based Granger-causality on QSER, with bootstrap Wald tests detecting the designed causality at lag $p$74 in the nonlinear mixture and lags $p$75 and $p$76 in the ARMA example (Li, 2024).

6. Comparative interpretation, recurrent issues, and limitations

Across these literatures, the appeal of QFR is consistent. Fourier representations impose smooth periodic or oscillatory structure, quantile formulations avoid restricting attention to conditional means, and the resulting models can target distributional shape, tails, or quantile-dependent dependence structures directly. In the cyclostationary setting this yields smooth periodic quantile curves and periodically varying transition laws; in forecasting it yields direct multi-step quantile extrapolation; in spectral analysis it yields a quantile analogue of the ordinary DFT and periodogram; and in functional regression it yields a basis-expansion view of subject-specific distributions (Yang et al., 2017, Khojaste et al., 2024, Hatalis et al., 2017, Li, 2022, Li, 2024).

The recurring technical issues are also shared. Fourier bases are fixed, global, and periodic. In quantile functional regression they may require larger $p$77 to capture non-Gaussian features such as skewness, heavy tails, or asymmetric shoulders, and they do not isolate a Gaussian quantile subspace; by contrast, quantlets are constructed to be near-lossless and to contain an exact Gaussian subspace (Yang et al., 2017). In cyclostationary models, daily–annual interactions increase the number of terms quadratically in $p$78, which is why the reported practical guidance starts with $p$79 or $p$80 (Khojaste et al., 2024). In QFNN, the time-only input makes extrapolation principled when periodicity and trend dominate, but it cannot accommodate exogenous regime changes or intermittent weather-driven variability, and explicit non-crossing constraints are absent (Hatalis et al., 2017). In quantile-spectrum estimation, computation is heavier because trigonometric quantile regressions must be solved across many frequencies and quantiles; the AR order $p$81, the smoothing parameter $p$82, and the treatment of extreme quantiles materially affect performance (Li, 2022, Li, 2024).

Monotonicity and non-crossing occupy different roles across the variants. In the quantile functional regression formulation, monotonicity of the predicted quantile function is not explicitly imposed, although quantlet-based reconstructions and shrinkage are reported to produce predicted quantile functions that are typically monotone and should be checked for $p$83-monotonicity in practice (Yang et al., 2017). In the cyclostationary and QFNN formulations, the papers state that quantile curves are fitted separately and do not impose explicit non-crossing constraints, though both note that constrained multi-quantile fitting or derivative-based enforcement is possible (Khojaste et al., 2024, Hatalis et al., 2017).

Taken together, these uses establish QFR as a technical label for quantile methods that borrow the geometry of Fourier analysis without collapsing to mean-based harmonic modeling. The precise meaning depends on whether the index variable is probability level $p$84, physical time $p$85, forecast horizon encoded through time, or frequency $p$86. The unifying principle is the same: quantile structure is represented through trigonometric basis functions so that periodicity, oscillation, and distributional heterogeneity can be modeled simultaneously (Yang et al., 2017, Khojaste et al., 2024, Hatalis et al., 2017, Li, 2022, Li, 2024).