Online Distributional Regression

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.

134,527 characters · 15 sections · 72 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

\orcid{0009-0008-1409-9677} \email{[email removed]} \email{[email removed]}

\orcid{0000-0002-4944-9074} \email{[email removed]} \email{[email removed]} \orcid{0000-0002-2974-2660}

abstractLarge-scale streaming data are common in modern machine learning applications and have led to the development of online learning algorithms. Many fields, such as supply chain management, weather and meteorology, energy markets, and finance, have pivoted toward probabilistic forecasting. {\color{black} This results in} the need not only for accurate learning of the expected value but also for learning the conditional heteroskedasticity and conditional moments. Against this backdrop, we present a methodology for online estimation of regularized, linear distributional models. The proposed algorithm combines recent developments in online estimation of LASSO models with the well-known GAMLSS framework. We provide a case study on day-ahead electricity price forecasting, in which we show the competitive performance of the incremental estimation combined with strongly reduced computational effort. Our algorithms are implemented in a computationally efficient Python package {\color{black} ondil}.

Working Paper -- This Version: January 14, 2026

Introduction

Large-scale streaming data are common in modern applications of machine learning and have led to the development of online learning algorithms cesa2021online. For processes driven by a high-dimensional covariate space, regularized algorithms have been presented by, e.g., angelosante2010online, yang2010online, monti2018adaptive, yang2023online. In many settings, online statistical algorithms are used to issue forecasts. The advent of probabilistic forecasting in many fields, such as supply chain management, weather and meteorology, energy markets, and finance, yields the need not only for accurate learning of the conditional expected value but also for learning the conditional moments gneiting2014probabilistic, nowotarski2018recent, alvarez2021probabilistic, petropoulos2022forecasting, ziel2022m5. However, online learning approaches for distributional regression remain sparse in the literature and can be grouped as:

enumerate• Adaptive tracking of the {\color{black} variance-controlling} scale parameter alvarez2021probabilistic, vilmarest2024viking, • Adaptive estimation of the conditional heteroskedasticity priouret2005recursive, dahlhaus2007recursive, hendrych2018self, cipra2018robust,werge2022adavol, wintenberger2024stochastic, • Adaptive conformal prediction approaches zaffran2022adaptive, bhatnagar2023improved, gibbs2024conformal, dutot2024adaptive, brusaferri2024line,

and, to the best of our knowledge, there are no regularized online distributional regression approaches suitable for high-dimensional covariate processes available so far.

In this paper, we provide a regularized online learning algorithm for the conditional distribution parameters of the response variable $Y$ based on a combination of the online coordinate descent (OCD) algorithm introduced by angelosante2010online and messner2019online and the generalized additive models for location, scale and shape (GAMLSS) introduced by rigby2005generalized, stasinopoulos2008generalized, stasinopoulos2018gamlss. Formally, the GAMLSS framework assumes that {\color{black} the $N$-dimensional vector of response variables $\ensuremath{\boldsymbol{\mathbf{y}}} = \ensuremath{(y_1, \ldots, y_N)^\top}$} {\color{black} follows} the probability density function (PDF) $$f_{y} \left(y_n \mid \mu_n, \sigma_n, \nu_n, \tau_n \right) \quad {\color{black} \Leftrightarrow \quad f_Y\left(y_n \mid \theta_{n,1}, \theta_{n,2}, \theta_{n,3}, \theta_{n,4} \right)} $$ with (up to) four distribution parameters, each of which can be a (linear) function of explanatory variables. The first two parameters, $\mu_n$ and $\sigma_n$, commonly characterize the location and scale of the distribution, while $\nu_n$ and $\tau_n$ are commonly denoted as the shape parameters describing the skewness and kurtosis. {\color{black} Based on } rigby2005generalized we define the parametric GAMLSS as follows:

equation[equation omitted — 262 chars of source]

where $\mathcal{F}$ is {\color{black} a user-chosen} distribution with {\color{black} $K$ distribution parameters $\theta_{n,k}$ and $\ensuremath{\boldsymbol{\mathbf{\theta}}}_{n, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}} = (\theta_{n,1}, \ldots, \theta_{n,K})$.} {\color{black} Generally, we use subscript $\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}$ for indexing a column/row of a matrix. Depending on the domain of $Y$, the distribution's $\mathcal{F}$ support should match $\mathbb{R}$, the positive line $\mathbb{R}^+$, or an appropriate subset.} Let $g_k(\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}})$ be a known, monotonic link function relating a distribution parameter to a predictor $\ensuremath{\boldsymbol{\mathbf{\eta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}, k}$ by {\color{black}

equation[equation omitted — 509 chars of source]

where $\ensuremath{\boldsymbol{\mathbf{x}}}_j$ are the columns of the covariate matrix $\ensuremath{\boldsymbol{\mathbf{X}}}_k$, $\mathcal{J}_k$ is the index set for the linear effects, $\beta_{k,j}$ the distributional regression coefficients to be estimated. $b_{k,j}(\ensuremath{\boldsymbol{\mathbf{x}}}_{j})$ defines non-linear additive effects that can be represented by a linear combination of basis functions, $$b_{k,j}(\ensuremath{\boldsymbol{\mathbf{x}}}_{j}) = \sum^{D_{k,j}}_{l=1} \beta_{k,j,l} \phi_{k,j,l}(\ensuremath{\boldsymbol{\mathbf{x}}}_{j})$$ for some basis function $\phi_{k,j,l}(\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}})$ and the chosen set of additive effects $\mathcal{B}_k$; yielding the basis functions $\ensuremath{\boldsymbol{\mathbf{\phi}}}_{k,j}= \{\phi_{k,j,1}, \phi_{k,j,2},\ldots,\phi_{k,j,D_{k,j}}\}$ with $ |\ensuremath{\boldsymbol{\mathbf{\phi}}}_{k,j}| = D_{k,j}$. This class of additive effects includes e.g., B-Splines, or (smooth) ReLU-based splines muschinski2022cholesky, kneib2023rage, klein2024distributional and allows us to remain linear in the coefficients. Therefore, without loss of generality, we can rewrite the model in Equation (ref) as

equation[equation omitted — 309 chars of source]

where the model matrix $\ensuremath{\boldsymbol{\mathbf{X}}}_k = (\ensuremath{\boldsymbol{\mathbf{1}}}, (\ensuremath{\boldsymbol{\mathbf{x}}}_{j})_{j\in \mathcal{J}_k}, (\ensuremath{\boldsymbol{\mathbf{\phi}}}_{kjl}(\ensuremath{\boldsymbol{\mathbf{x}}}_j) )_{j\in \mathcal{B}_k, l=1,\ldots, D_{k,j} } )$ is of size $N \times J_k$ where $J_k = 1+ |\mathcal{J}_k| + \sum_{j\in \mathcal{B}_k} D_{k,j}$ and $\ensuremath{\boldsymbol{\mathbf{\beta}}}_k$ is the according linear coefficient vector and $|\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}|$ is the cardinality.} The {\color{black} parametric GAMLSS model} therefore allows the modeling of all conditional distribution parameters as {\color{black} non-}linear, additive functions of the explanatory variables in $\ensuremath{\boldsymbol{\mathbf{X}}}_k$. {\color{black} rigby2005generalized introduced estimation via iteratively reweighted least squares (IRLS), which has been extended to LASSO-type penalties by groll2019lasso. Further regularized estimation approaches have been proposed by ziel2021gamlss.} {\color{black} Alternatively, estimation via Bayesian approaches (direct and via Bayesian-IRLS), direct likelihood minimization using stochastic gradient descent thielmann2024neural,umlauf2024scalable and gradient-boosted decision trees maerz2019xgboostlss, cevid2022distributional are popular. }

figure[figure omitted — 745 chars of source]

{\color{black} The aforementioned methods concern the batch estimation of distributional regression models. However, in many applications, we regularly receive new data. Formally,} after having received the first $N$ observations, we receive a new pair set of data $y_{N+1}$ and $\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1, k}$. In the online setting, we are interested in updating the {\color{black} coefficients} of the model $\ensuremath{\boldsymbol{\mathbf{\beta}}}_k$, without recalling all previous observations. We {\color{black} hence} consider an algorithm for a strict online setting in which we can discard all observations after updating the model coefficients. Our proposed algorithm {\color{black} builds on} the IRLS algorithm proposed by rigby2005generalized and exploits its agnosticism towards the estimation method {\color{black} within each iteration. This allows us to exchange the used batch estimation methods with online estimation methods --- including sparse and regularized estimation approaches --- and thereby propose an algorithm suitable for streaming data.}

