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.
134,482 characters
Online Distributional Regression
\date{January 14, 2026}
\title{Online Distributional Regression}
\author{Simon Hirsch}
\orcid{0009-0008-1409-9677}
\email{[email removed]}
\email{[email removed]}
\affiliation{
\institution{Statkraft Trading GmbH}
\country{Germany}
}
\affiliation{
\institution{Data Science in Energy and Environment, University of Duisburg-Essen}
\country{Germany}
}
\author{Jonathan Berrisch}
\orcid{0000-0002-4944-9074}
\email{[email removed]}
\author{Florian Ziel}
\email{[email removed]}
\orcid{0000-0002-2974-2660}
\affiliation{
\institution{Data Science in Energy and Environment, University of Duisburg-Essen}
\country{Germany}
}
\keywords{{online learning, time series, GAMLSS, LASSO, electricity price forecasting (EPF)}}
\date{\today}
\begin{abstract}
Large-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 \texttt{Python} package {\color{black} \texttt{ondil}}.
\end{abstract}
\maketitle
\textbf{Working Paper -- This Version:} January 14, 2026
\section{Introduction}\label{sec:introduction}
Large-scale streaming data are common in modern applications of machine learning and have led to the development of online learning algorithms \citep{cesa2021online}. For processes driven by a high-dimensional covariate space, regularized algorithms have been presented by, e.g., \cite{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 \citep[for a review on distributional regression see, e.g. \cite{kneib2023rage} and \cite{klein2024distributional}, for probabilistic forecasting see][]{gneiting2014probabilistic, nowotarski2018recent, alvarez2021probabilistic, petropoulos2022forecasting, ziel2022m5}. However, online learning approaches for distributional regression remain sparse in the literature and can be grouped as:
\begin{enumerate}
\item Adaptive tracking of the {\color{black} variance-controlling} scale parameter \citep[see e.g.][]{alvarez2021probabilistic, vilmarest2024viking},
\item Adaptive estimation of the conditional heteroskedasticity \cite[see e.g.][]{priouret2005recursive, dahlhaus2007recursive, hendrych2018self, cipra2018robust,werge2022adavol, wintenberger2024stochastic},
\item Adaptive conformal prediction approaches \citep[see e.g.][]{zaffran2022adaptive, bhatnagar2023improved, gibbs2024conformal, dutot2024adaptive, brusaferri2024line},
\end{enumerate}
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 \textit{online coordinate descent} (OCD) algorithm introduced by \cite{angelosante2010online} and \cite{messner2019online} and the \textit{generalized additive models for location, scale and shape} (GAMLSS) introduced by \cite{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 } \cite{rigby2005generalized} we define the parametric GAMLSS as follows:
\begin{equation}\label{eq:y_distributed_f_theta}
{\color{black} y_n \sim
\mathcal{F}\left(\theta_{n,1} ,\ldots, \theta_{n,K}\right) = \mathcal{F}\left(\ensuremath{\boldsymbol{\mathbf{\theta}}}_{n, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}}\right) }
\end{equation}
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}
\begin{equation}
g_k(\ensuremath{\boldsymbol{\mathbf{\theta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}, k}) = \ensuremath{\boldsymbol{\mathbf{\eta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}, k} = \beta_0 +
\underbrace{ \sum_{j\in \mathcal{J}_k} {\beta_{k,j} \ensuremath{\boldsymbol{\mathbf{x}}}_{j}}}_{\text{linear part}} +
\underbrace{ \sum_{j\in\mathcal{B}_k} b_{k,j}(\ensuremath{\boldsymbol{\mathbf{x}}}_{j})}_{
\text{spline part}}
\label{eq:distreg_model}
\end{equation}
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 \citep[see e.g.][]{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{eq:distreg_model} as
\begin{equation}
g_k(\ensuremath{\boldsymbol{\mathbf{\theta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}, k}) = \ensuremath{\boldsymbol{\mathbf{\eta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}, k} = \ensuremath{\boldsymbol{\mathbf{X}}}_k\ensuremath{\boldsymbol{\mathbf{\beta}}}_k
\end{equation}
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} \cite{rigby2005generalized} introduced estimation via iteratively reweighted least squares (IRLS), which has been extended to LASSO-type penalties by \cite{groll2019lasso}. Further regularized estimation approaches have been proposed by \cite{ziel2021gamlss}.}
{\color{black}
Alternatively, estimation via Bayesian approaches (direct and via Bayesian-IRLS), direct likelihood minimization using stochastic gradient descent \citep[SGD, see e.g.][]{thielmann2024neural,umlauf2024scalable} and gradient-boosted decision trees \citep[see e.g.][]{maerz2019xgboostlss, cevid2022distributional} are popular.
}
\begin{figure}[htb!]
\Description[Building blocks of the algorithm]{The paper combines two strands of literature by merging online estimation of linear models and iteratively reweighted least squares estimation for distributional regression models.}
\centering
\includegraphics[width=0.8\textwidth]{figures/excalidraw_paper_flow_tight.png}
\caption{The building blocks of our algorithm give the flow of the paper. Section \ref{sec:batch_gamlss} introduces batch distributional regression and estimation using IRLS, Section \ref{sec:online_regression} introduces incremental estimation methods for linear models. Finally, Section \ref{sec:method} combines both to yield the online distributional regression model.}
\label{fig:paper_flow}
\end{figure}
{\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 \cite{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{fig:paper_flow} gives an idea of the building blocks of our algorithm. In summary, we make the following contributions to the current literature:
\begin{enumerate}
\item We propose {\color{black} an \emph{online} and \emph{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{sec:method}).
\item 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{sec:simulation}).
\item 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{sec:forecasting}).
\item We provide an open-source, ready-to-use \texttt{Python} implementation of the online distributional regression model {\color{black} in the package \texttt{ondil}}.
Our implementation is based on \texttt{numpy} and \texttt{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{sec:python})}.
\end{enumerate}
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 \citep[as implemented in the batch case by][]{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{fig:paper_flow}.} The following two sections give the building blocks of our algorithm: First, Section \ref{sec:batch_gamlss} briefly reviews the algorithm for the parametric GAMLSS. Afterwards, Section \ref{sec:online_regression} reviews online estimation of regularized linear models. Section \ref{sec:method} presents our main contribution, the online distributional regression. The following Sections \ref{sec:simulation} and \ref{sec:forecasting} present a simulation study on the properties of the online algorithm and our real-world application for energy markets. Section \ref{sec:python} presents our open-source implementation. Finally, Section \ref{sec:conclusion} discusses our results and concludes the paper.}paper.}
{\color{black}
\section{Batch Distributional Regression} \label{sec:batch_gamlss}}
{\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 \citet[Appendix B]{rigby2005generalized} and \cite{stasinopoulos2008generalized}. The full algorithm} is outlined in Algorithm \ref{alg:batch_gamlss}. {\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{eq:y_distributed_f_theta} and \ref{eq:distreg_model}:
\begin{equation*}
y_n \sim \mathcal{F}(\ensuremath{\boldsymbol{\mathbf{\theta}}}_{n, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}}) \quad \text{and }\quad
g_k(\ensuremath{\boldsymbol{\mathbf{\theta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}, k}) = \ensuremath{\boldsymbol{\mathbf{\eta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}, k} = \ensuremath{\boldsymbol{\mathbf{X}}}_k\ensuremath{\boldsymbol{\mathbf{\beta}}}_k
\end{equation*}
}
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
\begin{equation}\label{eq:gamlss_vector_u}
\ensuremath{\boldsymbol{\mathbf{u}}}_k
= \dfrac{\partial \ensuremath{\boldsymbol{\mathbf{\ell}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k}
= \dfrac{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k}
\odot \left(
\dfrac{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k}
\right)^{-1}
\end{equation}
}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
\begin{equation}\label{eq:gamlss_score}
\ensuremath{\boldsymbol{\mathbf{z}}}_k = \ensuremath{\boldsymbol{\mathbf{\eta}}}_k + \ensuremath{\boldsymbol{\mathbf{W}}}_{k,k}^{-1}\ensuremath{\boldsymbol{\mathbf{u}}}_k
\end{equation}
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}
\begin{align}
\ensuremath{\boldsymbol{\mathbf{W}}}^{\text{NR}}_{k,s} &= - \text{diag} \left(\dfrac{
\partial^2 \ensuremath{\boldsymbol{\mathbf{\ell}}}
}{
\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k \partial \ensuremath{\ensuremath{\boldsymbol{\mathbf{\eta}}}^\top}_s
}\right), \quad
\ensuremath{\boldsymbol{\mathbf{W}}}^{\text{FS}}_{k,s} &= - \text{diag} \left(\mathbb{E}\left[\dfrac{
\partial^2 \ensuremath{\boldsymbol{\mathbf{\ell}}}
}{
\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k \partial \ensuremath{\ensuremath{\boldsymbol{\mathbf{\eta}}}^\top}_s
}\right] \right), \quad
\ensuremath{\boldsymbol{\mathbf{W}}}^{\text{qNR}}_{k,s} &= - \text{diag}\left( \dfrac{
\partial \ensuremath{\boldsymbol{\mathbf{\ell}}}
}{
\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k
} \ensuremath{\odot}
\dfrac{
\partial \ensuremath{\boldsymbol{\mathbf{\ell}}}
}{
\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_s
}
\right) \label{eq:gamlss_weight}
\end{align}
giving the observed (Newton-Raphson), expected (Fisher-Scoring) or product (quasi Newton-Raphson) score function \citep[see Appendix B of][]{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{eq:gamlss_weight}
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:
\begin{equation*}
\mathbb{E}\left[
\frac{\partial^2 \ensuremath{\boldsymbol{\mathbf{\ell}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k^2}
\right]
= \mathbb{E}\left[
\frac{\partial \ensuremath{\boldsymbol{\mathbf{\ell}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k} \odot
\frac{\partial \ensuremath{\boldsymbol{\mathbf{\ell}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k}
\right]
= \mathbb{E}\left[
\frac{\partial \ensuremath{\boldsymbol{\mathbf{\ell}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k} \odot
\frac{\partial \ensuremath{\boldsymbol{\mathbf{\ell}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k}
\right]
\ensuremath{\oslash} \left(\frac{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k}\right)^{2}.
\end{equation*}
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:
\begin{align}\label{eq:gamlss_working_vector}
\ensuremath{\boldsymbol{\mathbf{z}}}_k
&= \ensuremath{\boldsymbol{\mathbf{\eta}}}_k + \ensuremath{\boldsymbol{\mathbf{W}}}_{k,k}^{-1} \left(\dfrac{\partial
\ensuremath{\boldsymbol{\mathbf{\eta}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k}
\odot \left(
\dfrac{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k}
\right)^{-1}\right) \\
&= \ensuremath{\boldsymbol{\mathbf{\eta}}}_k +\left(
\mathbb{E}\left[
\frac{\partial \ensuremath{\boldsymbol{\mathbf{\ell}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k} \odot
\frac{\partial \ensuremath{\boldsymbol{\mathbf{\ell}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k}
\right]
\ensuremath{\oslash} \left(\frac{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k}\right)^{2}\right)^{-1} \ensuremath{\odot}
\left(\dfrac{\partial
\ensuremath{\boldsymbol{\mathbf{\eta}}}}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k} \odot \left(\dfrac{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k}
\right)^{-1}
\right),
\end{align}
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:
\begin{equation}
\dfrac{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k}{\partial \ensuremath{\boldsymbol{\mathbf{\theta}}}_k}
= \left(\dfrac{\partial g_k^{-1}(\ensuremath{\boldsymbol{\mathbf{\eta}}}_k)}{\partial \ensuremath{\boldsymbol{\mathbf{\eta}}}_k}\right)^{-1}
\end{equation}
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)$ \citep{rigby2005generalized, stasinopoulos2008generalized, stasinopoulos2024generalized}, however, \cite{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.
\begin{algorithm}[h]
\caption{{\color{black} Iteratively reweighted least squares (IRLS) algorithm for the estimation of distributional regression models, see \cite[][RS-algorithm]{rigby2005generalized}}.}
\label{alg:batch_gamlss}
\DontPrintSemicolon
\KwIn{\text{Distribution $\mathcal{F}$, response variable $\ensuremath{\boldsymbol{\mathbf{y}}}$ and covariates $\ensuremath{\boldsymbol{\mathbf{X}}}$}}
Initialize the fitted values $\widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_k^{[0, 0]}$ for $k = 1, \ldots, K$. \;
Evaluate the linear predictors ${\widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}^{[0, 0]}_k = g_k(\widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_k^{[0, 0]})}$ for $k = 1, \ldots, K$. \;
\For{$r = 1, \ldots, \text{\normalfont until convergence}$}{
\For{$k = 1, \ldots, K$}{
Start the inner cycle. \;
\For{$i = 0, 1, \ldots$}{
Evaluate $\ensuremath{\boldsymbol{\mathbf{u}}}_k^{[r,i]}$, $\ensuremath{\boldsymbol{\mathbf{W}}}_{k,k}^{[r,i]}$ and $\ensuremath{\boldsymbol{\mathbf{z}}}_k^{[r,i]}$. \;
Regress $\ensuremath{\boldsymbol{\mathbf{z}}}_k^{[r,i]}$ on $\ensuremath{\boldsymbol{\mathbf{X}}}$ using weights $\ensuremath{\boldsymbol{\mathbf{W}}}_{k,k}^{[r,i]}$ to obtain $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_k^{[r,i+1]}$. \;
Calculate the updated $\widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}_k^{[r,i+1]}$ and $\widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_k^{[r,i+1]}$. \;
Evaluate the convergence.\;
}
{\color{black} \textbf{break} if $| \ell(\ensuremath{\boldsymbol{\mathbf{y}}} \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}^{[r, i+1]}) - \ell(\ensuremath{\boldsymbol{\mathbf{y}}} \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}^{[r, i]} )|$ is sufficiently small.} \;
Set
$\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_k^{[r + 1, i]} \gets \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_k^{[r, i]}$, set
$\widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}_k^{[r + 1, i]} \gets \widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}_k^{[r, i]}$, set
$\widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_k^{[r + 1, i]} \gets \widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_k^{[r, i]}$. \;
}
{\color{black} \textbf{break} if $| \ell(\ensuremath{\boldsymbol{\mathbf{y}}} \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}^{[r+1]}) - \ell(\ensuremath{\boldsymbol{\mathbf{y}}} \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}^{[r]}) |$ is sufficiently small.} \;
}
\KwOut{$\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_k=\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{k}^{[r,i]}$, $\widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}_k=\widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}_{k}^{[r,i]}$ and $\widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_k = \widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_{k}^{[r,i]}$ for all $k=1,\ldots,K$}
\end{algorithm}
{\color{black} The iteratively reweighted least squares algorithm \citep[][]{rigby2005generalized} for the estimation of distributional regression models is given in Algorithm \ref{alg:batch_gamlss} 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]}$ \citep[][see p. 113]{stasinopoulos2024generalized}. This has facilitated the development of multiple extensions to the batch GAMLSS, such as \texttt{gamlss.boost} \citep{mayr2012generalized, hofner2014gamboostlss}, mixtures with Neural Networks \citep{rugamer2024semi}, and regularized approaches as in \cite{groll2019lasso} and \texttt{gamlss.lasso} \citep{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{remark:gaussian_distributional_regression} introduces a simple Gaussian example and gives the concrete working vectors and weights.
\begin{remark}\label{remark:gaussian_distributional_regression}
When 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}\label{eq:working_vector_remark}
\ensuremath{\boldsymbol{\mathbf{z}}}_1 &= \ensuremath{\boldsymbol{\mathbf{z}}}_{\mu} = \ensuremath{\boldsymbol{\mathbf{y}}} \quad \text{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 \text{and} &\ensuremath{\boldsymbol{\mathbf{W}}}_{2,2} &= \operatorname{diag}\left((2\ensuremath{\boldsymbol{\mathbf{\sigma}}}^4)^{-1}\right).
\end{align}
\end{remark}
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$ \cite[see e.g.][Chapter~6]{francq2019garch}. For this approach, consistent model selection using the adaptive LASSO and asymptotic normality of the estimators have been shown by \cite{dette2013least} for independent observations and \cite{ziel2016iteratively} for time series settings. The weights~$\ensuremath{\boldsymbol{\mathbf{W}}}_{2,2}$ in Remark \ref{remark:gaussian_distributional_regression} additionally correct for heteroskedasticity in the estimation of the conditional scale parameter. \cite{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{eq:working_vector_remark} in Remark \ref{remark:gaussian_distributional_regression}) 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 \citep{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 \citet{ziel2022m5}.}
\section{Online Estimation of Regularized Linear Models}\label{sec:online_regression}
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{sec:recursive_least_squares}) and the incremental estimation of $L_1$ and $L_2$ penalized linear models using online coordinate descent~(OCD, see Section \ref{sec:online_coordinate_descent}).
\subsection{Weighted Recursive Least Squares}\label{sec:recursive_least_squares}
The following is a condensed introduction \cite[e.g.][Chapter 10-12]{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{sec:method}.} The weighted least squares regression problem estimates the coefficients $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_N^\text{WLS}$ that minimize the loss:
\begin{equation}
\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}^\text{WLS}_N = \mathop{\rm arg\,min}_{\ensuremath{\boldsymbol{\mathbf{\beta}}}} \left\{ \left\lVert \ensuremath{\boldsymbol{\mathbf{W}}}_N^\frac{1}{2} (\ensuremath{\boldsymbol{\mathbf{y}}}_N - \ensuremath{\boldsymbol{\mathbf{X}}}_N\ensuremath{\boldsymbol{\mathbf{\beta}}}) \right\rVert^2_2 \right\} = (\ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}_N^\top}\ensuremath{\boldsymbol{\mathbf{W}}}_N\ensuremath{\boldsymbol{\mathbf{X}}}_N)^{-1} \ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}_N^\top}\ensuremath{\boldsymbol{\mathbf{W}}}_N \ensuremath{\boldsymbol{\mathbf{y}}}_N,
\end{equation}
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 \emph{sample} weights. We will return to these assumptions in Section \ref{sec:method} for the case of \emph{estimation} weights. In either case, the weights should be living on a consistent scale. By employing the Sherman-Morrison formula (see Appendix \ref{app:gram_update}), we can efficiently update the Gramian matrix by taking
\begin{align}
\ensuremath{\boldsymbol{\mathbf{G}}}_{N+1}^{-1} = \ensuremath{\boldsymbol{\mathbf{G}}}_{N}^{-1} - \dfrac{
w_{N+1}\ensuremath{\boldsymbol{\mathbf{G}}}_{N}^{-1}\ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}^\top}\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}\ensuremath{\boldsymbol{\mathbf{G}}}_{N}^{-1}
}{
1+ w_{N+1}\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}\ensuremath{\boldsymbol{\mathbf{G}}}_{N}^{-1}\ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}^\top}
}\label{eq:update_weighted_inverse_gramian}
\end{align}
and we receive
\begin{align}
\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}^{\text{WLS}}_{N+1}
= (\ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}_{N+1}^\top}\ensuremath{\boldsymbol{\mathbf{W}}}_{N+1}\ensuremath{\boldsymbol{\mathbf{X}}}_{N+1})^{-1} \ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}_{N+1}^\top}\ensuremath{\boldsymbol{\mathbf{W}}}_{N+1} \ensuremath{\boldsymbol{\mathbf{y}}}_{N+1}
= \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}^{\text{WLS}}_{N} + \ensuremath{\boldsymbol{\mathbf{G}}}_{N+1}^{-1} w_{N+1}\ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}^\top}\underbrace{\left(y_{N+1} - \ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}^{\text{WLS}}_{N} \right)}_{\text{Forecast error in }N+1} \label{eq:update_equation_wls}
\end{align}
where the latter Equation \ref{eq:update_equation_wls} 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}
\begin{equation}
\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_N = \mathop{\rm arg\,min}_{\ensuremath{\boldsymbol{\mathbf{\beta}}}} \left\{ \left\lVert
\ensuremath{\boldsymbol{\mathbf{\Gamma}}}_N
^\frac{1}{2}
(\ensuremath{\boldsymbol{\mathbf{y}}}_N- \ensuremath{\boldsymbol{\mathbf{X}}}_N\ensuremath{\boldsymbol{\mathbf{ \beta}}}) \right\rVert^2_2 \right\},
\end{equation}
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
\begin{equation}
\ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}^\top}_{N+1}\ensuremath{\boldsymbol{\mathbf{\Gamma}}}_{N+1}\ensuremath{\boldsymbol{\mathbf{X}}}_{N+1} = (1 -\gamma) (\ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}_{N}^\top}\ensuremath{\boldsymbol{\mathbf{\Gamma}}}_{N}\ensuremath{\boldsymbol{\mathbf{X}}}_N)+ \ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}^\top}_{N+1}\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}
\end{equation}
and the {\color{black} corresponding} update equation for the exponentially discounted inverted Gram matrix
\begin{equation}
\begin{aligned}
(\ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}^\top}_{N+1}\ensuremath{\boldsymbol{\mathbf{\Gamma}}}_{N+1}\ensuremath{\boldsymbol{\mathbf{X}}}_{N+1})^{-1}
= \frac{1}{1 -\gamma} \left((\ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}^\top}_{N}\ensuremath{\boldsymbol{\mathbf{\Gamma}}}_{N}\ensuremath{\boldsymbol{\mathbf{X}}}_{N})^{-1} -
\frac{
(\ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}^\top}_{N}\ensuremath{\boldsymbol{\mathbf{\Gamma}}}_{N}\ensuremath{\boldsymbol{\mathbf{X}}}_{N})^{-1} \ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}^\top}_{N+1}\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1} (\ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}^\top}_{N}\ensuremath{\boldsymbol{\mathbf{\Gamma}}}_{N}\ensuremath{\boldsymbol{\mathbf{X}}}_{N})^{-1}
}{
(1 - \gamma) + \ensuremath{\boldsymbol{\mathbf{X}}}_{N+1}(\ensuremath{\ensuremath{\boldsymbol{\mathbf{X}}}^\top}_{N}\ensuremath{\boldsymbol{\mathbf{\Gamma}}}_{N}\ensuremath{\boldsymbol{\mathbf{X}}}_{N})^{-1}\ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}^\top}_{N+1}
}\right)
\end{aligned}\label{eq:update_exponentially_discounted_inverse_gramian}
\end{equation}
which allows for an efficient online update of the OLS coefficients under exponential forgetting by using Equation \ref{eq:update_equation_wls} and setting $w_{N+1} = 1$. {\color{black} Note the subtle difference between Equation \ref{eq:update_weighted_inverse_gramian} for updating the inverse Gramian matrix with sample weights $\ensuremath{\boldsymbol{\mathbf{W}}}_N$ and Equation \ref{eq:update_exponentially_discounted_inverse_gramian} 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
\begin{equation} \label{eq:full_gramian_update}
\ensuremath{\boldsymbol{\mathbf{G}}}_{N+1} = (1 - \gamma) \ensuremath{\boldsymbol{\mathbf{G}}}_{N} + w_{N+1}\ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}^\top}\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}
\end{equation}
and accordingly the update of the inverted Gram matrix $\ensuremath{\boldsymbol{\mathbf{G}}}_{N+1}^{-1}$ can be written as:
\begin{equation} \label{eq:full_inverted_gramian_update}
\ensuremath{\boldsymbol{\mathbf{G}}}^{-1}_{N+1} = \dfrac{1}{1-\gamma}\left( \ensuremath{\boldsymbol{\mathbf{G}}}^{-1}_{N} -
\frac{
w_{N+1} \ensuremath{\boldsymbol{\mathbf{G}}}^{-1}_{N} \ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}^\top}_{N+1}\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1} \ensuremath{\boldsymbol{\mathbf{G}}}^{-1}_{N}
}{
(1 - \gamma) + w_{N+1}\left(\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}\ensuremath{\boldsymbol{\mathbf{G}}}^{-1}_{N}\ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}^\top}_{N+1}\right)
}
\right)
\end{equation}
which can be plugged into Equation \ref{eq:update_equation_wls} to update the coefficients $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}$, as in Algorithm \ref{alg:online_ols}.
\begin{algorithm}[h]
\caption{Recursive Weighted Exponentially Discounted Least-Squares}\label{alg:online_ols}
\DontPrintSemicolon
\KwIn{New observations $\ensuremath{\boldsymbol{\mathbf{X}}}_{N+1}, y_{N+1}, w_{N+1}$ and stored $\ensuremath{\boldsymbol{\mathbf{G}}}^{-1}_N, \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_N$.}
Update $\ensuremath{\boldsymbol{\mathbf{G}}}^{-1}_N \rightarrow \ensuremath{\boldsymbol{\mathbf{G}}}^{-1}_{N+1}$ according to Equation \ref{eq:update_weighted_inverse_gramian} resp. \ref{eq:update_exponentially_discounted_inverse_gramian}. \;
{\color{black} Update $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_N \rightarrow \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1}$ according to Equation \ref{eq:update_equation_wls}.} \;
\KwOut{$\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1}$.}
\end{algorithm}
\subsection{Online Coordinate Descent}\label{sec: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 \citep{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:
\begin{equation}
\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_N
= \mathop{\rm arg\,min}_{\ensuremath{\boldsymbol{\mathbf{\beta}}}} \left\{ \left\| \ensuremath{\boldsymbol{\mathbf{W}}}_N^\frac{1}{2} \left( \ensuremath{\boldsymbol{\mathbf{y}}}_N - \ensuremath{\boldsymbol{\mathbf{X}}}_N \ensuremath{\boldsymbol{\mathbf{\beta}}} \right) \right\|_2^2
+ \alpha \lambda\left\| \ensuremath{\boldsymbol{\mathbf{\beta}}} \right\|_1
+ (1 - \alpha) \lambda \left\| \ensuremath{\boldsymbol{\mathbf{\beta}}} \right\|_2
\right\}.
\end{equation}
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, }
\citet{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}
\begin{equation}
\widehat{\beta}_{N,j} \xleftarrow{} \frac{S\left( \sum_{n=1}^{N} w_n x_{n,j}(y_n - \tilde{y}_n^{(j)}), \alpha \lambda \right)}{\sum_{n=1}^{N} w_n x_{n,j}^2 + ( 1 - \alpha) \lambda }\label{eq:lasso_update_batch}
\end{equation}}
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, \cite{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 \cite{angelosante2010online} and \cite{messner2019online} for LASSO estimation, we rewrite the update Equation \ref{eq:lasso_update_batch} 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{sec:online_regression} and Equation \ref{eq:full_gramian_update}).
{\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
\begin{align*}
G_{N,j,j} = \sum_{n=1}^N w_n x_{n,j}^2 && \text{and} && {h}_{N,j} =\sum_{n}^{N} w_n x_{n,j}y_n .
\end{align*} Thus, writing the nominator inside the soft-thresholding function as
\begin{align*}
\sum_{n=1}^{N} w_n x_{n,j}\left(y_n - \tilde{y}_n^{(j)}\right) &= \sum_{n=1}^{N} w_n x_{n,j}\left(y_n - \widehat{y}_n + x_{nj}\widehat{\beta}_j\right)
= \sum_{n=1}^{N} w_n x_{n,j}\left(y_n - \sum_i^Jx_{ni}\widehat{\beta}_i + x_{nj}\widehat{\beta}_j\right) \\
&= \underbrace{{\sum_{n=1}^{N} w_n x_{n,j}y_n}}_{=h_{N,j} }
- \underbrace{\sum_i^J \sum_n^N w_n x_{n,j} x_{n,i}\beta_i}_{=\ensuremath{\boldsymbol{\mathbf{G}}}_{N,j, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}} \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N}}
+ \underbrace{\sum_{n=1}^{N} w_n x_{n,j}^2\widehat{\beta}_{N,j}}_{=G_{N, j, j}\widehat{\beta}_{N,j}}
\end{align*}
and} we receive the {\color{black} update equation for the elastic net as}
{\color{black}
\begin{equation}
\widehat{\beta}_{N,j} \xleftarrow{} \frac{
S\left(h_{N,j} - \ensuremath{\boldsymbol{\mathbf{G}}}_{N, j, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}} \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N} + G_{N,j, j}\widehat{\beta}_{N,j}, \alpha \lambda \right)
}{
G_{N,j,j}+ ( 1 - \alpha) \lambda
}. \label{eq:lasso_update_online}
\end{equation}}
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 \citep{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{alg:online_lasso} presents the algorithm schematically.
\begin{algorithm}[htb]
\caption{{\color{black} Online Elastic Net Regression using Online Coordinate Descent Estimation}, see \cite{angelosante2010online} and \cite{messner2019online}}\label{alg:online_lasso}
\DontPrintSemicolon
\KwIn{New observations $\ensuremath{\boldsymbol{\mathbf{x}}}_{N+1}, y_{N+1}, w_{N+1}$ and stored $\ensuremath{\boldsymbol{\mathbf{G}}}_N, \ensuremath{\boldsymbol{\mathbf{h}}}_N$.}
Update $\ensuremath{\boldsymbol{\mathbf{G}}}_N \rightarrow \ensuremath{\boldsymbol{\mathbf{G}}}_{N+1}$ and $\ensuremath{\boldsymbol{\mathbf{h}}}_N \rightarrow \ensuremath{\boldsymbol{\mathbf{h}}}_{N+1}$ according to Equation \ref{eq:full_gramian_update}. \;
Update $\lambda_{\text{max}} = \max{}\lvert \ensuremath{\boldsymbol{\mathbf{h}}}_{N+1} \rvert$ and initialize $\ensuremath{\boldsymbol{\mathbf{\lambda}}}$ as exponential grid. \;
\For{$\lambda \in \ensuremath{\boldsymbol{\mathbf{\lambda}}}$}{
Set starting coefficients $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1,\lambda} \gets \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N,\lambda}$ \;
\While{\text{\normalfont{not converged}}}{
\ForAll{$j \in 1, \ldots, J$}{
\text{Update $\widehat{\beta}_{N+1, \lambda, j}$ according to Equation \ref{eq:lasso_update_online}} \;
}
Check convergence for $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1, \lambda}$ and proceed to next $\lambda$ if converged. \;
}
}
\KwOut{
$\widehat{ \ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1, \lambda}$ for all $\lambda \in \ensuremath{\boldsymbol{\mathbf{\lambda}}}$}
\end{algorithm}
{\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 \texttt{glmnet} \texttt{R} package \citep{friedman2010regularization, tay2023elastic}, Equation \ref{eq:lasso_update_online} is then formulated as:
\begin{equation*}
\widehat{\beta}_{N,j} \xleftarrow{} \operatorname{clip}\left(
\frac{
S\left(h_{N,j} - \ensuremath{\boldsymbol{\mathbf{G}}}_{N, j, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}} \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N} + G_{N,j, j}\widehat{\beta}_{N,j}, \alpha \lambda \right)
}{
G_{N,j,j}+ ( 1 - \alpha) \lambda
}, L_j, U_j\right).
\end{equation*}}
\section{Online Distributional Regression} \label{sec:method}
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{alg:online_gamlss} 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.}
\subsection{Update Algorithm}
{\color{black} We start again with the distributional regression model from the introduction given in Equation \ref{eq:y_distributed_f_theta} and \ref{eq:distreg_model}:
\begin{equation*}
y_n \sim \mathcal{F}(\ensuremath{\boldsymbol{\mathbf{\theta}}}_{n, \mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}}) \quad \text{and }\quad
g_k(\ensuremath{\boldsymbol{\mathbf{\theta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}, k}) = \ensuremath{\boldsymbol{\mathbf{\eta}}}_{\mathbin{\vcenter{\hbox{\scalebox{0.5}{$\bullet$}}}}, k} = \ensuremath{\boldsymbol{\mathbf{X}}}_k\ensuremath{\boldsymbol{\mathbf{\beta}}}_k
\end{equation*}}
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
\begin{align}
\ensuremath{\boldsymbol{\mathbf{G}}}_{N+1, k}^{[r,i]} &= \left(1-\gamma\right)\ensuremath{\boldsymbol{\mathbf{G}}}_{N, k} + w_{N+1, k,k}^{[r,i]}\left( \ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}^\top}_{N+1} \ensuremath{\boldsymbol{\mathbf{x}}}_{N+1} \right) \label{eq:update_gram_gamlss} \\
\ensuremath{\boldsymbol{\mathbf{h}}}_{N+1, k}^{[r,i]} &= \left(1-\gamma\right)\ensuremath{\boldsymbol{\mathbf{h}}}_{N, k} + w_{N+1, k,k}^{[r,i]}\left( \ensuremath{\ensuremath{\boldsymbol{\mathbf{x}}}^\top}_{N+1} z_{N+1, k}^{[r,i]} \right) \label{eq:update_h_gamlss}
\end{align}}
and subsequently use an online estimation method from Section \ref{sec:online_regression} to update the coefficients inside the iteration. As in the batch algorithm (see Algorithm \ref{alg:batch_gamlss}), we run the outer and inner cycles until convergence.
\begin{algorithm}[!ht]
\caption{Online update for (regularized) distributional regression models.}\label{alg:online_gamlss}
\small
\DontPrintSemicolon
{\color{black} \KwIn{\text{Distribution $\mathcal{F}$, new data $({y}_{N+1}, \ensuremath{\boldsymbol{\mathbf{x}}}_{N+1})$, stored $\ensuremath{\boldsymbol{\mathbf{G}}}_{N, k}, \ensuremath{\boldsymbol{\mathbf{h}}}_{N, k}$ }}}
Initialize the fitted values $\widehat{\theta}_{N+1, k}^{[0, 0]} = \ensuremath{\boldsymbol{\mathbf{x}}}_{N+1, k}\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N, k}$ for $k = 1, \ldots, p$. \;
Evaluate the linear predictors ${\widehat{\eta}^{[0, 0]}_{N+1, k} = g_k(\widehat{\theta}_{N+1, k}^{[0, 0]})}$ for $k = 0, \ldots, p$. \;
\For{$r = 1, \ldots, \text{\normalfont until convergence}$}{
\For{$k = 1, \ldots, K$}{
\For{$i = 0, 1, \ldots$}{
Evaluate $u_{N+1,k}^{[r,i]}$, $w_{N+1,k,k}^{[r,i]}$ and $z_{N+1, k}^{[r,i]}$ (Eq. \ref{eq:gamlss_score}, \ref{eq:gamlss_weight} and \ref{eq:gamlss_working_vector}). \;
{\color{black} Update $\ensuremath{\boldsymbol{\mathbf{G}}}_{N+1, k}^{[r,i]}$ and $\ensuremath{\boldsymbol{\mathbf{h}}}_{N+1, k}^{[r,i]}$ using $z_{N+1, k}^{[r,i]}$, weights $w_{N+1, k,k}^{[r,i]}$ (Eq. \ref{eq:update_gram_gamlss} and \ref{eq:update_h_gamlss})}. \;
{\color{black} Use OCD to update $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1,k, \lambda} \rightarrow \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1, k, \lambda}^{[r,i+1]}$ based on $\ensuremath{\boldsymbol{\mathbf{G}}}_{N+1, k}^{[r,i]}$ and $\ensuremath{\boldsymbol{\mathbf{h}}}_{N+1, k}^{[r,i]}$.}\;
Select the optimal $\lambda$ using IC and set $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1,k}^{[r,i+1]} = \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{ N+1, k, \lambda^\text{opt}}^{[r,i+1]}$. \;
Calculate the updated $\widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}_{N+1, k}^{[r, i+1]}$ and $\widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_{N+1, k}^{[i,r+1]}$ \;
Evaluate the convergence.\;
}
{\color{black} \textbf{break} if $| \ell(\ensuremath{\boldsymbol{\mathbf{y}}} \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}^{[r, i+1]}) - \ell(\ensuremath{\boldsymbol{\mathbf{y}}} \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}^{[r, i]} )|$ is sufficiently small.} \;
Set
$\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1, k}^{[r + 1, 0]} \gets \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1, k}^{[r, i]}$, set
$\widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}_{N+1, k}^{[r + 1, 0]} \gets \widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}_{N+1, k}^{[r, i]}$ and set
$\widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_{N+1, k}^{[r + 1, 0]} \gets \widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_{N+1, k}^{[r, i]}$.
}
{\color{black} \textbf{break} if $| \ell(\ensuremath{\boldsymbol{\mathbf{y}}} \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}^{[r+1]}) - \ell(\ensuremath{\boldsymbol{\mathbf{y}}} \mid \ensuremath{\boldsymbol{\mathbf{\theta}}}^{[r]} )|$ is sufficiently small.} \;
}
\KwOut{$\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1, k}= \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_{N+1, k}^{[r,i]}$, $\widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}_{N+1, k} = \widehat{\ensuremath{\boldsymbol{\mathbf{\eta}}}}_{N+1, k}^{[r,i]}$ and $\widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_{N+1, k}=\widehat{\ensuremath{\boldsymbol{\mathbf{\theta}}}}_{N+1, k}^{[r,i]}$ for all $k=1,\ldots, K$}
\end{algorithm}
{\color{black} The following paragraphs discuss some points of the algorithm in detail.} To make the discussion clearer, we call an \emph{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 \emph{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{sec:online_model_selection}) and the difference between batch and online estimation (Section~\ref{sec:batch_online_theory}) in detail.}
\begin{itemize}
\item 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.
\item 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 \emph{underestimation} of the conditional heteroskedasticity and potentially higher moments.
\item {\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.}
\item {\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{alg:online_gamlss}, 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).}
\item 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.
\item {\color{black} Algorithm \ref{alg:online_gamlss} 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).}
\end{itemize}
\subsection{Online Model Selection}\label{sec:online_model_selection}
{\color{black} As already hinted at in Section \ref{sec:online_coordinate_descent}, 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
\begin{equation}
\operatorname{GIC}(\widehat{\mathcal{L}}, \nu) = - 2 \log(\widehat{\mathcal{L}}) + \nu_0 P + \nu_1 P \log(\tilde{N}) + \nu_2 P \log(\log(\tilde{N}))
\end{equation}
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 \citep{kim2012consistent, kock2016consistent}. {\color{black} From here, we have two options to proceed: we can use \emph{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{alg:batch_gamlss} and Line 9 in Algorithm \ref{alg:online_gamlss}) or \emph{global} model selection, based on the log-likelihood under the full model.
\begin{enumerate}
\item \emph{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}\label{eq:rss_likelihood}
\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).$$
\item \emph{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).$$
\end{enumerate}
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 \texttt{R} \texttt{gamlss.lasso} package \citep{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 \emph{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.} \cite{groll2019lasso} analyze the impact of different shrinkage parameters on batch distributional regression. Also \cite{marcjasz2023distributional} see the advantages of different regularization for the distribution parameters when using distributional neural networks for time series forecasting.
\subsection{Comparison of Online and Batch Estimation}\label{sec:batch_online_theory}
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 \emph{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{sec:online_regression}). Therefore, we set
\begin{equation}\label{eq:online_weight_matrx}
\ensuremath{\boldsymbol{\mathbf{W}}}_{N+1, k,k}^{[r,i]} =
\begin{pmatrix}
(1-\gamma)\ensuremath{\boldsymbol{\mathbf{W}}}_{N,k,k} & \ensuremath{\boldsymbol{\mathbf{0}}} \\
\ensuremath{\boldsymbol{\mathbf{0}}} & w_{N+1, k,k}^{[r,i]}
\end{pmatrix}
\end{equation}
and can only update $w_{N+1, k,k}^{[r,i]}$. Note that due to the update rule in Equation \ref{eq:full_gramian_update}, 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 \emph{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{eq:gamlss_vector_u} to \ref{eq:gamlss_working_vector}. The online weight matrix is retrieved by iteratively inserting the weight matrix in Equation \ref{eq:online_weight_matrx}.
For illustration, building on the standard linear-Gaussian case for the distributional regression in Remark~\ref{remark:gaussian_distributional_regression}, we can examine the difference explicitly for the mean equation:
{\color{black}
\begin{remark} \label{remark:weights_online}
Assume a linear-Gaussian distributional regression with $\theta_1 = \mu$ and $\theta_2 = \sigma^2$ of Remark \ref{remark:gaussian_distributional_regression}. 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}}
}_\text{\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]}}
}_\text{\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}\label{eq:weight_difference}
\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)}_\text{\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)}_\text{\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]}}}_\text{\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 \underline{not} identical, since the online estimation is path-dependent (see the recursive insertion of $w_{N+1, k,k}$ in Equation \ref{eq:online_weight_matrx}).
\end{remark}
}
From Remark \ref{remark:weights_online} 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{eq:weight_difference} 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{eq:weight_difference} 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
\section{Simulation Study}\label{sec:simulation}
{\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
\begin{equation}
y_{n}
\sim \mathcal{N}\left( \ensuremath{\boldsymbol{\mathbf{X}}}_{n}\ensuremath{\boldsymbol{\mathbf{\beta}}}_1 , \;g^{-1}(\ensuremath{\boldsymbol{\mathbf{X}}}_{n}\ensuremath{\boldsymbol{\mathbf{\beta}}}_2)\right)
\end{equation}
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:
\begin{align}
\ensuremath{\boldsymbol{\mathbf{\beta}}}_k \sim \mathcal{U}(-1, 1).
\end{align}
Unless mentioned otherwise, we set $J=10$, use $N=10000$ observations, use $2000$ observations for the initial fit, and execute 10 simulation runs.
\begin{enumerate}
\item 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{fig:simulation_study_results_init_batch} 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.
\item 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{fig:simulation_study_results_timing}. 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.
\item 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{fig:results_chaning_betas} 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{fig:simulation_study_params_forget} in Appendix~\ref{app:simulation_study_additional_results} presents the parameter estimates for selected forget values for all 20 parameters over time.
\end{enumerate}}
\begin{figure}[htbp]
\Description[Difference between the estimated and true coefficients for the distributional regression model in the first simulation study.]{
Difference between the estimated and true coefficients for the distributional regression model in the first simulation study. We compare estimation with OLS and OCD and see that the error for OCD on the location parameter is somewhat similar, but is lower for OCD on small sample sizes for the scale parameter.
}
\centering
\includegraphics[width=\linewidth]{figures/simulation_study_beta_error.pdf}
\caption{The difference $\| \ensuremath{\boldsymbol{\mathbf{\beta}}}_k - \widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_k \|_1$ relative to the selected initial batch size $n_0$. The main panel shows box-plots of all $M=10$ simulations. The top plot gives results for $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_0$, the bottom plot gives results for $\widehat{\ensuremath{\boldsymbol{\mathbf{\beta}}}}_1$.}
\label{fig:simulation_study_results_init_batch}
\end{figure}
\begin{figure}[htbp]
\Description[Difference in computation times between OLS and OCD based estimation.]{
Difference in computation times between OLS and OCD based estimation. We see that the estimation time increases faster for OLS than for OCD due to the quadratic complexity in the IRLS update. For a very small number of features, OLS is faster.
}
\centering
\includegraphics[width=\linewidth]{figures/simulations/rls_vs_ocd_timings_bar.pdf}
\caption{Average estimation time in Milliseconds, conditional on the number of features, using RLS and OCD.}
\label{fig:simulation_study_results_timing}
\end{figure}
\begin{figure}[htb]
\Description[
Various scoring rules, averaged over all runs for different values of the forget in the simulation study with randomly changing structural breaks.
]{
Various scoring rules, averaged over all runs for different values of the forget in the simulation study with randomly changing structural breaks. We see the lowest scores around a forget of 2^8, which corresponds to an effective sample size of 256 observations. The computation time increases with increasing forget.
}
\centering
\includegraphics[width=\linewidth]{figures/simulation_study_changing_betas_performance.pdf}
\caption{Average performance scores (left axis) and estimation time (right axis) across all 10 runs, for different forget values.}
\label{fig:results_chaning_betas}
\end{figure}
\FloatBarrier
\section{Forecasting Study for Electricity Prices}\label{sec:forecasting}
\subsection{Data, Study Design and Models}
The following section presents an application in electricity price forecasting (EPF). We employ the same setting as \cite{marcjasz2023distributional} and \cite{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{tab:results_epf_batch_vs_online}). Furthermore, we compare our approach to established benchmark models from the literature and conduct several ablation studies (see Appendix \ref{app:forecasting_study_additional_results}, Tables \ref{tab:results_epf_benchmark} and \ref{tab:results_epf_ablation}). The following paragraphs describe the German day-ahead electricity market, data set, and models in detail.
\begin{figure}[htb]
\Description[
Time Series plots of electricity prices, residual load and fundamental fuel prices.
]{
Time Series plots of electricity prices, residual load and fundamental fuel prices. Electricity prices are characterized by their high volatility, while the residual load is somewhat more seasonal effects and slightly skewed. The fuel price time series behave more stochastic and show less seasonal pattern.
}
\centering
\includegraphics[width=\linewidth]{figures/electricity_price/time_series_overview.pdf}
\caption{Overview of the time series data for the electricity price forecasting application. The top shows the EPEX day-ahead spot price, the middle panel shows the residual load, i.e., the system load minus renewable energy production, and the bottom panel shows the fundamental fuel prices.}
\label{fig:epf_time_series_overview}
\end{figure}
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} \citep[further information can be found in, e.g.][]{viehmann2017state, marcjasz2023distributional, ziel2018day, lago2018forecasting, epftoolbox}. {\color{black} Figure \ref{fig:epf_time_series_overview} 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 \cite{marcjasz2023distributional}\footnote{See \url{https://github.com/gmarcjasz/distributionalnn}. \cite{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 \cite{brusaferri2024line} and \cite{marcjasz2023distributional}.
We employ a so-called expert-type model as it is common in electricity price forecasting \citep{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:
\begin{align}\label{eq:model}
g_k(\theta_{k,d,h})
&= \beta_{k, 0, h} + \sum_{s \in \{0, 1, \ldots, H\}} \beta_{k, s, h} P_{d-1, s} + \sum_{l=2}^{L=7} \beta_{k, 24+l, h} P_{d-l, h}\\ \nonumber
&+ \beta_{k,31, h}\operatorname{min}(\ensuremath{\boldsymbol{\mathbf{P}}}_{d-1})
+ \beta_{k,32,h}\max(\ensuremath{\boldsymbol{\mathbf{P}}}_{d-1})
+ \beta_{k,33,h}\mathcal{Q}_{10}(\ensuremath{\boldsymbol{\mathbf{P}}}_{d-1})
+ \beta_{k,34,h}\mathcal{Q}_{90}(\ensuremath{\boldsymbol{\mathbf{P}}}_{d-1}) \\ \nonumber
&+ \beta_{k,35,h}\widehat{\operatorname{ResLoad}}_{d}
+ \sum_{\operatorname{WD} \in \text{\{Mon,Tue,Thu,Fri,Sat,Sun,Hol\}}} \beta_{k, 35+w, h} \operatorname{WD}_{d} \\ \nonumber
&+ b_{k,1,h}(\widehat{\operatorname{ResLoad}}_{d, h})
+ b_{k,2}(\operatorname{EUA}_{d-2})
+ b_{k,3}(\operatorname{Gas}_{d-2})
+ b_{k,4}(\operatorname{Coal}_{d-2})
+ b_{k,5}(\operatorname{Oil}_{d-2})
\end{align}
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 \texttt{Python} \texttt{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 \citep[see][for more information on B-Splines]{eilers1996flexible}.
We benchmark forecasts using established, strictly proper probabilistic scoring rules \citep{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 \citep[IS, also known as Winkler Score, see ][]{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 \texttt{scoringrules} package \citep{zanetta_scoringrules_2024}. We evaluate the statistical significance of the difference in predictive accuracy using the Diebold-Mariano test \citep{diebold2002comparing, diebold2015comparing}. {\color{black} For the coverage, we evaluate the Kupiec-test (see \cite{kupiec1995techniques}) per hour $h$. Detailed definitions of the scoring rules and the statistical tests can be found in Appendix~\ref{app:scoring_rules}}. Lastly, we evaluate the computation time for all forecasting studies.
\subsection{Results and Discussion}\label{sec:results}
The following section presents and discusses the results of the forecasting study. Figure~\ref{fig:illustrative_forecast} gives an exemplary forecast. The time-varying width of the prediction intervals illustrates the importance of modelling the conditional volatility. Figure~\ref{fig:results_epf_nonlinear_effect_scale} 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{app:forecasting_study_additional_results}). 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.
\begin{figure}[!ht]
\Description[
Illustrative forecasts and prediction intervals for the week around Christmas 2020. The prediction intervals are time-varying and adjust to the price level and uncertainty.
]{
Illustrative forecasts and prediction intervals for the week around Christmas 2020. The prediction intervals are time-varying and adjust to the price level and uncertainty. On 24th, 26th prediction intervals are wider. On 28th, prediction intervals are very wide and prices turn out negative.
}
\centering
\includegraphics[width=0.75\textwidth]{figures/illustration_intervals.pdf}
\caption{Illustrative Forecast. Prediction Intervals derived from the online probabilistic model. Here, we assume the power prices to follow Johnson's $S_U$ distribution and model all moments conditionally on Equation \ref{eq:model}. The prediction intervals correspond to the $\{0.005, 0.05, 0.1, \ldots, 0.95, 0.995\}$ quantiles of the predictive distribution. Note that the period with extremely low prices corresponds to Christmas.}
\label{fig:illustrative_forecast}
\end{figure}
\begin{figure}
\Description[
Non-linear effects of the residual load on the location and scale parameter for the $t$-distribution.
]{
Non-linear effects of the residual load on the location and scale parameter for the $t$-distribution. The scale parameter (i.e., the price volatility) increases if the residual load is very high or very low, leading to a $u$-shaped non-linear effect, while the strength of the effect changes through the out-of-sample set.
}
\centering
\includegraphics[width=\linewidth]{figures/electricity_price/residual_load_effects_scale.pdf}
\caption{Non-linear effects of the residual load on the scale (volatility) parameter for the $t$-distribution. The background coloring gives the non-linear effect for the value on the $x$-axis, and the $y$-axis denotes the update steps on the test set. The white line corresponds to the actual value. For the location, the effect on the location parameter is steeper towards the lower and upper end of the variable's domain. The scale parameter (i.e., the price volatility) increases if the residual load isS very high or very low, leading to a $u$-shaped non-linear effect, while the strength of the effect changes through the out-of-sample set.}
\label{fig:results_epf_nonlinear_effect_scale}
\end{figure}
\begin{table}[!ht]
\resizebox{\textwidth}{!}{
\begin{tabular}{llllrrrrrrrrr}
\toprule
Setting & Distribution & Method & K & MAE & RMSE & COV50 & COV80 & IS50 & IS80 & CRPS & LS & Time (min) \\
\midrule
\color{black} Batch & \color{black} Normal & \color{black} OLS & \color{black} 1 & {\cellcolor[HTML]{9DB7CA}} \color[HTML]{000000} \color{black} 4.402 & {\cellcolor[HTML]{74AEDA}} \color[HTML]{F1F1F1} \color{black} 6.314 & {\cellcolor[HTML]{8CB3D1}} \color[HTML]{000000} \color{black} 0.525 & {\cellcolor[HTML]{8CB3D1}} \color[HTML]{000000} \color{black} 0.816 & {\cellcolor[HTML]{9FB7CA}} \color[HTML]{000000} \color{black} 25.109 & {\cellcolor[HTML]{AABAC6}} \color[HTML]{000000} \color{black} 21.663 & {\cellcolor[HTML]{A0B7C9}} \color[HTML]{000000} \color{black} 3.252 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 3.256 & {\cellcolor[HTML]{9EB7CA}} \color[HTML]{000000} \color{black} 5.6 \\
\color{black} Batch & \color{black} Normal & \color{black} OLS & \color{black} 2 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} \itshape 3.948 & {\cellcolor[HTML]{62AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 6.083 & {\cellcolor[HTML]{8CB3D1}} \color[HTML]{000000} \color{black} 0.474 & {\cellcolor[HTML]{C2AEAE}} \color[HTML]{000000} \color{black} 0.762 & {\cellcolor[HTML]{63AAE1}} \color[HTML]{F1F1F1} \color{black} 22.231 & {\cellcolor[HTML]{71ADDB}} \color[HTML]{F1F1F1} \color{black} \itshape 19.297 & {\cellcolor[HTML]{6DACDD}} \color[HTML]{F1F1F1} \color{black} \itshape 2.934 & {\cellcolor[HTML]{87B2D3}} \color[HTML]{000000} \color{black} 3.009 & {\cellcolor[HTML]{ACBAC5}} \color[HTML]{000000} \color{black} 9.0 \\
\color{black} Batch & \color{black} Normal & \color{black} LASSO & \color{black} 1 & {\cellcolor[HTML]{97B5CD}} \color[HTML]{000000} \color{black} 4.345 & {\cellcolor[HTML]{73AEDA}} \color[HTML]{F1F1F1} \color{black} 6.302 & {\cellcolor[HTML]{A5B9C7}} \color[HTML]{000000} \color{black} 0.539 & {\cellcolor[HTML]{9DB7CB}} \color[HTML]{000000} \color{black} 0.821 & {\cellcolor[HTML]{9AB6CC}} \color[HTML]{000000} \color{black} 24.845 & {\cellcolor[HTML]{ACBAC5}} \color[HTML]{000000} \color{black} 21.736 & {\cellcolor[HTML]{9DB7CB}} \color[HTML]{000000} \color{black} 3.230 & {\cellcolor[HTML]{D26263}} \color[HTML]{F1F1F1} \color{black} 3.248 & {\cellcolor[HTML]{AFBBC4}} \color[HTML]{000000} \color{black} 10.3 \\
\color{black} Batch & \color{black} Normal & \color{black} LASSO & \color{black} 2 & {\cellcolor[HTML]{65AAE0}} \color[HTML]{F1F1F1} \color{black} 3.928 & {\cellcolor[HTML]{63AAE0}} \color[HTML]{F1F1F1} \color{black} \itshape 6.101 & {\cellcolor[HTML]{8CB3D1}} \color[HTML]{000000} \color{black} 0.473 & {\cellcolor[HTML]{A5B9C7}} \color[HTML]{000000} \color{black} 0.767 & {\cellcolor[HTML]{64AAE0}} \color[HTML]{F1F1F1} \color{black} 22.281 & {\cellcolor[HTML]{66ABE0}} \color[HTML]{F1F1F1} \color{black} 18.867 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} 2.898 & {\cellcolor[HTML]{78AFD9}} \color[HTML]{F1F1F1} \color{black} 2.985 & {\cellcolor[HTML]{C88E8F}} \color[HTML]{F1F1F1} \color{black} 95.7 \\
\color{black} Batch & \color{black} T & \color{black} OLS & \color{black} 1 & {\cellcolor[HTML]{6DACDD}} \color[HTML]{F1F1F1} \color{black} \itshape 3.991 & {\cellcolor[HTML]{5BA8E4}} \color[HTML]{F1F1F1} \color{black} 5.999 & {\cellcolor[HTML]{CE7576}} \color[HTML]{F1F1F1} \color{black} 0.421 & {\cellcolor[HTML]{CA8586}} \color[HTML]{F1F1F1} \color{black} 0.750 & {\cellcolor[HTML]{A8B9C6}} \color[HTML]{000000} \color{black} \itshape 25.553 & {\cellcolor[HTML]{84B1D4}} \color[HTML]{000000} \color{black} \itshape 20.094 & {\cellcolor[HTML]{74AEDA}} \color[HTML]{F1F1F1} \color{black} \itshape 2.975 & {\cellcolor[HTML]{A6B9C7}} \color[HTML]{000000} \color{black} \itshape 3.057 & {\cellcolor[HTML]{C59D9E}} \color[HTML]{F1F1F1} \color{black} 55.5 \\
\color{black} Batch & \color{black} T & \color{black} OLS & \color{black} 2 & {\cellcolor[HTML]{63AAE0}} \color[HTML]{F1F1F1} \color{black} \itshape 3.913 & {\cellcolor[HTML]{60A9E2}} \color[HTML]{F1F1F1} \color{black} \itshape 6.062 & {\cellcolor[HTML]{C0B6B6}} \color[HTML]{000000} \color{black} 0.453 & {\cellcolor[HTML]{C79596}} \color[HTML]{F1F1F1} \color{black} 0.751 & {\cellcolor[HTML]{6AABDE}} \color[HTML]{F1F1F1} \color{black} 22.527 & {\cellcolor[HTML]{66ABE0}} \color[HTML]{F1F1F1} \color{black} \itshape 18.877 & {\cellcolor[HTML]{65AAE0}} \color[HTML]{F1F1F1} \color{black} \itshape 2.888 & {\cellcolor[HTML]{6CACDD}} \color[HTML]{F1F1F1} \color{black} 2.967 & {\cellcolor[HTML]{C59C9D}} \color[HTML]{F1F1F1} \color{black} 56.7 \\
\color{black} Batch & \color{black} T & \color{black} OLS & \color{black} 3 & {\cellcolor[HTML]{62AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 3.902 & {\cellcolor[HTML]{5FA9E2}} \color[HTML]{F1F1F1} \color{black} 6.058 & {\cellcolor[HTML]{C3A5A6}} \color[HTML]{000000} \color{black} 0.451 & {\cellcolor[HTML]{C79596}} \color[HTML]{F1F1F1} \color{black} 0.749 & {\cellcolor[HTML]{6AACDE}} \color[HTML]{F1F1F1} \color{black} 22.593 & {\cellcolor[HTML]{63AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 18.751 & {\cellcolor[HTML]{62AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 2.871 & {\cellcolor[HTML]{79AFD8}} \color[HTML]{F1F1F1} \color{black} \itshape 2.988 & {\cellcolor[HTML]{C69B9B}} \color[HTML]{F1F1F1} \color{black} 59.9 \\
\color{black} Batch & \color{black} T & \color{black} LASSO & \color{black} 1 & {\cellcolor[HTML]{6CACDD}} \color[HTML]{F1F1F1} \color{black} \itshape 3.983 & {\cellcolor[HTML]{5CA8E3}} \color[HTML]{F1F1F1} \color{black} 6.007 & {\cellcolor[HTML]{CE7576}} \color[HTML]{F1F1F1} \color{black} 0.431 & {\cellcolor[HTML]{C3A5A6}} \color[HTML]{000000} \color{black} 0.757 & {\cellcolor[HTML]{A3B8C8}} \color[HTML]{000000} \color{black} 25.322 & {\cellcolor[HTML]{84B1D4}} \color[HTML]{000000} \color{black} 20.065 & {\cellcolor[HTML]{73AEDA}} \color[HTML]{F1F1F1} \color{black} 2.969 & {\cellcolor[HTML]{A3B8C8}} \color[HTML]{000000} \color{black} 3.051 & {\cellcolor[HTML]{C79394}} \color[HTML]{F1F1F1} \color{black} 78.9 \\
\color{black} Batch & \color{black} T & \color{black} LASSO & \color{black} 2 & {\cellcolor[HTML]{61AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 3.900 & {\cellcolor[HTML]{61AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 6.079 & {\cellcolor[HTML]{C59E9E}} \color[HTML]{F1F1F1} \color{black} 0.448 & {\cellcolor[HTML]{C98D8E}} \color[HTML]{F1F1F1} \color{black} 0.750 & {\cellcolor[HTML]{6EACDD}} \color[HTML]{F1F1F1} \color{black} 22.727 & {\cellcolor[HTML]{5FA9E2}} \color[HTML]{F1F1F1} \color{black} \itshape 18.642 & {\cellcolor[HTML]{62AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 2.869 & {\cellcolor[HTML]{63AAE0}} \color[HTML]{F1F1F1} \color{black} 2.955 & {\cellcolor[HTML]{CB8284}} \color[HTML]{F1F1F1} \color{black} 145.3 \\
\color{black} Batch & \color{black} T & \color{black} LASSO & \color{black} 3 & {\cellcolor[HTML]{61AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 3.897 & {\cellcolor[HTML]{61AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 6.081 & {\cellcolor[HTML]{C3A5A6}} \color[HTML]{000000} \color{black} 0.449 & {\cellcolor[HTML]{C98D8E}} \color[HTML]{F1F1F1} \color{black} 0.749 & {\cellcolor[HTML]{6EACDD}} \color[HTML]{F1F1F1} \color{black} 22.717 & {\cellcolor[HTML]{5EA9E3}} \color[HTML]{F1F1F1} \color{black} \itshape 18.581 & {\cellcolor[HTML]{61AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 2.863 & {\cellcolor[HTML]{63AAE0}} \color[HTML]{F1F1F1} \color{black} 2.954 & {\cellcolor[HTML]{CD7879}} \color[HTML]{F1F1F1} \color{black} 207.8 \\
\color{black} Batch & \color{black} JSU & \color{black} OLS & \color{black} 1 & {\cellcolor[HTML]{72ADDB}} \color[HTML]{F1F1F1} \color{black} \itshape 4.032 & {\cellcolor[HTML]{5FA9E2}} \color[HTML]{F1F1F1} \color{black} 6.047 & {\cellcolor[HTML]{BEBEBE}} \color[HTML]{000000} \color{black} 0.461 & {\cellcolor[HTML]{A5B9C7}} \color[HTML]{000000} \color{black} 0.770 & {\cellcolor[HTML]{9AB6CC}} \color[HTML]{000000} \color{black} \itshape 24.867 & {\cellcolor[HTML]{87B2D3}} \color[HTML]{000000} \color{black} \itshape 20.163 & {\cellcolor[HTML]{78AFD9}} \color[HTML]{F1F1F1} \color{black} \itshape 2.998 & {\cellcolor[HTML]{AFBBC4}} \color[HTML]{000000} \color{black} \itshape 3.072 & {\cellcolor[HTML]{C4A3A4}} \color[HTML]{F1F1F1} \color{black} 44.3 \\
\color{black} Batch & \color{black} JSU & \color{black} OLS & \color{black} 2 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} \itshape 3.946 & {\cellcolor[HTML]{AEBAC4}} \color[HTML]{000000} \color{black} 7.099 & {\cellcolor[HTML]{C79596}} \color[HTML]{F1F1F1} \color{black} 0.435 & {\cellcolor[HTML]{D06D6F}} \color[HTML]{F1F1F1} \color{black} 0.738 & {\cellcolor[HTML]{76AED9}} \color[HTML]{F1F1F1} \color{black} \itshape 23.115 & {\cellcolor[HTML]{61AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 18.708 & {\cellcolor[HTML]{65AAE0}} \color[HTML]{F1F1F1} \color{black} \itshape 2.885 & {\cellcolor[HTML]{6AACDE}} \color[HTML]{F1F1F1} \color{black} \itshape 2.966 & {\cellcolor[HTML]{C4A1A1}} \color[HTML]{F1F1F1} \color{black} 47.9 \\
\color{black} Batch & \color{black} JSU & \color{black} OLS & \color{black} 3 & {\cellcolor[HTML]{7BAFD7}} \color[HTML]{F1F1F1} \color{black} \itshape 4.108 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 14.976 & {\cellcolor[HTML]{D16566}} \color[HTML]{F1F1F1} \color{black} 0.419 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 0.704 & {\cellcolor[HTML]{9FB7CA}} \color[HTML]{000000} \color{black} \itshape 25.094 & {\cellcolor[HTML]{80B0D5}} \color[HTML]{F1F1F1} \color{black} \itshape 19.915 & {\cellcolor[HTML]{7CAFD7}} \color[HTML]{F1F1F1} \color{black} \itshape 3.024 & {\cellcolor[HTML]{70ADDC}} \color[HTML]{F1F1F1} \color{black} \itshape 2.973 & {\cellcolor[HTML]{C69898}} \color[HTML]{F1F1F1} \color{black} 67.0 \\
\color{black} Batch & \color{black} JSU & \color{black} OLS & \color{black} 4 & {\cellcolor[HTML]{C4A1A1}} \color[HTML]{F1F1F1} \color{black} 4.999 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 8584.146 & {\cellcolor[HTML]{D16566}} \color[HTML]{F1F1F1} \color{black} 0.419 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 0.705 & {\cellcolor[HTML]{CE7475}} \color[HTML]{F1F1F1} \color{black} 31.168 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 33.018 & {\cellcolor[HTML]{D35F60}} \color[HTML]{F1F1F1} \color{black} 4.217 & {\cellcolor[HTML]{74AEDA}} \color[HTML]{F1F1F1} \color{black} \itshape 2.979 & {\cellcolor[HTML]{C69898}} \color[HTML]{F1F1F1} \color{black} 67.1 \\
\color{black} Batch & \color{black} JSU & \color{black} LASSO & \color{black} 1 & {\cellcolor[HTML]{6BACDD}} \color[HTML]{F1F1F1} \color{black} \itshape 3.978 & {\cellcolor[HTML]{5CA9E3}} \color[HTML]{F1F1F1} \color{black} 6.021 & {\cellcolor[HTML]{C2AEAE}} \color[HTML]{000000} \color{black} 0.469 & {\cellcolor[HTML]{8CB3D1}} \color[HTML]{000000} \color{black} 0.780 & {\cellcolor[HTML]{91B4CF}} \color[HTML]{000000} \color{black} \itshape 24.418 & {\cellcolor[HTML]{80B0D6}} \color[HTML]{F1F1F1} \color{black} \itshape 19.890 & {\cellcolor[HTML]{71ADDB}} \color[HTML]{F1F1F1} \color{black} \itshape 2.958 & {\cellcolor[HTML]{A2B8C9}} \color[HTML]{000000} \color{black} \itshape 3.050 & {\cellcolor[HTML]{CA8889}} \color[HTML]{F1F1F1} \color{black} 116.1 \\
\color{black} Batch & \color{black} JSU & \color{black} LASSO & \color{black} 2 & {\cellcolor[HTML]{60A9E2}} \color[HTML]{F1F1F1} \color{black} \itshape 3.888 & {\cellcolor[HTML]{7AAFD8}} \color[HTML]{F1F1F1} \color{black} 6.397 & {\cellcolor[HTML]{CE7576}} \color[HTML]{F1F1F1} \color{black} 0.432 & {\cellcolor[HTML]{D06D6F}} \color[HTML]{F1F1F1} \color{black} 0.741 & {\cellcolor[HTML]{7AAFD8}} \color[HTML]{F1F1F1} \color{black} \itshape 23.304 & {\cellcolor[HTML]{5EA9E3}} \color[HTML]{F1F1F1} \color{black} \itshape 18.589 & {\cellcolor[HTML]{5FA9E2}} \color[HTML]{F1F1F1} \color{black} \itshape 2.848 & {\cellcolor[HTML]{61AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 2.951 & {\cellcolor[HTML]{CE7577}} \color[HTML]{F1F1F1} \color{black} 224.8 \\
\color{black} Batch & \color{black} JSU & \color{black} LASSO & \color{black} 3 & {\cellcolor[HTML]{5BA8E4}} \color[HTML]{F1F1F1} \color{black} \itshape 3.851 & {\cellcolor[HTML]{A6B9C7}} \color[HTML]{000000} \color{black} 6.982 & {\cellcolor[HTML]{D06D6F}} \color[HTML]{F1F1F1} \color{black} 0.431 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 0.728 & {\cellcolor[HTML]{75AEDA}} \color[HTML]{F1F1F1} \color{black} \itshape 23.038 & {\cellcolor[HTML]{5AA8E4}} \color[HTML]{F1F1F1} \bfseries \color{black} \itshape 18.420 & {\cellcolor[HTML]{5AA8E4}} \color[HTML]{F1F1F1} \color{black} \itshape 2.822 & {\cellcolor[HTML]{5BA8E4}} \color[HTML]{F1F1F1} \color{black} \itshape 2.942 & {\cellcolor[HTML]{D16668}} \color[HTML]{F1F1F1} \color{black} 389.2 \\
\color{black} Batch & \color{black} JSU & \color{black} LASSO & \color{black} 4 & {\cellcolor[HTML]{5AA8E4}} \color[HTML]{F1F1F1} \bfseries \color{black} \itshape 3.839 & {\cellcolor[HTML]{A5B8C8}} \color[HTML]{000000} \color{black} 6.968 & {\cellcolor[HTML]{CC7D7E}} \color[HTML]{F1F1F1} \color{black} 0.431 & {\cellcolor[HTML]{D16566}} \color[HTML]{F1F1F1} \color{black} 0.733 & {\cellcolor[HTML]{72ADDB}} \color[HTML]{F1F1F1} \color{black} \itshape 22.950 & {\cellcolor[HTML]{5CA9E3}} \color[HTML]{F1F1F1} \color{black} \itshape 18.513 & {\cellcolor[HTML]{5AA8E4}} \color[HTML]{F1F1F1} \bfseries \color{black} \itshape 2.821 & {\cellcolor[HTML]{5AA8E4}} \color[HTML]{F1F1F1} \bfseries \color{black} \itshape 2.940 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 551.6 \\
\midrule
\color{black} Online & \color{black} Normal & \color{black} OLS & \color{black} 1 & {\cellcolor[HTML]{98B6CC}} \color[HTML]{000000} \color{black} \itshape 4.355 & {\cellcolor[HTML]{72ADDB}} \color[HTML]{F1F1F1} \color{black} \itshape 6.285 & {\cellcolor[HTML]{ADBAC4}} \color[HTML]{000000} \color{black} 0.532 & {\cellcolor[HTML]{94B5CE}} \color[HTML]{000000} \color{black} 0.815 & {\cellcolor[HTML]{9BB6CB}} \color[HTML]{000000} \color{black} \itshape 24.895 & {\cellcolor[HTML]{A7B9C7}} \color[HTML]{000000} \color{black} \itshape 21.530 & {\cellcolor[HTML]{9BB6CB}} \color[HTML]{000000} \color{black} \itshape 3.222 & {\cellcolor[HTML]{D35F60}} \color[HTML]{F1F1F1} \color{black} 3.254 & {\cellcolor[HTML]{5BA8E4}} \color[HTML]{F1F1F1} \color{black} \itshape 0.5 \\
\color{black} Online & \color{black} Normal & \color{black} OLS & \color{black} 2 & {\cellcolor[HTML]{6DACDD}} \color[HTML]{F1F1F1} \color{black} 3.990 & {\cellcolor[HTML]{65AAE0}} \color[HTML]{F1F1F1} \color{black} 6.125 & {\cellcolor[HTML]{73AEDA}} \color[HTML]{F1F1F1} \bfseries \color{black} 0.498 & {\cellcolor[HTML]{94B5CE}} \color[HTML]{000000} \color{black} 0.785 & {\cellcolor[HTML]{5CA9E3}} \color[HTML]{F1F1F1} \color{black} \itshape 21.937 & {\cellcolor[HTML]{79AFD8}} \color[HTML]{F1F1F1} \color{black} 19.618 & {\cellcolor[HTML]{73AEDA}} \color[HTML]{F1F1F1} \color{black} 2.972 & {\cellcolor[HTML]{87B2D3}} \color[HTML]{000000} \color{black} 3.008 & {\cellcolor[HTML]{5AA8E4}} \color[HTML]{F1F1F1} \color{black} \itshape 0.5 \\
\color{black} Online & \color{black} Normal & \color{black} LASSO & \color{black} 1 & {\cellcolor[HTML]{98B6CC}} \color[HTML]{000000} \color{black} 4.355 & {\cellcolor[HTML]{73AEDA}} \color[HTML]{F1F1F1} \color{black} 6.305 & {\cellcolor[HTML]{A5B9C7}} \color[HTML]{000000} \color{black} 0.536 & {\cellcolor[HTML]{9DB7CB}} \color[HTML]{000000} \color{black} 0.819 & {\cellcolor[HTML]{9CB6CB}} \color[HTML]{000000} \color{black} 24.947 & {\cellcolor[HTML]{ACBAC5}} \color[HTML]{000000} \color{black} 21.730 & {\cellcolor[HTML]{9DB7CA}} \color[HTML]{000000} \color{black} 3.234 & {\cellcolor[HTML]{D26062}} \color[HTML]{F1F1F1} \color{black} 3.251 & {\cellcolor[HTML]{64AAE0}} \color[HTML]{F1F1F1} \color{black} \itshape 0.8 \\
\color{black} Online & \color{black} Normal & \color{black} LASSO & \color{black} 2 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} 3.941 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} 6.145 & {\cellcolor[HTML]{6AACDE}} \color[HTML]{F1F1F1} \color{black} 0.497 & {\cellcolor[HTML]{7BAFD7}} \color[HTML]{F1F1F1} \color{black} 0.788 & {\cellcolor[HTML]{5AA8E4}} \color[HTML]{F1F1F1} \bfseries \color{black} \itshape 21.822 & {\cellcolor[HTML]{65AAE0}} \color[HTML]{F1F1F1} \color{black} 18.842 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} 2.898 & {\cellcolor[HTML]{75AEDA}} \color[HTML]{F1F1F1} \color{black} 2.981 & {\cellcolor[HTML]{75AEDA}} \color[HTML]{F1F1F1} \color{black} \itshape 1.4 \\
\color{black} Online & \color{black} T & \color{black} OLS & \color{black} 1 & {\cellcolor[HTML]{74AEDA}} \color[HTML]{F1F1F1} \color{black} 4.048 & {\cellcolor[HTML]{5CA8E3}} \color[HTML]{F1F1F1} \color{black} 6.008 & {\cellcolor[HTML]{CE7576}} \color[HTML]{F1F1F1} \color{black} 0.413 & {\cellcolor[HTML]{CE7576}} \color[HTML]{F1F1F1} \color{black} 0.743 & {\cellcolor[HTML]{B2BBC3}} \color[HTML]{000000} \color{black} 26.077 & {\cellcolor[HTML]{8CB3D1}} \color[HTML]{000000} \color{black} 20.401 & {\cellcolor[HTML]{7BAFD7}} \color[HTML]{F1F1F1} \color{black} 3.018 & {\cellcolor[HTML]{B0BBC3}} \color[HTML]{000000} \color{black} 3.073 & {\cellcolor[HTML]{66ABE0}} \color[HTML]{F1F1F1} \color{black} \itshape 0.8 \\
\color{black} Online & \color{black} T & \color{black} OLS & \color{black} 2 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} 3.950 & {\cellcolor[HTML]{62AAE1}} \color[HTML]{F1F1F1} \color{black} 6.081 & {\cellcolor[HTML]{8CB3D1}} \color[HTML]{000000} \color{black} 0.477 & {\cellcolor[HTML]{ADBAC4}} \color[HTML]{000000} \color{black} 0.774 & {\cellcolor[HTML]{62AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 22.193 & {\cellcolor[HTML]{6EACDD}} \color[HTML]{F1F1F1} \color{black} 19.185 & {\cellcolor[HTML]{6BACDD}} \color[HTML]{F1F1F1} \color{black} 2.925 & {\cellcolor[HTML]{6EACDC}} \color[HTML]{F1F1F1} \color{black} 2.971 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} \itshape 0.8 \\
\color{black} Online & \color{black} T & \color{black} OLS & \color{black} 3 & {\cellcolor[HTML]{69ABDE}} \color[HTML]{F1F1F1} \color{black} 3.959 & {\cellcolor[HTML]{62AAE1}} \color[HTML]{F1F1F1} \color{black} 6.081 & {\cellcolor[HTML]{73AEDA}} \color[HTML]{F1F1F1} \color{black} 0.486 & {\cellcolor[HTML]{8CB3D1}} \color[HTML]{000000} \bfseries \color{black} 0.791 & {\cellcolor[HTML]{61AAE1}} \color[HTML]{F1F1F1} \color{black} \itshape 22.141 & {\cellcolor[HTML]{6EACDC}} \color[HTML]{F1F1F1} \color{black} 19.206 & {\cellcolor[HTML]{6BACDD}} \color[HTML]{F1F1F1} \color{black} 2.922 & {\cellcolor[HTML]{84B1D4}} \color[HTML]{000000} \color{black} 3.003 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} \itshape 0.8 \\
\color{black} Online & \color{black} T & \color{black} LASSO & \color{black} 1 & {\cellcolor[HTML]{6FADDC}} \color[HTML]{F1F1F1} \color{black} 4.011 & {\cellcolor[HTML]{5AA8E4}} \color[HTML]{F1F1F1} \bfseries \color{black} 5.985 & {\cellcolor[HTML]{C98D8E}} \color[HTML]{F1F1F1} \color{black} 0.431 & {\cellcolor[HTML]{C3A5A6}} \color[HTML]{000000} \color{black} 0.757 & {\cellcolor[HTML]{A6B9C7}} \color[HTML]{000000} \color{black} 25.447 & {\cellcolor[HTML]{86B2D3}} \color[HTML]{000000} \color{black} 20.154 & {\cellcolor[HTML]{75AEDA}} \color[HTML]{F1F1F1} \color{black} 2.986 & {\cellcolor[HTML]{A6B9C7}} \color[HTML]{000000} \color{black} 3.057 & {\cellcolor[HTML]{6EACDC}} \color[HTML]{F1F1F1} \color{black} \itshape 1.1 \\
\color{black} Online & \color{black} T & \color{black} LASSO & \color{black} 2 & {\cellcolor[HTML]{64AAE0}} \color[HTML]{F1F1F1} \color{black} 3.922 & {\cellcolor[HTML]{65AAE0}} \color[HTML]{F1F1F1} \color{black} 6.127 & {\cellcolor[HTML]{9DB7CB}} \color[HTML]{000000} \color{black} 0.472 & {\cellcolor[HTML]{94B5CE}} \color[HTML]{000000} \color{black} 0.776 & {\cellcolor[HTML]{64AAE0}} \color[HTML]{F1F1F1} \color{black} \itshape 22.310 & {\cellcolor[HTML]{66ABE0}} \color[HTML]{F1F1F1} \color{black} 18.891 & {\cellcolor[HTML]{66ABE0}} \color[HTML]{F1F1F1} \color{black} 2.892 & {\cellcolor[HTML]{63AAE1}} \color[HTML]{F1F1F1} \color{black} 2.954 & {\cellcolor[HTML]{7CAFD7}} \color[HTML]{F1F1F1} \color{black} \itshape 1.7 \\
\color{black} Online & \color{black} T & \color{black} LASSO & \color{black} 3 & {\cellcolor[HTML]{63AAE0}} \color[HTML]{F1F1F1} \color{black} 3.918 & {\cellcolor[HTML]{63AAE0}} \color[HTML]{F1F1F1} \color{black} 6.109 & {\cellcolor[HTML]{8CB3D1}} \color[HTML]{000000} \color{black} 0.473 & {\cellcolor[HTML]{94B5CE}} \color[HTML]{000000} \color{black} 0.776 & {\cellcolor[HTML]{65AAE0}} \color[HTML]{F1F1F1} \color{black} \itshape 22.339 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} 18.923 & {\cellcolor[HTML]{66ABE0}} \color[HTML]{F1F1F1} \color{black} 2.893 & {\cellcolor[HTML]{62AAE1}} \color[HTML]{F1F1F1} \color{black} 2.953 & {\cellcolor[HTML]{8DB3D1}} \color[HTML]{000000} \color{black} \itshape 3.1 \\
\color{black} Online & \color{black} JSU & \color{black} OLS & \color{black} 1 & {\cellcolor[HTML]{79AFD8}} \color[HTML]{F1F1F1} \color{black} 4.095 & {\cellcolor[HTML]{5FA9E2}} \color[HTML]{F1F1F1} \color{black} 6.053 & {\cellcolor[HTML]{C0B6B6}} \color[HTML]{000000} \color{black} 0.454 & {\cellcolor[HTML]{BEBEBE}} \color[HTML]{000000} \color{black} 0.763 & {\cellcolor[HTML]{A5B8C8}} \color[HTML]{000000} \color{black} 25.396 & {\cellcolor[HTML]{8FB4D0}} \color[HTML]{000000} \color{black} 20.503 & {\cellcolor[HTML]{80B0D6}} \color[HTML]{F1F1F1} \color{black} 3.045 & {\cellcolor[HTML]{BEBDBD}} \color[HTML]{000000} \color{black} 3.096 & {\cellcolor[HTML]{6BACDD}} \color[HTML]{F1F1F1} \color{black} \itshape 1.0 \\
\color{black} Online & \color{black} JSU & \color{black} OLS & \color{black} 2 & {\cellcolor[HTML]{71ADDB}} \color[HTML]{F1F1F1} \color{black} 4.022 & {\cellcolor[HTML]{90B4CF}} \color[HTML]{000000} \color{black} 6.683 & {\cellcolor[HTML]{D16566}} \color[HTML]{F1F1F1} \color{black} 0.428 & {\cellcolor[HTML]{CE7576}} \color[HTML]{F1F1F1} \color{black} 0.734 & {\cellcolor[HTML]{84B1D4}} \color[HTML]{000000} \color{black} 23.736 & {\cellcolor[HTML]{6EACDD}} \color[HTML]{F1F1F1} \color{black} 19.180 & {\cellcolor[HTML]{6FADDC}} \color[HTML]{F1F1F1} \color{black} 2.947 & {\cellcolor[HTML]{8BB3D2}} \color[HTML]{000000} \color{black} 3.014 & {\cellcolor[HTML]{6CACDD}} \color[HTML]{F1F1F1} \color{black} \itshape 1.0 \\
\color{black} Online & \color{black} JSU & \color{black} OLS & \color{black} 3 & {\cellcolor[HTML]{C59E9E}} \color[HTML]{F1F1F1} \color{black} 5.033 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 14.133 & {\cellcolor[HTML]{D16566}} \color[HTML]{F1F1F1} \color{black} 0.401 & {\cellcolor[HTML]{D06D6F}} \color[HTML]{F1F1F1} \color{black} 0.681 & {\cellcolor[HTML]{D06D6F}} \color[HTML]{F1F1F1} \color{black} 31.619 & {\cellcolor[HTML]{C98D8E}} \color[HTML]{F1F1F1} \color{black} 24.982 & {\cellcolor[HTML]{C69899}} \color[HTML]{F1F1F1} \color{black} 3.735 & {\cellcolor[HTML]{C4A1A1}} \color[HTML]{F1F1F1} \color{black} 3.142 & {\cellcolor[HTML]{6DACDD}} \color[HTML]{F1F1F1} \color{black} \itshape 1.0 \\
\color{black} Online & \color{black} JSU & \color{black} OLS & \color{black} 4 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} $>10^4$ & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} $>10^4$ & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 0.396 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 0.674 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} $>10^4$ & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} $>10^4$ & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} $>10^4$ & {\cellcolor[HTML]{CA8889}} \color[HTML]{F1F1F1} \color{black} 3.184 & {\cellcolor[HTML]{6EACDD}} \color[HTML]{F1F1F1} \color{black} \itshape 1.0 \\
\color{black} Online & \color{black} JSU & \color{black} LASSO & \color{black} 1 & {\cellcolor[HTML]{73AEDA}} \color[HTML]{F1F1F1} \color{black} 4.041 & {\cellcolor[HTML]{5CA8E3}} \color[HTML]{F1F1F1} \color{black} 6.010 & {\cellcolor[HTML]{C2AEAE}} \color[HTML]{000000} \color{black} 0.458 & {\cellcolor[HTML]{C0B6B6}} \color[HTML]{000000} \color{black} 0.768 & {\cellcolor[HTML]{9DB7CB}} \color[HTML]{000000} \color{black} 24.987 & {\cellcolor[HTML]{87B2D3}} \color[HTML]{000000} \color{black} 20.208 & {\cellcolor[HTML]{79AFD8}} \color[HTML]{F1F1F1} \color{black} 3.003 & {\cellcolor[HTML]{ADBAC4}} \color[HTML]{000000} \color{black} 3.068 & {\cellcolor[HTML]{75AEDA}} \color[HTML]{F1F1F1} \color{black} \itshape 1.4 \\
\color{black} Online & \color{black} JSU & \color{black} LASSO & \color{black} 2 & {\cellcolor[HTML]{68ABDF}} \color[HTML]{F1F1F1} \color{black} 3.954 & {\cellcolor[HTML]{72ADDB}} \color[HTML]{F1F1F1} \color{black} 6.295 & {\cellcolor[HTML]{D06D6F}} \color[HTML]{F1F1F1} \color{black} 0.430 & {\cellcolor[HTML]{CC7D7E}} \color[HTML]{F1F1F1} \color{black} 0.738 & {\cellcolor[HTML]{83B1D5}} \color[HTML]{000000} \color{black} 23.710 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} 18.920 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} 2.900 & {\cellcolor[HTML]{79AFD8}} \color[HTML]{F1F1F1} \color{black} 2.986 & {\cellcolor[HTML]{82B1D5}} \color[HTML]{000000} \color{black} \itshape 2.1 \\
\color{black} Online & \color{black} JSU & \color{black} LASSO & \color{black} 3 & {\cellcolor[HTML]{66ABE0}} \color[HTML]{F1F1F1} \color{black} 3.932 & {\cellcolor[HTML]{7DB0D7}} \color[HTML]{F1F1F1} \color{black} 6.422 & {\cellcolor[HTML]{D06D6F}} \color[HTML]{F1F1F1} \color{black} 0.421 & {\cellcolor[HTML]{D06D6F}} \color[HTML]{F1F1F1} \color{black} 0.718 & {\cellcolor[HTML]{86B2D3}} \color[HTML]{000000} \color{black} 23.848 & {\cellcolor[HTML]{67ABDF}} \color[HTML]{F1F1F1} \color{black} 18.938 & {\cellcolor[HTML]{66ABE0}} \color[HTML]{F1F1F1} \color{black} 2.891 & {\cellcolor[HTML]{78AFD9}} \color[HTML]{F1F1F1} \color{black} 2.986 & {\cellcolor[HTML]{91B4CF}} \color[HTML]{000000} \color{black} \itshape 3.5 \\
\color{black} Online & \color{black} JSU & \color{black} LASSO & \color{black} 4 & {\cellcolor[HTML]{5FA9E2}} \color[HTML]{F1F1F1} \color{black} 3.877 & {\cellcolor[HTML]{D35D5F}} \color[HTML]{F1F1F1} \color{black} 8.985 & {\cellcolor[HTML]{D06D6F}} \color[HTML]{F1F1F1} \color{black} 0.419 & {\cellcolor[HTML]{D16566}} \color[HTML]{F1F1F1} \color{black} 0.716 & {\cellcolor[HTML]{82B1D5}} \color[HTML]{000000} \color{black} 23.673 & {\cellcolor[HTML]{63AAE1}} \color[HTML]{F1F1F1} \color{black} 18.764 & {\cellcolor[HTML]{5FA9E2}} \color[HTML]{F1F1F1} \color{black} 2.855 & {\cellcolor[HTML]{6EACDD}} \color[HTML]{F1F1F1} \color{black} 2.971 & {\cellcolor[HTML]{9DB7CB}} \color[HTML]{000000} \color{black} \itshape 5.4 \\
\bottomrule
\end{tabular}}
\centering\includegraphics[width=0.75\textwidth]{tables/legend.pdf}
\caption{
Results for the main experiment: Repeated batch estimation vs online estimation. Distribution refers to the assumed parametric form of the response distribution. Setting refers to whether models are estimated incrementally or retrained on the increasing window training set. The best value in each column is marked \textbf{bold}. We conduct the pairwise DM-test (Online, vs. Batch) for the RMSE, MAE, IS, CRPS, and LS and mark results that cross the $p<0.05$ significance threshold \textit{italic}. For the coverage of the 50\% and 80\% prediction intervals, the background color corresponds to the number of hours where the \cite{kupiec1995techniques}-test is passed at the 5\%-level. Note that the timing corresponds to a full forecasting study, i.e., the estimation of $534 \times 24 = 12.816$ models.
}
\label{tab:results_epf_batch_vs_online}
\end{table}
Table \ref{tab:results_epf_batch_vs_online} 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:
\begin{itemize}
\item Comparing the distributional regression models in Table \ref{tab:results_epf_batch_vs_online} 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.
\item 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{tab:results_epf_ablation}, Panel~A). The LASSO-estimated JSU model is the strongest in our study.
\end{itemize}
Furthermore, we make the following observations when comparing to established benchmarks and conducting ablation studies:
\begin{itemize}
\item 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.
\item Compared with state-of-the art results of \cite{brusaferri2024line}, our approach delivers slightly better performance in terms of the CRPS, but slightly worse coverage. The DDNN approach of \cite{marcjasz2023distributional} yields better results than the distributional regression approaches and the conformalized DDNNs of \cite{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.
\item 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{sec:method}.
\item Table \ref{tab:results_epf_ablation}, 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.
\item Further discussion of the ablation studies can be found in Appendix~\ref{app:forecasting_study_additional_results}. 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.
\end{itemize}
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.
\section{Implementation in the Python Package \texttt{ondil}}\label{sec:python}
We provide an open-source implementation of our algorithm in the \texttt{Python} package \texttt{ondil}. To the best of the authors' knowledge, this package is the first native \texttt{Python} implementation of IRLS-based distributional regression and can, therefore, provide a basis for future extensions.\footnote{The package \texttt{pyNM} implements GAMLSS for \texttt{Python} as binding to the \texttt{R} library \url{https://github.com/ppsp-team/PyNM}.} We provide the following features:
\begin{itemize}
\item Our package is written in an online-first fashion and provides an API compatible with \texttt{scikit-learn} \cite{pedregosa2011scikit}, the major Python machine learning package. Beyond the \texttt{estimator.fit(X,~y)} and \texttt{estimator.predict(X,~y)} methods, our estimator classes provide the \texttt{estimator.update(X,~y)} method to allow for incremental updating.
\item The implementation relies on few dependencies \citep[only \texttt{numpy}, \texttt{numba}, and \texttt{scipy}, see][]{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 \citep[see e.g.][]{shi2016primer, wright2015coordinate}.
\item We implement standard distributions and facilitate extensions based on \texttt{scipy}-distributions.
\end{itemize}
The code for our implementation is open source {\color{black} on \texttt{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 \texttt{ondil} version 0.3.2.}
\section{Discussion and Conclusion}\label{sec: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 \citep[for batch settings see, e.g.][]{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 \texttt{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.}
\bibliographystyle{ACM-Reference-Format}
\bibliography{references}
\begin{acks}
Simon 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.
\end{acks}
\newpage