Figure (ref) gives an idea of the building blocks of our algorithm. In summary, we make the following contributions to the current literature:

enumerate• We propose {\color{black} an online and regularized algorithm enabling} scalable distributional regression for {\color{black} streaming data. We implement LASSO regularization ($L_1$-regularization), ridge-based ($L_2$-regularization) and the elastic-net for combined $L_1$ and $L_2$ penalties.} We discuss the implementation of the incremental update step, including regularization and online model selection in-depth (Section (ref)). • We formally analyze the relationship between the online estimation and repeated batch estimation and discuss the influence of the exponential discounting, the running time, and the initial sample size. Due to the generic setting of our algorithm, exact results depend on the chosen distributional assumption and link functions. Nevertheless, we aim to provide some intuitive guidance and validate our results in an extensive simulation study (Section (ref)). • We validate the proposed method in a forecasting study for electricity prices and demonstrate the competitive performance of our model, combined with a strongly reduced estimation time compared to batch estimation. {\color{black} Furthermore, we conduct ablation studies to analyze the impact of several algorithmic and hyperparameters in a real-world example} (Section (ref)). • We provide an open-source, ready-to-use Python implementation of the online distributional regression model {\color{black} in the package ondil}. Our implementation is based on numpy and scipy and employs \texttt{numba} just-in-time compilation for high computational efficiency. The code can be accessed on \texttt{GitHub},\footnote{See: \url{https://github.com/simon-hirsch/ondil}.} and the package is available on PyPI\footnote{See: \url{https://pypi.org/project/ondil}.} {\color{black} (Section (ref))}.

Our work opens multiple avenues for future research. First, our implementation is constrained to {\color{black} parametric GAMLSS models, and hence only allows for non-linear additive effects that can be expressed as a sum of basis functions and linear coefficients. While this class includes a wide array of non-linear effects and can be estimated in a regularized way by the proposed sparse estimation methods,} the inclusion of smooth terms, random effects, and regularized P-splines seems like a worthwhile extension rigby2005generalized, eilers1996flexible. Furthermore, the inclusion of autoregressive and cross-moments effects can provide valuable tools for modeling highly complex {\color{black} time series} processes such as weather, electricity markets, or supply chain applications. Lastly, a thorough theoretical analysis of the error bounds of the proposed online approximation compared to the batch setting should further increase the trust in the presented methods.

{\color{black} The remainder of this paper is structured by the {\color{black} flow given in Figure (ref).} The following two sections give the building blocks of our algorithm: First, Section (ref) briefly reviews the algorithm for the parametric GAMLSS. Afterwards, Section (ref) reviews online estimation of regularized linear models. Section (ref) presents our main contribution, the online distributional regression. The following Sections (ref) and (ref) present a simulation study on the properties of the online algorithm and our real-world application for energy markets. Section (ref) presents our open-source implementation. Finally, Section (ref) discusses our results and concludes the paper.}paper.}

{\color{black}

Batch Distributional Regression

} {\color{black} We denote scalar values as lowercase letters $a$, vectors as bold lowercase letters $\ensuremath{\boldsymbol{\mathbf{a}}} = (a_1, \ldots, a_N)$ and matrices as bold uppercase letters $\ensuremath{\boldsymbol{\mathbf{A}}} = (\ensuremath{\boldsymbol{\mathbf{a}}}_1, \ldots, \ensuremath{\boldsymbol{\mathbf{a}}}_M)$. For an arbitrary matrix $\ensuremath{\boldsymbol{\mathbf{A}}}$, $\ensuremath{\boldsymbol{\mathbf{A}}}_{n, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}}$ denotes accessing the $n$-th row, while $\ensuremath{\boldsymbol{\mathbf{A}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}, i}$ denotes accessing the $i$-th column. We denote element-wise or Hadamard multiplication using $\odot$ and element-wise division using $\oslash$, that is $\ensuremath{\boldsymbol{\mathbf{a}}} \oslash \ensuremath{\boldsymbol{\mathbf{b}}} = \ensuremath{\boldsymbol{\mathbf{a}}} \odot \ensuremath{\boldsymbol{\mathbf{b}}}^{-1}$ for vectors.} {\color{black} The following exposition of the IRLS algorithm for fitting distributional regression models largely follows rigby2005generalized and stasinopoulos2008generalized. The full algorithm} is outlined in Algorithm (ref). {\color{black} Our algorithm follows the RS-Algorithm (Rigby & Stasinopoulos Algorithm), although an implementation of the GC-Algorithm (Cole & Green Algorithm) is possible as well.} {\color{black} {\color{black} We start with the distributional regression model from the introduction given in Equations (ref) and (ref):

equation*[equation* omitted — 466 chars of source]

} and, without loss of generality, we assume that the same set of covariates is used for all distribution parameters $k \in 1, \ldots, K$ and therefore $\ensuremath{\boldsymbol{\mathbf{X}}}_k = (\ensuremath{\boldsymbol{\mathbf{1}}}, (\ensuremath{\boldsymbol{\mathbf{x}}}_{j})_{j\in \mathcal{J}_k}, (\ensuremath{\boldsymbol{\mathbf{\phi}}}_{k,j,l}(\ensuremath{\boldsymbol{\mathbf{x}}}_j) )_{j\in \mathcal{B}_k, l=1,\ldots, D_{kj} } )$. Let $\ensuremath{\boldsymbol{\mathbf{\theta}}}_k = \ensuremath{\boldsymbol{\mathbf{\theta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}},k}$, $\ensuremath{\boldsymbol{\mathbf{\eta_k}}} = \ensuremath{\boldsymbol{\mathbf{\eta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}},k} = \ensuremath{\boldsymbol{\mathbf{\eta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}},k}(\ensuremath{\boldsymbol{\mathbf{\theta}}}_k)$ and let $\ell(y_n \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}_{n, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}}) = \log(\mathcal{L}(y_n \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}_{n, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}}))$ denote the log-likelihood function of the assumed distribution $\mathcal{F}$. We use the shorthand $\ensuremath{\boldsymbol{\mathbf{\ell}}} = \ensuremath{\boldsymbol{\mathbf{\ell}}}(\ensuremath{\boldsymbol{\mathbf{\theta}}}) =\ell(\ensuremath{\boldsymbol{\mathbf{y}}} \mid \ensuremath{\boldsymbol{\mathbf{\theta}}})$. Then

equation[equation omitted — 471 chars of source]

}tion} } denotes the score vectors. Generally, partial derivatives are evaluated element-wise, hence $\partial \ensuremath{\boldsymbol{\mathbf{\ell}}} / \partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k$ is a vector of length $N$ with the elements $\partial \ell(y_n \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}_{n, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}}) / \partial \eta_{nk}$. Let

equation[equation omitted — 211 chars of source]

be the adjusted dependent variables, also called the working vector. $\ensuremath{\boldsymbol{\mathbf{W}}}_{k,s}$ is a diagonal, iterative weight matrix, which can have one of the following forms: {\color{black}

align[align omitted — 1,087 chars of source]

giving the observed (Newton-Raphson), expected (Fisher-Scoring) or product (quasi Newton-Raphson) score function rigby2005generalized. Note that differentiation is again evaluated element-wise and the full Hessian tensor $\ensuremath{\boldsymbol{\mathbf{H}}}$ is of size $N \times K \times K$. Combining Equations (ref) and the RS and CG algorithms gives six possible combinations of the RS and CG algorithms and the three score types. In the following, we focus on the RS algorithm with Fisher scoring, but the other combinations can be implemented similarly. Using the chain rule and Fisher's identity, we have that:

equation*[equation* omitted — 866 chars of source]

The Fisher information matrix can be recovered as the variance of the score function or as the second derivative of the log-likelihood with respect to the predictor. Hence, the adjusted observation vector or working vector $\ensuremath{\boldsymbol{\mathbf{z}}}_k$ reads:

align[align omitted — 1,284 chars of source]

where we note that the weight matrix $\ensuremath{\boldsymbol{\mathbf{W}}}_{k,k}$ evaluated at the diagonal only. Therefore its inverse corresponds to $\operatorname{diag}(\ensuremath{\boldsymbol{\mathbf{W}}}_{k,k})^{-1}$ and the multiplication with the score vector is element-wise.} In the implementation, we regularly employ that:

equation[equation omitted — 273 chars of source]

where $g_k^{-1}(\ensuremath{\boldsymbol{\mathbf{\eta}}}_k) = \ensuremath{\boldsymbol{\mathbf{\theta}}}_k$ is the inverse of the link function. For pseudo-Newton methods, similar results hold. This allows for a mix-and-match implementation of link functions and distributions, since we only need the Jacobian, potentially the Fisher Information and the first derivative of $g_k^{-1}(\cdot)$ rigby2005generalized, stasinopoulos2008generalized, stasinopoulos2024generalized, however, hirsch2025online show that true Newton-Raphson methods additionally need the second derivative of $g_k^{-1}(\cdot)$ to allow for a similar mix-and-match implementation of the link functions and distributions.

algorithm[algorithm omitted — 2,922 chars of source]

{\color{black} The iteratively reweighted least squares algorithm rigby2005generalized for the estimation of distributional regression models is given in Algorithm (ref) and consists of two nested cycles. Let $r$ and $i$ be the indices for the outer and inner cycle. The outer cycle maximizes the penalized likelihood with respect to $\ensuremath{\boldsymbol{\mathbf{\beta}}}_k$. Within each outer cycle, the algorithm iterates through the distribution parameters $k = 1, \ldots, K$. For each distribution parameter $\ensuremath{\boldsymbol{\mathbf{\theta}}}_k$ the inner cycle consists of repeated, weighted regression of the score {vector $\ensuremath{\boldsymbol{\mathbf{z}}}_k^{[r,i]}$} on the design matrix $\ensuremath{\boldsymbol{\mathbf{X}}}_k$ using the iterative weights $\ensuremath{\boldsymbol{\mathbf{W}}}_{k,k}^{[r,i]}$. We denote iterations as superscript, i.e. $\ensuremath{\boldsymbol{\mathbf{\theta}}}^{[r, i]}_k$ is the value of $\ensuremath{\boldsymbol{\mathbf{\theta}}}_k$ at the outer iteration $r$ and the inner iteration $i$. The weighted fit in Line 8 is at the core of the Newton step for the estimation of the conditional distribution parameters. It is important to see that the general algorithm is agnostic to which (weighted) statistical learning method is employed to regress $\ensuremath{\boldsymbol{\mathbf{z}}}_k^{[r,i]}$ on the design matrix $\ensuremath{\boldsymbol{\mathbf{X}}}_k$ using the iterative weights $\ensuremath{\boldsymbol{\mathbf{W}}}_{k,k}^{[r,i]}$ stasinopoulos2024generalized. This has facilitated the development of multiple extensions to the batch GAMLSS, such as gamlss.boost mayr2012generalized, hofner2014gamboostlss, mixtures with Neural Networks rugamer2024semi, and regularized approaches as in groll2019lasso and gamlss.lasso ziel2021gamlss. We do not include the backfitting iteration here as we constrain ourselves to non-linear additive effects that can be represented as a combination of basis functions of $\ensuremath{\boldsymbol{\mathbf{X}}}_k$.}

{\color{black} Lastly, let us discuss two important aspects of using iteratively reweighted least squares for the estimation of distributional regression models using a concrete example. The calculation of the score vector, working vector, and weights depends on the distributional assumption and the chosen link function. Remark (ref) introduces a simple Gaussian example and gives the concrete working vectors and weights.

remarkWhen considering a distributional regression model for the normal distribution with density $$ f(x| \mu, \sigma) = 1 / \sqrt{2 \pi \sigma^2} \exp\left(-1/2 (x-\mu)^2/ \sigma^2 \right), $$ and using a parameterization of $\theta_1 = \mu$, $\theta_2= \sigma^2$ and identity link functions $g_1(\theta)=\theta$ and $g_2(\theta)=\theta$, the iteratively reweighted least squares algorithm uses the following working vectors and weights \begin{align} \ensuremath{\boldsymbol{\mathbf{z}}}_1 &= \ensuremath{\boldsymbol{\mathbf{z}}}_{\mu} = \ensuremath{\boldsymbol{\mathbf{y}}} \quad and & \ensuremath{\boldsymbol{\mathbf{W}}}_{1,1} &= \operatorname{diag} \left((\ensuremath{\boldsymbol{\mathbf{\sigma}}}^2)^{-1}\right) \\ \ensuremath{\boldsymbol{\mathbf{z}}}_2 &= \ensuremath{\boldsymbol{\mathbf{z}}}_{\sigma^2} = (\ensuremath{\boldsymbol{\mathbf{y}}} - \ensuremath{\boldsymbol{\mathbf{\mu}}})^2 \quad and &\ensuremath{\boldsymbol{\mathbf{W}}}_{2,2} &= \operatorname{diag}\left((2\ensuremath{\boldsymbol{\mathbf{\sigma}}}^4)^{-1}\right). \end{align}

The working vectors $\ensuremath{\boldsymbol{\mathbf{z}}}_1, \ensuremath{\boldsymbol{\mathbf{z}}}_2$ and weights $\ensuremath{\boldsymbol{\mathbf{W}}}_{1,1}$ equal the working vectors used to estimate mean-variance models for conditional heteroskedasticity using least-squares estimation, where $\ensuremath{\boldsymbol{\mathbf{W}}}_{2,2}$ is often omitted and set as constant $\ensuremath{\boldsymbol{\mathbf{W}}}_{2,2} = 1$ francq2019garch. For this approach, consistent model selection using the adaptive LASSO and asymptotic normality of the estimators have been shown by dette2013least for independent observations and ziel2016iteratively for time series settings. The weights $\ensuremath{\boldsymbol{\mathbf{W}}}_{2,2}$ in Remark (ref) additionally correct for heteroskedasticity in the estimation of the conditional scale parameter. groll2019lasso discuss LASSO estimation in the general setting of distributional regression in an application and simulation study; however, general results on consistent model selection are, to the best knowledge of the authors, not available. }

The working vectors (Equation (ref) in Remark (ref)) also underscore an important fact in distributional modeling: {\color{black} The models for the conditional distribution parameters depend on each other. Therefore, errors due to low-quality models for one distribution parameter can propagate through the distributional model and adversely affect the estimation of other distribution parameters. This can lead to hard-to-predict effects in the estimation, if the distribution parameters are not information-orthogonal cox1987parameter. There is a directional component to this issue, as there is no hope in estimating the conditional scale correctly if the model for the location is misspecified, while the impact of a misspecified scale model on the location is not as severe. For models with conditional skewness and tail behavior, the issue can be even more severe, as the relationships might be highly non-linear. This issue is discussed in detail in the context of the M5 competition in ziel2022m5.}

Online Estimation of Regularized Linear Models

Let us now move to the online or streaming setting. For the sake of simplicity, we assume we are in a regular regression {\color{black} to the mean} setting and omit the subscript $k$ for the distribution parameter used in the previous subsection. The issue at hand can be summarized as follows: given some data $\ensuremath{\boldsymbol{\mathbf{y}}}_N=\ensuremath{(y_1,\ldots, y_N)^\top}$, $\ensuremath{\boldsymbol{\mathbf{X}}}_N = (\ensuremath{\boldsymbol{\mathbf{x}}}_1, \ldots, \ensuremath{\boldsymbol{\mathbf{x}}}_N)$ and weights $\ensuremath{\boldsymbol{\mathbf{W}}}_N = \operatorname{Diag}\left(w_1, \ldots, w_N \right)$ we have estimated a set of coefficients $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_N$. We are now interested in recovering the coefficients $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1}$, given a new observation for $y_{N+1}$, a new row of $\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}$ and a new weight $w_{N+1}$. {\color{black} In slight abuse of notation, we employ the subscripts $N$, and $N+1$ to indicate which data is available to the forecaster by now, thereby $\ensuremath{\boldsymbol{\mathbf{y}}}_{N+1} = (y_1, \ldots, y_N, y_{N+1})$ and $\ensuremath{\boldsymbol{\mathbf{X}}}_{N+1} = (\ensuremath{\boldsymbol{\mathbf{x}}}_1, \ldots, \ensuremath{\boldsymbol{\mathbf{x}}}_N, \ensuremath{\boldsymbol{\mathbf{x}}}_{N+1})$. We employ this notation exclusively around the updating step $N \rightarrow N+1$.} The following two sections review the estimation of weighted linear models using recursive least squares (Section (ref)) and the incremental estimation of $L_1$ and $L_2$ penalized linear models using online coordinate descent (OCD, see Section (ref)).

Weighted Recursive Least Squares

The following is a condensed introduction haykin2014adaptive {\color{black} to recursive estimation for weighted least squares estimation}. {\color{black} Nevertheless, it is useful to emphasize a few points that will be of interest for the online distributional learning algorithm introduced in Section (ref).} The weighted least squares regression problem estimates the coefficients $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_N^\text{WLS}$ that minimize the loss:

equation[equation omitted — 644 chars of source]

by noting that $\ensuremath{\boldsymbol{\mathbf{W}}}_N$ is invertible as long as all weights are positive and non-zero. {\color{black} Define the weighted Gramian $\ensuremath{\boldsymbol{\mathbf{G}}}_N = \ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}_N^\top}\ensuremath{\boldsymbol{\mathbf{W}}}_N\ensuremath{\boldsymbol{\mathbf{X}}}_N$.} Under the assumptions that (1) the weights are known and (2) do not change for past observations, the weighted least-squares problem can be re-formulated as a recursive algorithm for the next observations $y_{N+1}$ and covariates. {\color{black} These assumptions are in general not an issue if one is using sample weights. We will return to these assumptions in Section (ref) for the case of estimation weights. In either case, the weights should be living on a consistent scale. By employing the Sherman-Morrison formula (see Appendix (ref)), we can efficiently update the Gramian matrix by taking

align[align omitted — 549 chars of source]

and we receive

align[align omitted — 793 chars of source]

where the latter Equation (ref) gives a convenient update equation for streaming settings.} An important special case is the recursive least squares with exponential forget, which minimizes the loss {\color{black}

equation[equation omitted — 383 chars of source]

where $$\ensuremath{\boldsymbol{\mathbf{\Gamma}}}_N = \operatorname{Diag}\left((1-\gamma)^{N-1}, (1-\gamma)^{N-2}, \ldots,(1-\gamma)^{1}, (1-\gamma)^{0} \right)$$ and $\gamma$ denotes the forget factor. We} note that the exponential discounting of older observations leads to

equation[equation omitted — 429 chars of source]

and the {\color{black} corresponding} update equation for the exponentially discounted inverted Gram matrix

equation[equation omitted — 1,178 chars of source]

which allows for an efficient online update of the OLS coefficients under exponential forgetting by using Equation (ref) and setting $w_{N+1} = 1$. {\color{black} Note the subtle difference between Equation (ref) for updating the inverse Gramian matrix with sample weights $\ensuremath{\boldsymbol{\mathbf{W}}}_N$ and Equation (ref) using discounting weights $\ensuremath{\boldsymbol{\mathbf{\Gamma}}}_N$, as the sample weights are constant, but the unique structure of exponential discounting allows the weights $(1-\gamma)^n$ to depend on the distance to the current newest observation.}

Lastly, we combine sample weights and exponential discounting: For the exponentially discounted, weighted Gram matrix $\ensuremath{\boldsymbol{\mathbf{G}}}_N = (\ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}_N^\top} \ensuremath{\boldsymbol{\mathbf{\Gamma}}}_N \ensuremath{\boldsymbol{\mathbf{W}}}_N \ensuremath{\boldsymbol{\mathbf{X}}}_N)$, we have

equation[equation omitted — 259 chars of source]

and accordingly the update of the inverted Gram matrix $\ensuremath{\boldsymbol{\mathbf{G}}}_{N+1}^{-1}$ can be written as:

equation[equation omitted — 613 chars of source]

which can be plugged into Equation (ref) to update the coefficients $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}$, as in Algorithm (ref).

algorithm[algorithm omitted — 853 chars of source]

Online Coordinate Descent

{\color{black} The recursive least squares estimator can be numerically unstable due to its reliance on the inverted Gramian matrix $\ensuremath{\boldsymbol{\mathbf{G}}}_N^{-1}$. Furthermore, for regression problems with large sets of potential covariates, some regularization is commonly employed to keep a parsimonious model. The elastic net friedman2010regularization, friedman2007pathwise minimizes a combination of $L_1$ (or LASSO) and $L_2$ (or Ridge) penalties on the coefficient vector and is commonly written as:

equation[equation omitted — 548 chars of source]

where the ratio between $L_1$ and $L_2$ penalties is governed by $\alpha$. From the elastic net, we can recover ordinary least squares ($\lambda = 0$), the LASSO ($\alpha=1$), and the ridge regression ($\alpha=0$). In the batch setting, } friedman2007pathwise, friedman2010regularization introduced path-wise cyclic coordinate descent (CCD) for elastic net regression problems, {\color{black} and due to its computational efficiency, it remains the state-of-the-art method for coefficient estimation}. For given regularization parameters {\color{black} $\lambda_1$ and $\lambda_2$}, we repeatedly update the coefficient vector $\ensuremath{\boldsymbol{\mathbf{\beta}}}_N = (\beta_{N,1}, \ldots, \beta_{N,J})$ of length $J$ and update {\color{black}

equation[equation omitted — 232 chars of source]

} where $S(\beta, \lambda) = \text{sign}(\beta) \max(\lvert \beta \rvert - \lambda, 0)$ is the soft-thresholding function {\color{black} and $y_n - \tilde{y}_n^{(j)} = y_n - \widehat{y}_n + x_{nj}\widehat{\beta}_j$ is the partial residual excluding the contribution of $x_{nj}$}. Note that, implicitly, friedman2010regularization assume that $\sum_{n=1}^N w_n = 1$. We use a numerical convergence criterion to {\color{black} terminate} the algorithm after convergence.

As proposed in angelosante2010online and messner2019online for LASSO estimation, we rewrite the update Equation (ref) to calculate the updated parameter from the weighted Gramian matrix $\ensuremath{\boldsymbol{\mathbf{G}}}_N = \ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}^\top}_N\ensuremath{\boldsymbol{\mathbf{W}}}_N\ensuremath{\boldsymbol{\mathbf{X}}}_N$ and the {\color{black} vector $\ensuremath{\boldsymbol{\mathbf{h}}}_N = \ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}^\top}_N\ensuremath{\boldsymbol{\mathbf{W}}}_N\ensuremath{\boldsymbol{\mathbf{y}}}_N$} and their potentially exponentially discounted counterparts (see Section (ref) and Equation (ref)). {\color{black} Since $\ensuremath{\boldsymbol{\mathbf{G}}}_N$ and $\ensuremath{\boldsymbol{\mathbf{h}}}_N $ can be updated online, this allows an efficient incremental estimation on a data stream.} {\color{black} We observe for the elements $G_{N,j,j}$ and $h_{N,j}$ of $\ensuremath{\boldsymbol{\mathbf{G}}}_{N}$ and $\ensuremath{\boldsymbol{\mathbf{h}}}_{N}$ that

align*[align* omitted — 117 chars of source]

Thus, writing the nominator inside the soft-thresholding function as

align*[align* omitted — 670 chars of source]

and} we receive the {\color{black} update equation for the elastic net as} {\color{black}

equation[equation omitted — 383 chars of source]

} Therefore, the algorithm only needs to store the $J \times J$ {{matrix $\ensuremath{\boldsymbol{\mathbf{G}}}_N$}} and the {{$J$-dimensional vector\;$\ensuremath{\boldsymbol{\mathbf{h}}}_N$}}. {\color{black} Once a new pair of observations arrives, we can update $\ensuremath{\boldsymbol{\mathbf{G}}}_{N+1}$ and $\ensuremath{\boldsymbol{\mathbf{h}}}_{N+1}$ and re-run the cyclic coordinate descent based on the previous estimated coefficients.} {\color{black} As in batch} coordinate descent, we can run the algorithm for a decreasing sequence of regularization strengths $\lambda$, starting with $\lambda_\text{max} = \max{}\lvert \ensuremath{\boldsymbol{\mathbf{h}}}_N \rvert$ as the element-wise maximum in $\ensuremath{\boldsymbol{\mathbf{h}}}_N$ and using an exponential grid towards $\lambda_\text{min} = \epsilon \lambda_\text{max}$ with $\epsilon_\lambda = 0.001$ as typical values friedman2010regularization, friedman2007pathwise, messner2019online. {\color{black} Note that, as $\ensuremath{\boldsymbol{\mathbf{h}}}_N$ changes in the online setting with each new pair of observations, the regularization grid varies slowly.} Algorithm (ref) presents the algorithm schematically.

algorithm[algorithm omitted — 1,594 chars of source]

{\color{black} We extend the formulation to allow for the possibility to box-constrain $L_j \leq \beta_j \leq U_j$ with user-chosen bounds $L_j$ and $U_j$. As in the batch implementation in the glmnet R package friedman2010regularization, tay2023elastic, Equation (ref) is then formulated as:

equation*[equation* omitted — 408 chars of source]

}

Online Distributional Regression

Having established the batch IRLS algorithm for the estimation of distributional regression models and the online coordinate descent for regularized linear models, this section proceeds by putting the pieces together and presents the online linear {\color{black} distributional regression} model. {\color{black} Algorithm (ref) outlines the main procedure and we will subsequently describe the details with respect to model selection, mini-batch updates, stopping criteria, and the possible approximation error towards the batch case. We present the algorithm using a path-based, regularized estimation method in each inner iteration, and therefore include model selection; however, using RLS or OCD for unregularized estimation is possible as well.}

Update Algorithm

{\color{black} We start again with the distributional regression model from the introduction given in Equation (ref) and (ref):

equation*[equation* omitted — 458 chars of source]

} After having seen $N$ observations, we have estimated coefficients $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N, k}$ for $k = 1, .., K$ distribution parameters. Given some new data $(y_{N+1}, \ensuremath{\boldsymbol{\mathbf{x}}}_{N+1})$, we aim to update our coefficients towards $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1, k}$. On a high level, for updating the coefficients in the online distributional regression, we do the following steps in the outer cycle: we initialize the predictors ${\eta}_{N+1, k}^{[0, 0]} = g_k\left(\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1, k} \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N, k}\right)$ and this allows us to evaluate $\ensuremath{\boldsymbol{\mathbf{u}}}_{N+1, k}^{[0, 0]}$, ${w}_{N+1, k,k}^{[0, 0]}$ and ${z}_{N+1, k}^{[0, 0]}$. From there, we can update the weighted, exponentially discounted Gramian matrices $\ensuremath{\boldsymbol{\mathbf{G}}}_{N, k} \rightarrow \ensuremath{\boldsymbol{\mathbf{G}}}_{N+1, k}^{[0, 0]}$ and the vectors $\ensuremath{\boldsymbol{\mathbf{h}}}_{N, k} \rightarrow \ensuremath{\boldsymbol{\mathbf{h}}}_{N+1, k}^{[0, 0]}$ {\color{black} by taking

align[align omitted — 592 chars of source]

} and subsequently use an online estimation method from Section (ref) to update the coefficients inside the iteration. As in the batch algorithm (see Algorithm (ref)), we run the outer and inner cycles until convergence.

algorithm[algorithm omitted — 3,550 chars of source]

{\color{black} The following paragraphs discuss some points of the algorithm in detail.} To make the discussion clearer, we call an update step a full update for all distribution parameters $\widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_{N, k} \rightarrow \widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_{N+1, k}$ and associated regression coefficients $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N, k} \rightarrow \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1, k}$ for $k = 1, \ldots, K$, while step might refer to any arbitrary step in the algorithm. {\color{black} We start with a few general points and subsequently discuss the issue of model selection (Section (ref)) and the difference between batch and online estimation (Section (ref)) in detail.}

itemize• Note that in Line 8, in each inner iteration, we start at the Gramian matrix of the previous full fit, not at the previous iterations $r$ or $i$ of the algorithm, since this would imply adding the $\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1, k}$ and $z_{N+1, k}^{[r,i]}$ multiple times to the Gramian matrices within one update step. However, we can (and should) warm-start the OCD algorithm using the coefficient path from the previous iterations within each inner iteration of the update step. • We can potentially use different forget factors for $\gamma_k$ for each distribution parameter. Generally, a higher forget leads to faster adaptation of the coefficients. However, since the estimation of higher moments depends on the estimation of the first moment(s), we note that too aggressive adaptation of the coefficients of the location will lead to an underestimation of the conditional heteroskedasticity and potentially higher moments. • {\color{black} The algorithm is able to process one-step and batch updates of size $b$. In the case of mini-batch updates, all scalar values $u_{N+1, k}^{[r,i]}$, $w_{N+1,k,k}^{[r,i]}$ and $z_{N+1, k}^{[r,i]}$ need to be replaced by their vector counterparts. The Gramian $\ensuremath{\boldsymbol{\mathbf{G}}}_{ N+b,k}$ and $\ensuremath{\boldsymbol{\mathbf{h}}}_{N+b,k}$ can be updated in the same fashion as above, and the inverse $\ensuremath{\boldsymbol{\mathbf{G}}}_{N+b,k}^{-1}$ can be updated using the Woodbury matrix identity. After the update of $\ensuremath{\boldsymbol{\mathbf{G}}}_{N+b,k}$ and $\ensuremath{\boldsymbol{\mathbf{h}}}_{N+b,k}$, the algorithm remains unchanged.} • {\color{black} As a related point, for one-step updates, this can lead to identification issues: Assume we observe a new, somewhat extreme new observation $y_{N+1}$. Whether this new observation warrants an update of the scale, tail, or skewness parameter is unclear. For the presented Algorithm (ref), we observe that the largest updates are generally placed on the distribution parameters that are updated first, and therefore some caution with respect to the ordering should be taken (especially with respect to the parameters for skewness and (symmetric) tail behavior).} • The convergence criteria is the change of the deviance $\mathcal{d} = -2 \sum^N_{n1=} \ell(y_n \mid \hat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_n)$ given the current iteration's fitted values, both for the batch and online case. We track the (exponentially discounted) deviance in the online case. We employ a relative and absolute stopping criteria. • {\color{black} Algorithm (ref) is constant in memory use with respect to the number of observations $N, N+1, N+2, \ldots$ seen. In the non-regularized case, the algorithm needs to store $\ensuremath{\boldsymbol{\mathbf{G}}}_{N,k}$ of size $J \times J$, the vector $\ensuremath{\boldsymbol{\mathbf{h}}}_{N,k}$ of size $J \times 1$ and the coefficient vectors $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N,k}$ of size $J \times 1$ for each distribution parameter $k = 1, \ldots, K$, which gives a total space complexity of $\mathcal{O}(K(J^2+2J))$. In the regularized case, we additionally store the estimation path (for warm-starting) and data for the model selection (see the following paragraph).}

Online Model Selection

{\color{black} As already hinted at in Section (ref), the estimation of the coefficients $\ensuremath{\boldsymbol{\mathbf{\beta}}}_k$ on a path of regularization parameters yields the need for model selection.} Since distributional regression is a likelihood-based approach, we propose the use of information criteria (IC) for this task. We define the generalized information criterion (GIC) as

equation[equation omitted — 166 chars of source]

where $P$ denotes the number of estimated {\color{black} non-zero coefficients}, $\nu = (\nu_0, \nu_1, \nu_2)$ denotes a triplet of parameters, $\tilde{N}_k = \frac{(1 - \gamma_k^N)}{1 - \gamma_l}$ is the effective training length, and $\widehat{\mathcal{L}}$ denotes the maximized value of the likelihood function of the model. The most commonly used information criteria Akaike's Information Criterion {(AIC, $\nu = (2, 0, 0)$)}, Bayesian Information Criterion {(BIC, $\nu = (0, 1, 0)$)} and the Hannan-Quinn Criterion {(HQC, $\nu = (0, 0, 2)$)} can be recovered from the GIC kim2012consistent, kock2016consistent. {\color{black} From here, we have two options to proceed: we can use local model selection on the level of the distribution parameter $k$, based on the residual sum of squares (RSS) of each inner iteration (Line 8 in Algorithm (ref) and Line 9 in Algorithm (ref)) or global model selection, based on the log-likelihood under the full model.

enumerate• Local Model Selection: Under the Gaussian assumption, a local likelihood for each IRLS iteration can be formulated as a function of the residual sum of squares of the local regression of $\ensuremath{\boldsymbol{\mathbf{z_k}}}$ on the design matrix \begin{equation} \log(\widehat{\mathcal{L}}_k) = - \frac{\tilde{N}_k}{2}\log\left(\frac{\sum_{n=1}^N(1-\gamma_k)^{N-n-1}(\ensuremath{\boldsymbol{\mathbf{z}}}_k - \ensuremath{\boldsymbol{\mathbf{X}}}_k\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_k)^2}{\tilde{N}_k}\right) + C_k \end{equation} where $N$ is the number of observations received so far and $C_k = - \frac{\tilde{N}_k}{2} \left(1 + \log(2 \pi) \right)$ is a constant which only depends on the data and hence can be neglected if the data underlying the model selection is the same for all models. This allows an efficient online update of the information criterion for each distribution parameter $\ensuremath{\boldsymbol{\mathbf{\theta}}}_k$ individually. We then choose $$\lambda_k^\ensuremath{\text{opt}} = \mathop{\rm arg\,min}_{\ensuremath{\boldsymbol{\mathbf{\lambda}}}_k} \operatorname{GIC}(\widehat{\mathcal{L}_k}(\lambda_k), \nu_k).$$ • Global Model Selection: Here, we employ the log-likelihood of the distribution $\mathcal{F}$ directly in the GIC and select the optimal regularization parameter for each by taking: $$\lambda_k^\ensuremath{\text{opt}} = \mathop{\rm arg\,min}_{\ensuremath{\boldsymbol{\mathbf{\lambda}}}_k} \operatorname{GIC}(\widehat{\mathcal{L}}(\lambda_k), \nu).$$

The local model selection is used, e.g., for the selection of regularization parameters in the backfitting step in the classic GAMLSS algorithm and in the R gamlss.lasso package ziel2021gamlss. It has the advantage of avoiding expensive likelihood evaluations, but relies on the assumptions that the local fit's residuals are approximately Gaussian distributed. Furthermore, the local model selection allows specifying different parameterizations of the GIC for different distribution parameters.} We suggest selecting a {\color{black} more conservative parameterization of the GIC} for higher $k$ to avoid overfitting in modeling the conditional scale, kurtosis, and skewness {\color{black} when using local model selection. We also note that since the coefficients and estimates for different distribution parameters live on very different scales, it can be the case that $\varepsilon_{k, \lambda}$ needs to be smaller, otherwise the lower end of the coefficient path might not approach the OLS solution.} groll2019lasso analyze the impact of different shrinkage parameters on batch distributional regression. Also marcjasz2023distributional see the advantages of different regularization for the distribution parameters when using distributional neural networks for time series forecasting.

Comparison of Online and Batch Estimation

An important relationship to the batch version of our proposed algorithm becomes apparent when comparing how weights are updated. In the batch setting, with repeated fits of the {\color{black} model using all observations}, the weight matrix $\ensuremath{\boldsymbol{\mathbf{W}}}_{k,k}^{[r,i]} = \operatorname{Diag}\left(w_{0,k,k},\ldots, w_{N,k,k}\right)$ is updated in every iteration $r,i$ and again, in the next batch fit for $N+1$, the weight matrix is updated for all $w_{n,k,k}$ from $n = 1,\ldots, N+1$ in all iterations. In the online setting, $\ensuremath{\boldsymbol{\mathbf{W}}}_{N,k,k}$ is the weight matrix after convergence of the update step $N-1 \rightarrow N$. In the update step $N \rightarrow N+1$, we cannot update $\ensuremath{\boldsymbol{\mathbf{W}}}_{N,k,k}$ anymore (see also the Assumptions noted in Section (ref)). Therefore, we set

equation[equation omitted — 304 chars of source]

and can only update $w_{N+1, k,k}^{[r,i]}$. Note that due to the update rule in Equation (ref), this effect is counter-weighted by the exponential forget $\gamma$. However, there is a delicate balance to strike to trade off the beneficial effect of the forget and increased instability in the coefficient estimation {\color{black} due to smaller effective sample sizes}. This might lead to slower convergence compared to batch learning if the data is drawn from a stationary process. The online distributional regression is, therefore, an approximation of the repeated batch estimation, contrary to the case for, e.g., RLS or online coordinate descent in the case of regression to the mean, where the update step leads to the equivalent results. When choosing a specific distribution assumption for $\mathcal{F}$ in the IRLS framework, the weight matrix for the {\color{black} distribution parameter} $\theta_k$ in each update step can be recovered with the help of Equations (ref) to (ref). The online weight matrix is retrieved by iteratively inserting the weight matrix in Equation (ref). For illustration, building on the standard linear-Gaussian case for the distributional regression in Remark (ref), we can examine the difference explicitly for the mean equation: {\color{black}

remarkAssume a linear-Gaussian distributional regression with $\theta_1 = \mu$ and $\theta_2 = \sigma^2$ of Remark (ref). Assume that we have seen $n = 1, \ldots., N, N+1, \ldots, M$ observations, where $N$ denotes the initial fit and the steps $N+1, \ldots, M$ are updated online. In the batch and online setting, the weight matrices for $\theta_1=\mu$ for the step $(M-1) \rightarrow M$ are: \begin{small} \begin{align*} \ensuremath{\boldsymbol{\mathbf{W}}}_{M,1,1}^{[r,i], \operatorname{online}} &= \operatorname{diag}\Biggl( \underbrace{ \dfrac{(1- \gamma)^{M-1}}{\widehat{\sigma}^2_{1|N}}, \ldots, \dfrac{(1-\gamma)^{M-N-1}}{\widehat{\sigma}^2_{N|N}} }_\normalfont Initial batch $n=1, \ldots, N$., \underbrace{ \dfrac{(1-\gamma)^{M-N-2}}{\widehat{\sigma}^2_{N+1|N+1}} , \ldots, \dfrac{(1-\gamma)}{\widehat{\sigma}^2_{M-1|M-1}}, \dfrac{1}{(\widehat{\sigma}^2_{M|M})^{[r, i]}} }_\normalfont Online updates $n = N+1, \ldots, N+M$. \Biggl) \\ \ensuremath{\boldsymbol{\mathbf{W}}}_{M,1,1}^{[r,i], \operatorname{batch}} &= \operatorname{diag}\left( \dfrac{(1 - \gamma)^{M-1}}{(\widehat{\sigma}^2_{1|M})^{[r,i]}}, \ldots, \dfrac{(1- \gamma)^{M-N-1}}{(\widehat{\sigma}^2_{N|M})^{[r, i]}}, \ldots, \dfrac{(1 - \gamma)}{(\widehat{\sigma}^2_{M-1|M})^{[r, i]}} , \dfrac{1}{(\widehat{\sigma}^2_{M|M})^{[r,i]}} \right) \end{align*} \end{small} where $\gamma$ is the exponential discounting factor, which is applied to the repeated batch estimation and the online estimation, and the subscript $n|i$ denotes the last observation $i$ available. The difference can be summarized in three terms \begin{equation} \begin{split} \ensuremath{\boldsymbol{\mathbf{W}}}_{M,1,1}^{[r, i], \operatorname{online}} - \ensuremath{\boldsymbol{\mathbf{W}}}_{M,1,1}^{[r, i], \operatorname{batch}} = \underbrace{\sum_{n=1}^{N} (1-\gamma)^{M-n} \left(\dfrac{1}{\widehat{\sigma}^{2,\operatorname{batch}}_{n|N}} -\dfrac{1}{(\widehat{\sigma}^{2,\operatorname{online}}_{n|N})^{[r, i]}}\right)}_\normalfont Initial batch. + \\ \underbrace{\sum_{n=N+1}^{M-1} (1-\gamma)^{M-n} \left(\dfrac{1}{\widehat{\sigma}^{2,\operatorname{batch}}_{n|n}} -\dfrac{1}{(\widehat{\sigma}^{2,\operatorname{online}}_{n|M})^{[r,i]}}\right)}_\normalfont Online fit $N+1, \ldots, M-1$ + \underbrace{\dfrac{1}{(\widehat{\sigma}^{2,\operatorname{online}}_{M|M})^{[r,i]}} - \dfrac{1}{(\widehat{\sigma}^{2,\operatorname{batch}}_{M|M})^{[r,i]}}}_\normalfont Current step $M-1 \rightarrow M$. \end{split} \end{equation} where it is important to note that even in the last term $(\widehat{\sigma}^{2,\operatorname{online}}_{M|M})^{[r,i]}$ and $(\widehat{\sigma}^{2,\operatorname{batch}}_{M|M})^{[r,i]}$ from the online and batch estimation are not identical, since the online estimation is path-dependent (see the recursive insertion of $w_{N+1, k,k}$ in Equation (ref)).

}

From Remark (ref) we see that the difference in the online and batch setting depends on the length of the initial batch {\color{black} $N$, how long the algorithm is running online $M-N$, and the strength of the exponential discounting. Looking at Equation (ref) and keeping in mind that for large sample sizes, $(1-\gamma)^{M}, \ldots, (1-\gamma)^{M-n}$ will eventually approach values close to zero, it is tempting to conclude that for sufficiently high discounting and after a sufficiently large number of online updates $M$ the first term of Equation (ref) ceases to influence the estimation.\footnote{Precisely, when $(1-\gamma)^{(N-M)} \rightarrow 0$ the observations that are in the initial sample receive a weight of zero.} However, due to the path-dependency of the online algorithm, where each update step's weight depends on the previous update step, this argument does not hold. For different combinations of the distribution $\mathcal{F}$, link functions $g_k(\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}})$, and distribution parameter $k$ the exact form of the above weights might be significantly more complex. This result shows a main difference to the regression to the mean} setting, where the batch and online results coincide {\color{black} for OCD and RLS algorithms}.

\FloatBarrier

Simulation Study

{\color{black} This section presents three distinct simulation studies on different aspects of the online distributional regression model. We analyze timings, the impact of the length of the initial training set, and the forgetting factor. Here, the data-generating process (DGP) is similar for all three studies. We consider a setting where $\ensuremath{\boldsymbol{\mathbf{y}}}$ is normally distributed. The location and scale parameters depend on the $N \times J$ covariate matrix $\ensuremath{\boldsymbol{\mathbf{X}}} =\ensuremath{\boldsymbol{\mathbf{X}}}_{1} = \ensuremath{\boldsymbol{\mathbf{X}}}_{2}$. We simulate the entries of $\ensuremath{\boldsymbol{\mathbf{X}}}$ by a $J$-dimensional normal distribution ${\ensuremath{\boldsymbol{\mathbf{X}}}_{i} \sim \mathcal{N}( 0_J, I_J/2)}$ and define

equation[equation omitted — 240 chars of source]

with $g(x, \alpha) = \log\left(1 + e^{x}\right) + \alpha$, a shifted softplus link function. Accordingly, we use a shifted softplus link function for estimating the scale parameter of the normal distribution. We set $\alpha=0.5$ to ensure sufficiently large variances. This greatly improves the stability of the estimations. The coefficient vectors $\ensuremath{\boldsymbol{\mathbf{\beta}}}_k$ of length $J$ are drawn from the uniform distribution:

align[align omitted — 86 chars of source]

Unless mentioned otherwise, we set $J=10$, use $N=10000$ observations, use $2000$ observations for the initial fit, and execute 10 simulation runs.

enumerate• The first simulation examines the convergence of the parameters to their true values in a batch setting. For that, we compute the $\ell_1$-norm $\| \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}} \|_1$ of the difference between $\ensuremath{\boldsymbol{\mathbf{\beta}}}_k$ and $\ensuremath{\boldsymbol{\mathbf{\widehat{\beta}}}}_k$ for different initial sample sizes. Thereby, we consider sample sizes of an exponential grid running from $2^7$ to $2^{13}$. We consider model estimations using RLS and OCD. Figure (ref) illustrates the results. The results are as expected. Both parameters converge to their true values as $N$ grows. There are no notable differences between RLS and OCD in this setting. • The second study compares the estimation times between RLS and OCD. Thereby, we distinguish between the initial fit and the subsequent update steps. Furthermore, we vary the number of features $J$. We consider an exponential grid ranging from $J=2^2$ to $J=2^8$. RLS has a per-iteration time complexity of $\mathcal{O}(JN^2)$. Online coordinate descent has a complexity of $\mathcal{O}(JN)$, which means we expect RLS computation times to be lower for smaller feature sets. On the other hand, we expect OCD to be faster if the feature set is sufficiently large. Our measurements are depicted in Figure (ref). Surprisingly, RLS and OCD need more time for the initial fit of the smallest model ($J=4$) compared to the larger ones. We barely observe differences in the timings for smaller data sets. For more than $J=2^6=64$ features, online updates using OCD are faster compared to RLS. • Additionally, we conduct a third study. This demonstrates how forgetting affects the speed of parameter adjustment. For this, we slightly adjust the data-generating process. We define change points at which $\theta_{j}$ changes. These change points are drawn from a geometric distribution with a mean of $1000$. Then, we estimate distributional models using LASSO, and we consider a grid of forget values $\gamma = \{0, 2^{-11}, 2^{-10}, \ldots, 2^{-3}\}$. Figure (ref) presents the results. In this setting, a positive forget is beneficial. More precisely, the forget value of $2^{-8}$ produces the best results across all considered scoring rules. This forget corresponds to an effective sample size of $2^8 = 256$. Furthermore, we observe an increase in estimation time as the forgetting increases. This is expected. As the forgetting increases, the effective sample size decreases, and the algorithm becomes less stable. This means more iterations are needed for the algorithm to converge. Figure (ref) in Appendix (ref) presents the parameter estimates for selected forget values for all 20 parameters over time.

}

figure[figure omitted — 1,043 chars of source]

figure[figure omitted — 613 chars of source]

figure[figure omitted — 792 chars of source]

\FloatBarrier

Forecasting Study for Electricity Prices

Data, Study Design and Models

The following section presents an application in electricity price forecasting (EPF). We employ the same setting as marcjasz2023distributional and brusaferri2024line in forecasting day-ahead electricity prices for the German short-term electricity market. We directly compare the forecasting performance of online estimation vs repeated batch estimation directly for models using the Gaussian, Student-$t$, and Johnson $S_U$ (JSU) distribution and for OLS/RLS and LASSO estimation (see Table (ref)). Furthermore, we compare our approach to established benchmark models from the literature and conduct several ablation studies (see Appendix (ref), Tables (ref) and (ref)). The following paragraphs describe the German day-ahead electricity market, data set, and models in detail.

figure[figure omitted — 880 chars of source]

Let $h \in \{0, 1, \ldots, H\}$ and $H=23$ denote the 24 delivery hours and $d$ denote the delivery day. The day-ahead electricity market serves as the main trading venue for electricity and as reference price for the long-term futures market. It is organized as a daily auction at $d-1$, 12:00 hours for all 24 delivery hours of the following day, {\color{black} where electricity generators, consumers, and traders can place price-volume bids to buy and sell electricity in each delivery hour. The market is cleared in a pay-as-cleared auction, taking into account the available transfer capacity between the European countries} viehmann2017state, marcjasz2023distributional, ziel2018day, lago2018forecasting, epftoolbox. {\color{black} Figure (ref) gives an overview of the electricity price time series, the residual load (system load minus renewable energy production), and fuel prices for conventional generation assets. The residual load is the amount of electricity production that has to be covered by conventional generation sources (e.g., gas, coal, and oil). From an economic point of view, the renewable generation is often assumed to operate at zero or even negative marginal costs, while conventional power plants need to cover fuel and emission certificate costs, and therefore, the residual load is an important factor in the price formation.} The data is retrieved from marcjasz2023distributional\footnote{See \url{https://github.com/gmarcjasz/distributionalnn}. marcjasz2023distributional provide data from the European Network of Transmission System Operators (ENTSO-E) and Refinitiv.} The initial training set comprises 2015-01-15 to 2019-06-26 (1631 days), and the test set comprises 2019-06-26 to 2020-12-31 (534 days). This is in line with the test and training split of brusaferri2024line and marcjasz2023distributional.

We employ a so-called expert-type model as it is common in electricity price forecasting lago2018forecasting, ziel2018day, marcjasz2023distributional. The model consists {\color{black} of 43 linear terms and of 5 non-linear additive effects modeled by B-Spline bases. The covariates capture} the autoregressive price effects, seasonal effects, the fundamental effects of renewable generation, and the influence of fuel prices. We employ the {\color{black} Gaussian}, Student-$t$, and Johnson's $S_U$ (JSU) distributions. {\color{black} The Gaussian distribution serves as a benchmark, and the latter two} have been used for electricity price forecasting and other financial applications. In the distributional framework, we model each distribution parameter $k = 1, \ldots, p$ for each delivery hour $h = \{0, 1, \ldots, H\}$ individually:

align[align omitted — 907 chars of source]

where $\beta_{k, 0, h}$ is the intercept or bias, $\beta_{k, 1, h}$ to $\beta_{k, 30, h}$ capture autoregressive effects for the previous day $d-1$, all delivery hours, and for the last 7 days of the same delivery hour $h$. $\beta_{k, 31, h}$ to $\beta_{k, 34, h}$ capture the influence of very large and low prices for the previous day. $\beta_{k, 35, h}$ models the influence of the daily residual load. $\beta_{k, 36, h}$ to $\beta_{k, 42, h}$ capture the effects of weekly seasonality and holidays.\footnote{German public holidays are taken from the Python holidays package. On a public holiday, we set all other weekday dummies to zero. The holiday days include New Year's, Good Friday, Easter Monday, Labor Day, Ascension Day, Whit Monday, German Unity Day, and Christmas.} Lastly, the terms $b_{k,1,h}(\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}})$ to $b_{k,5}$ capture the non-linear effects of hourly residual load, European Emission Allowances (EUAs), natural gas, coal and oil prices by B-Splines of degree 2 and using 4 knots, placed according to the quantiles of the covariates in the training set eilers1996flexible.

We benchmark forecasts using established, strictly proper probabilistic scoring rules gneiting2008probabilistic, gneiting2011making, gneiting2014probabilistic, nowotarski2018recent. For the mean prediction, we employ the Root Mean Squared Error (RMSE), and {\color{black} for the median prediction, we use the} Mean Absolute Error (MAE). We evaluate the Coverage (COV) and the Interval Score bracher2021evaluating for the {\color{black} 50% and 80%} prediction intervals (PI). For the full predictive distribution, we evaluate the Log Score (LS) and the continuous ranked probability score (CRPS) using the approximation via the Pinball Score (PS) on a dense grid of quantiles $\mathcal{Q} = \{0.01, 0.02, 0.03, \ldots, 0.99\}$. The implementation of the scoring rules is provided by the scoringrules package zanetta_scoringrules_2024. We evaluate the statistical significance of the difference in predictive accuracy using the Diebold-Mariano test diebold2002comparing, diebold2015comparing. {\color{black} For the coverage, we evaluate the Kupiec-test (see kupiec1995techniques) per hour $h$. Detailed definitions of the scoring rules and the statistical tests can be found in Appendix (ref)}. Lastly, we evaluate the computation time for all forecasting studies.

Results and Discussion

The following section presents and discusses the results of the forecasting study. Figure (ref) gives an exemplary forecast. The time-varying width of the prediction intervals illustrates the importance of modelling the conditional volatility. Figure (ref) shows the non-linear effects of residual load on the scale parameter (A similar Figure for the location $\mu$ can be found in Appendix (ref)). The effect is $u$-shaped, indicating that the price volatility increases for very high and very low residual load levels. This is in line with economic theory, as extreme residual load levels indicate scarcity or oversupply situations, which lead to more volatile prices.

figure[figure omitted — 1,082 chars of source]

figure[figure omitted — 1,280 chars of source]

table[table omitted — 27,770 chars of source]

Table (ref) presents our main results, the direct comparison of repeated batch and online estimation for the distributional regression model. We make the following observations about the predictive performance in our forecasting study:

itemize• Comparing the distributional regression models in Table (ref) like-for-like, the predictive performance of the repeated batch estimation and online estimation is very close and, for various models and scoring rules, not statistically significantly different. At the same time, the computational effort of the online estimation is 2-3 orders of magnitude lower than that of the repeated batch estimation. This result underscores the robustness and efficiency of the proposed algorithm. • We see a general improvement in forecasting performance when using regularized estimation approaches compared to ordinary least squares. Especially for the JSU distribution, the classic OLS estimation struggles with the heavy-tailed and skewed electricity price data in the batch and online estimation. The issue is more pronounced in the online case, but can be successfully alleviated by employing regularized estimation methods such as the LASSO or elastic net (see also Table (ref), Panel A). The LASSO-estimated JSU model is the strongest in our study.

Furthermore, we make the following observations when comparing to established benchmarks and conducting ablation studies:

itemize• In the comparison to established time-series benchmarks and state-of-the-art, neural-network-based EPF approaches, we see the competitive performance of the online distributional regression models. We see that the AutoARIMA generally struggles to adapt to the highly volatile data, and our investigation showed that stepwise model selection seems not robust enough. • Compared with state-of-the art results of brusaferri2024line, our approach delivers slightly better performance in terms of the CRPS, but slightly worse coverage. The DDNN approach of marcjasz2023distributional yields better results than the distributional regression approaches and the conformalized DDNNs of brusaferri2024line. However, the DDNN approaches necessitate specialized computing infrastructure, and even without hyperparameter tuning, need more than 100 times the computing time of the online distributional regression on a standard desktop PC. • In the comparison between the linear online distributional regression and the model with B-Spline bases for the fundamental covariates, we see that the linear model yields worse performance for the Normal and Student-$t$ distribution, but slightly improved performance for the JSU distribution. We attribute this to an increased risk of overfitting and the issue of identifiability in distributional regression, especially for the skewness and tail behavior in the non-linear model, as already discussed in Section (ref). • Table (ref), Panel B gives results for the estimation using different forget factors, using an effective sample size of 1 to 4 years in the exponential discounting. The intuitive rationale behind exponential forgetting, that newer observations contribute more to forecasts, is confirmed in the results. For all three distributional assumptions, the forecasting can be improved using a small forget. However, we also note that increasing $\gamma$ yields higher computational costs to a small extent. • Further discussion of the ablation studies can be found in Appendix (ref). In summary, we note that the local model selection approach is preferable due to its computational efficiency, and that batch updates can be used to speed up the estimation further, at the cost of some predictive accuracy.

Wrapping up, our study demonstrates that online distributional regression models achieve competitive predictive performance compared to both traditional time-series benchmarks and advanced neural network-based approaches. Moreover, they offer significantly lower computational costs - often by two to three orders of magnitude - making them highly efficient for real-time forecasting.

Implementation in the Python Package ondil

We provide an open-source implementation of our algorithm in the Python package ondil. To the best of the authors' knowledge, this package is the first native Python implementation of IRLS-based distributional regression and can, therefore, provide a basis for future extensions.\footnote{The package pyNM implements GAMLSS for Python as binding to the R library \url{https://github.com/ppsp-team/PyNM}.} We provide the following features:

itemize• Our package is written in an online-first fashion and provides an API compatible with scikit-learn pedregosa2011scikit, the major Python machine learning package. Beyond the estimator.fit(X, y) and estimator.predict(X, y) methods, our estimator classes provide the estimator.update(X, y) method to allow for incremental updating. • The implementation relies on few dependencies 2020NumPy-Array, 2020SciPy-NMeth, lam2015numba. We employ just-in-time compilation using \texttt{numba} to achieve a high-performance implementation and employ various computational tricks such as active-set coordinate descent, different warm-starting options (using the previous fit or a mixture of previous fit and previous $\beta_\lambda$ on the same coefficient path), and allow for random selection during coordinate descent to speed up the convergence of the coordinate descent shi2016primer, wright2015coordinate. • We implement standard distributions and facilitate extensions based on \texttt{scipy}-distributions.

The code for our implementation is open source {\color{black} on GitHub} and the package is available at the Python Package Index.\footnote{See: \url{https://pypi.org/project/ondil} and \url{https://github.com/simon-hirsch/ondil}.} {\color{black} All studies in this paper have been run using ondil version 0.3.2.}

Discussion and Conclusion

{\color{black} This paper presents an efficient, scalable approach to distributional regression by the implementation of an online estimation algorithm} for the well-known, parametric GAMLSS rigby2005generalized, klein2024distributional. {\color{black} We discuss regularized estimation using ridge, LASSO, and the elastic net, and the associated issue of online model selection. Furthermore, we formally analyze the approximation quality of online estimation and batch estimation.} Lastly, we provide an open-source Python implementation in the ondil package. We validate our approach in a forecasting study for the German day-ahead electricity prices, a highly volatile data set. The online distributional regression delivers competitive forecasting accuracy in the CRPS compared to repeated batch estimation and other distributional forecasting approaches while reducing the computation time by some 2-3 orders of magnitude. Furthermore, we provide simulation results on the impact of the initial training size and the exponential discounting on the approximation quality of the online models compared to repeated batch fitting. Our simulation results are in line with the aforementioned theoretical results.

Our research opens up multiple avenues for future research. First, while our empirical results are promising, further theoretical results on the error bounds of the online algorithm would further increase the trust in the proposed method. Secondly, we present a regularized method for the online estimation, but the (online) model selection in distributional regression is a relatively untapped field, and advances here will directly benefit practical applications. {\color{black} While our approach is capable of the estimation of non-linear models using linear basis functions, the implementation of regularized splines and other additive effects is a practically relevant addition. Lastly, online learning for distributional models is naturally confronted with issues of identifiability, and our approach is no exception. We note that regularization plays a critical role in controlling the issue in a practical study; a rigorous formal analysis seems like a promising avenue for further research, benefiting all research on online distributional models.}

acksSimon Hirsch is employed as an industrial PhD student by Statkraft Trading GmbH and gratefully acknowledges the support and funding received. Simon Hirsch is grateful to Daniel Gruhlke for many helpful discussions. This work contains the author's opinions and does not necessarily reflect Statkraft position. The authors declare no conflict of interest.