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.
68,161 characters
Neural networks for nonlinear regression with serially correlated disturbances: Evidence from cloud cover
\title{Neural networks for nonlinear regression with serially correlated disturbances: Evidence from cloud cover}
\date{}
\author[a]{Sebastian Jensen\thanks{
Corresponding author:
Department of Economics and Business Economics and CoRE,
Aarhus University,
Universitetsbyen 51,
8000 Aarhus C, Denmark.
E-mail:\
[email removed]}}}
\author[b]{Siem Jan Koopman}
\affil[a]{Aarhus University, Denmark}
\affil[b]{Vrije Universiteit Amsterdam, The Netherlands}
{\let\newpage\relax\maketitle}
\begin{abstract}
We propose a new treatment of nonlinear regression with serially correlated disturbances that incorporates autoregressive moving average structures into feedforward neural networks. The resulting model provides an alternative to modeling temporal dependence using lagged variables. In simulations, the proposed method accurately recovers regression functions of varying complexity and the underlying error dynamics across a range of time-series lengths and signal-to-noise ratios. Finite-sample properties and out-of-sample predictive performances are shown to be robust to model misspecification induced by omitted lagged variables and incorrect specification of the error dynamics. Cloud cover is an important factor in climate projections. In an empirical study of cloud cover prediction for a grid of locations within and around the Mediterranean Sea, our proposed model yields more accurate predictions than existing methods, including long short-term memory networks. Improvements are observed broadly and are particularly pronounced in mountain areas relative to linear models with serially correlated errors, consistent with the presence of stronger nonlinear effects in cloud composure in such regions.
\end{abstract}
\noindent \emph{Keywords:} Neural networks; Time series; Nonlinear regressions; Serially correlated disturbances; Autoregressive moving average model; Durbin-Levinson; Cloud cover prediction. \\
\newpage
\section{Introduction}
Machine learning methods, including neural networks, are increasingly used in data science, econometrics and statistics \citep{hastie2009elements,Athey2019,Efron2020}. These methods offer considerable flexibility for modeling nonlinear relationships in regression problems, even in high-dimensional settings. Temporal dependence is often accommodated by including a set of lagged dependent (output) and lagged regressor (input) variables in the regression component. The disturbances are typically left unmodeled or treated as coming from a white noise process, under the implicit assumption that all predictable temporal dependence is absorbed by a sufficiently rich set of lagged variables. Any remaining serial correlation in the disturbance process therefore reflects misspecification. This approach of expanding the predictor space may, however, unnecessarily increase the number of predictors, particularly in a setting that is already high-dimensional. Furthermore, in a machine learning context, it is difficult to use statistical inference to guide the choice of variables and their lag lengths. An alternative approach to modeling temporal dependence via lagged variables in the regression component is to capture it through a dynamic stochastic process for the disturbances. Although presented here as alternatives, the two approaches can be combined. A classical time series model, such as the autoregressive moving average (ARMA) model, provides a parsimonious and interpretable representation of temporal dependence. In the linear framework, the use of such structured stochastic models for the disturbance process is well established in the literature. However, integrating structured stochastic models for the disturbances into flexible machine learning frameworks remains largely unexplored and is the focus of this paper.
In this paper, we develop a methodology for flexible nonlinear regression and prediction with serially correlated disturbances by integrating ARMA error structures into feedforward neural networks. The resulting model is denoted as NNARMA. We focus on ARMA dynamics for the disturbances as they provide a parsimonious representation of a wide variety of serial correlation patterns, see \cite{boxjenkins2015} for a textbook treatment. We focus on feedforward neural networks for a number of reasons. They are universal approximators \citep{HORNIK1989, Cybenko1989, LESHNO1993} and are known to perform well in practice on a number of regression and predictions tasks \citep{Kelly2020,Jensen2023}. Theoretical justification is starting to be developed \citep{Farrell2021}. \cite{Schmidt-Hieber2020} shows that feedforward neural networks are superior to classical nonparametric methods, such as kernel-based estimators and splines, if the unknown target function is a composition of simpler functions. To reduce overfitting, neural networks incorporate various regularization techniques \citep{Goodfellow-et-al-2016}. \cite{Bach2017} and \cite{Bauer2019} show that feedforward neural networks solve the curse of dimensionality that affects classical nonparametric methods \citep{Stone1982}.
The proposed estimation procedure is a conceptually straightforward generalization of well-known techniques for linear regression with serially correlated disturbances. We compute a Cholesky decomposition of the disturbance covariance structure using the Durbin-Levinson algorithm \citep{Levinson1946,Durbin1960isr} and pre-whiten the model to obtain serially uncorrelated disturbances. Estimation proceeds by minimizing the residual sum of squares of the pre-whitened model, corresponding to what is also known as feasible generalized least squares. Complications arise from the nonlinear regression function. Our main contribution is the development of a neural network framework that is able to jointly estimate the unknown regression function and ARMA parameters.
We conduct an empirical study of cloud cover prediction, an important factor in climate projections. Using daily measures of cloud cover for a grid of locations within and around the Mediterranean Sea (Southern Europe and Northern Africa), the proposed NNARMA model leads to more accurate predictions than existing methods, including the long short-term memory (LSTM) network of \cite{Hochreiter1997}. The empirical results indicate improvements across locations, with particularly pronounced gains in mountain areas relative to linear models with ARMA errors, consistent with stronger nonlinear effects in such regions.
Simulation evidence suggests that the proposed estimation procedure accurately recovers regression functions of varying complexity and the underlying ARMA dynamics across a range of time-series lengths and signal-to-noise ratios. The finite-sample properties of the estimation procedure and the out-of-sample predictive performance of the NNARMA model demonstrate robustness to ARMA misspecification, and the latter further demonstrates robustness to dynamic misspecification induced by omitted lagged variables. Both remain stable except under substantial underspecification of the ARMA disturbance structure or the network architecture.
The remaining part of the paper is structured as follows. Section \ref{sec: literature review} discusses the literature on regression with serially correlated disturbances and our contribution. In Section \ref{sec: methodology}, we present our methodology for neural network regression with ARMA errors, including estimation and model selection. Section \ref{sec: Monte Carlo (main text)} demonstrates the finite-sample properties and out-of-sample predictive performance of our methodology, including under various forms of misspecification, by summarizing key results from extensive Monte Carlo experiments reported in Appendix \ref{sec: Monte Carlo (appendix)}. In Section \ref{sec: CFC application}, we present our empirical study of cloud cover prediction. Section \ref{sec: Conclusion} briefly summarizes and concludes.
\section{Nonlinear regression with serially correlated disturbances}\label{sec: literature review}
Consider the univariate time series process for some variable $y$ as denoted by $\{y_t\}_{t\in\mathbb{Z}}$, for time index $t$ with regular spacing. We define the nonlinear regression model with serially correlated disturbances for $y_t$ as
\begin{equation}\label{eq: model intro}
y_t = f(x_t;\beta) + u_t,
\end{equation}
where $f(\cdot)$ is the scalar nonlinear regression function, $x_t$ is a column vector of known regressors (explanatory variables), $\beta$ is a column vector of regression coefficients, and the disturbance sequence $\{u_t\}_{t\in\mathbb{Z}}$ is a weakly stationary process with zero mean, $\mathbb{E}(u_t)=0$, finite variance, $\text{Var}(u_t)=\omega ^2>0$, and serial correlation $\text{Cov}(u_t, u_{t-h}) = \gamma _h$, for $h=1,2,3,\ldots$, where $\omega ^2$ and $\gamma _h$ can be treated as unknown coefficients. For notational convenience, we omit an intercept in \eqref{eq: model intro}; its inclusion is straightforward. In the linear case, re-centering $y_t$ and each element in $x_t$ provides an equivalent representation.
\subsection{Linear regression with ARMA disturbances}
It is common practice for the statistical treatment of model \eqref{eq: model intro} to impose parametric restrictions on both the functional form of $f(\cdot)$ and the dynamics of the serially corrrelated disturbances $u_t$. For example, one may consider model \eqref{eq: model intro} with
\begin{equation}\label{eq: model linear}
f(x_t;\beta) = x_t\, ' \beta ,
\hspace*{0.8cm} u_t = \phi u_{t-1} + e_t,
\end{equation}
where both $x_t$ and $\beta$ are column vectors, $\phi$ is an autoregressive coefficient, and $\{e_t\}_{t\in\mathbb{Z}}$ is a white noise sequence of innovations with zero mean, $\mathbb{E}(e_t)=0$, and finite variance, $\text{Var}(e_t)=\sigma ^2>0$. For this specification, model \eqref{eq: model intro} reduces to a linear regression with autocorrelated disturbances of order one, and with $\omega ^2= \sigma ^2 / (1-\phi ^2)$ and $\gamma _h = \omega ^2 \, \phi ^h$. To ensure a stationary disturbance process $\{u_t\}$, one typically imposes the restriction $|\, \phi \, | < 1$ such that both $\omega ^2<\infty$ and $\gamma _h<\infty$, for all $h$. The resulting linear regression model with serially correlated errors is the framework from which the Durbin-Watson test is developed for the null hypothesis $H_0 : \phi = 0$ \citep{Durbin50,Durbin51}. The best linear unbiased estimate (BLUE) of $\beta$, for a given and known value of $\phi$, can be obtained by least squares after pre-whitening the disturbances using the data transformations
\[
y_t^\ast = y_t - \phi y_{t-1},\hspace*{0.8cm}
x_t^\ast = x_t - \phi x_{t-1},\hspace*{0.8cm}
t=2,3,4,\ldots .
\]
The least squares method applied to $y_t^\ast=x_t^\ast \beta + e_t$ delivers the BLUE of $\beta$.
In this paper, we consider the autoregressive moving average process, denoted by ARMA$(p,q)$, for the disturbance sequence $\{u_t\}_{t\in \mathbb{Z}}$. It is specified as
\begin{equation}\label{eq: model ARMA}
u_t = \phi _1 u_{t-1} + \ldots + \phi _p u_{t-p} + e _t + \omega _1 e_{t-1} + \ldots + \omega _q e_{t-q},
\end{equation}
with autoregressive coefficients $\phi _1, \ldots , \phi _p$ and moving average coefficients $\omega _1,\ldots,\omega _q$. Particular stability conditions can be imposed on the autoregressive coefficients to ensure a stationary process $\{u_t\}$ \citep{hamilton1994series}. The least squares method applied to the regression model \eqref{eq: model linear} with ARMA disturbances from \eqref{eq: model ARMA} is originally treated by \cite{Pierce1971}. The exact maximum likelihood approach of \cite{HarveyPhillips1979} entails the joint estimation of $\beta$ and the ARMA coefficients by adopting data transformations and using the Kalman filter to numerically compute the likelihood function. We adopt a similar approach but rely on the pre-whitening transformations carried out by the Durbin--Levinson algorithm, which is tailored for the case of ARMA disturbances. Also, we adopt this approach for nonlinear neural network regression.
\subsection{Parametric and nonparametric nonlinear regression}
The nonlinear regression model \eqref{eq: model intro} with a known parametric functional form for $f(\cdot)$ and serially uncorrelated disturbances $\{u_t\}$ is treated in textbooks such as \cite{Gallant1987} and \cite{BatesWatts1988}, which establish the asymptotic theory of the nonlinear least squares (NLS) estimator. When disturbances follow an ARMA process \eqref{eq: model ARMA}, NLS remains consistent but is inefficient. This inefficiency has motivated likelihood-based methods that jointly estimate the nonlinear regression function $f(\cdot)$ and the ARMA parameters \citep{harvey1990, hamilton1994series}. These approaches build on the pre-whitening ideas in \cite{HarveyPhillips1979} and yield efficiency gains and more reliable inference.
In practice, it can be challenging to justify a known parametric functional form for $f(\cdot)$, motivating the use of nonparametric and semiparametric regression methods. Kernel-based estimators such as the Nadaraya–Watson estimator \citep{Nadaraya1964, Watson1964} and local polynomial estimators, together with spline methods, have been shown to be consistent under mixing or other weak‑dependence conditions \citep{Robinson1983, Robinson1988, HardleTsybakov1997}; see also \cite{PaganUllah1999} and \cite{LiRacine2007}. However, as in the parametric case, ignoring serial correlation in the disturbances may lead to inefficiency. Accordingly, methods based on pre-whitening transformations have also been proposed in classical nonparametric regression settings \citep{Xiao2003, Su2006, LIU2010}.
\subsection{Neural network regression}
Neural networks have emerged as powerful tools for nonlinear regression in data rich environments. They can handle high-dimensional input vectors $x_t$, whereas classical nonparametric methods suffer from the curse of dimensionality. Early econometric work of \cite{WhiteDomowitz1984} established consistency of feedforward networks as flexible nonlinear regressors, while the universal approximation results of \cite{HORNIK1989} provided the theoretical foundation for their use in regression; see also \cite{White1989}. Extensions to time series settings have been considered, for example, by \cite{Chen1999}, who establish results under mixing conditions.
In the context of time series forecasting, \cite{Zhang2003} considers hybrid neural network methods for ARMA models, while
\cite{Terasvirta2006} show how neural network regressions can be embedded in variable selection methods. Finally, the research of \cite{sun2021} is most closely related to our study. It considers a neural network regression model with serially correlated disturbances. The key differences are that \cite{sun2021} focus on a system of multiple time series variables where single regressions take place on the other variables, assume autoregressive disturbances of order one, $u_t = \phi u_{t-1}+e_t$, for each regression, and rely on approximate methods for pre-whitening the data. We treat the single equation \eqref{eq: model intro} with output variable $y_t$ and input vector $x_t$, assume ARMA disturbances $u_t$ as in \eqref{eq: model ARMA}, integrate the Durbin-Levinson pre-whitening within the neural network method, and jointly estimate the nonlinear function and ARMA parameters.
As in other nonparametric techniques, to control for the bias-variance tradeoff, estimation relies on a set of smoothing parameters. For example, \cite{Altman, Hart1991, HERRMANN1992,RayTsay2007, Krivobokova2007,BRABANTER2018} discuss the choice of smoothing parameters in nonparametric regression with serially correlated disturbances. In our neural network approach, the bias-variance tradeoff is controlled by both the number of layers through which the data are transformed and the size of each layer. We find that early stopping \cite{Prechelt2012}, a regularization technique, can alleviate the need for extensive model selection, making it practical to use a sufficiently flexible network architecture.
Finally, an alternative to the explicit modeling of serial correlation in the disturbances of model \eqref{eq: model intro} is to assume that disturbances are uncorrelated after lagged values of the dependent variable and the regressors are included in the regression function $f(\cdot)$. In the linear framework, this strategy corresponds to autoregressive distributed lag (ADL) models \citep{HendryNielsen2007}. Within nonlinear and nonparametric settings, this strategy is standard in machine learning methods \citep{ChenTsay1993NAAR,Medeiros2021}. Modern machine learning methods such as recurrent neural networks and their specialized form, long short-term memory (LSTM) networks, are designed to exploit lagged information to capture complex nonlinear dynamics \citep{rumelhart1986,Hochreiter1997}. The LSTM variant has been shown to achieve state-of-the-art performance for a wide range of problems \citep{Jozefowicz2015, Schmidhuber2017} and thus serves as a natural benchmark in our empirical analysis.
\section{Neural networks with ARMA disturbances}\label{sec: methodology}
We consider a time series length $T$, an observed univariate time series $y_1,\ldots,y_T$ (output), a number of regressor variables $K$, an observed $K\times 1$ vector time series $x_1,\ldots, x_T$ (input), and a univariate series of disturbances $u_1,\ldots,u_T$. In matrix notation, our proposed NNARMA model is given by
\begin{equation}\label{eq: NNARMA vector}
y=f(X;\theta _2) + u,
\end{equation}
where
\[
y := (y_1,\ldots,y_T)', \quad X := [x_1,\ldots,x_T]', \quad f(X) := (f(x_1),\ldots,f(x_T))', \quad u := (u_1,\ldots,u_T)',
\]
such that $y,f(X),u\in \mathbb{R}^T$ and $X\in \mathbb{R}^{T\times K}$, and where $f(\cdot;\theta _2)$ represents a feedforward neural network with some set of parameters $\theta_2$. We assume that the disturbance sequence $u_1,\ldots,u_T$ is a sample from the weakly stationary ARMA process as specified in \eqref{eq: model ARMA} with set of parameters $ \theta_1= \{ \phi_1,\ldots,\phi_p, \omega_1,\ldots,\omega_q \} $ and innovation variance $ \sigma^2$. It follows that $u$ has zero mean and covariance matrix $\text{Var}(u) = \Omega(\theta_1, \sigma^2) = \sigma^2 \Psi(\theta_1) $, where $ \Psi(\theta_1) $ is symmetric and has a Toeplitz (or band) structure. We compute a Cholesky decomposition of the scaled covariance matrix $ \Psi(\theta_1) = C(\theta_1)^{-1} [C(\theta_1)^{-1}]' $, where $ C(\theta_1) $ is a $ T \times T $ lower triangular matrix computed via the Durbin--Levinson algorithm; details are provided in Appendix \ref{sec: cholesky}. The (inverse) Cholesky factor $C(\theta_1)$ is used for pre-whitening the data. Pre-multiplying both sides of \eqref{eq: NNARMA vector} by $ C(\theta_1) $ yields a regression model with serially uncorrelated disturbances with mean zero and variance $\sigma^2$, we obtain
\begin{equation}
y^\ast = C(\theta_1)f(X;\theta_2) + u^\ast, \label{eq: C formulation}
\end{equation}
where $ y^\ast = C(\theta_1)y $ and $ u^\ast = C(\theta_1) u $, such that we have mean vector $\mathbb{E}(u^\ast)=0$ and covariance matrix $\text{Var}(u^\ast)=\sigma^2 I_T$. In case $f(\cdot)$ is linear, we simply have $\theta_2 = \{\beta\}$ and the matrix factor $C(\theta_1)$ can be pre-multiplied directly onto $X$ to obtain $X^\ast=C(\theta _1) X$. Then, the least squares estimator is obtained by minimizing the residual sum of squares in the pre-whitened regression model $ \lVert y^\ast - X^\ast \beta \rVert_2^2 $ with respect to $\beta$, subject to a given set of values for $\theta_1$. The joint estimation of $\theta_1$ and $\beta$ can be obtained by NLS, that is numerically minimizing the least squares objective with respect to $\{\theta _1,\beta\}$. Next, we extend this approach to the nonlinear case where $f(\cdot)$ is represented by a neural network.
For the transformed NNARMA model \eqref{eq: C formulation} where $ f(\cdot) $ is represented by a neural network, the complication arises that $f(\cdot)$ does not commute. However, we can exploit the compositional structure induced by the layered network architecture and implement the regression component in \eqref{eq: C formulation} through the scheme
\begin{equation} \label{eq: output NN}
C(\theta_1)f(X;\theta_2) = Z^{\ast (h)} \beta, \qquad Z^{\ast (h)} = C(\theta_1) Z^{(h)} ,
\end{equation}
where $\beta$ remains a regression coefficient vector and the $h\in\mathbb{N}$ hidden layers jointly construct the recursive system of matrix equations
\begin{align}
Z^{(h)} &= g \left( \kappa^{(h)} + Z^{(h-1)} \Gamma^{(h)} \right), \nonumber \\
&\vdots \label{eq: hidden last} \\
Z^{(1)} &= g \left( \kappa^{(1)} + X \Gamma^{(1)} \right),\nonumber
\end{align}
with element-wise nonlinear activation function $g(\cdot)$, intercept vector $\kappa ^{(\ell)}$ and coefficient matrix $\Gamma ^{(\ell)}$, for $\ell=1,\ldots,h$, and parameter set $\theta_2 = \{\beta, \kappa^{(1)},\ldots,\kappa^{(h)}, \Gamma^{(1)},\ldots, \Gamma^{(h)} \}$. The pre-whitening mechanism in \eqref{eq: output NN} is novel in a neural network context. Instead of applying the pre-whitening mechanism directly to $X$ as in the linear model, it is applied to the variables derived by the network $Z^{(h)}$ to obtain $ Z^{\ast (h)} $. The output of the network in \eqref{eq: output NN} is linear in $ Z^{\ast (h)}$, such that vector $Z^{\ast (h)} \beta \in \mathbb{R}^T$. The variables derived by the network constitute nonlinear transformations of the regressors $ X $, learned through a sequence of feedforward hidden layers in \eqref{eq: hidden last}. The output of each layer is a $ T \times m_\ell $ matrix $ Z^{(\ell)} $, $ \ell = 1,\ldots,h $, obtained from an affine transformation of its inputs followed by elementwise application of the nonlinear (activation) function $ g(\cdot) $. The function $ g(\cdot) $ induces nonlinear input transformations in each layer, introducing flexibility into the regression function and allowing the network to represent nonlinear relationships between $X$ and $y$. We use the Swish function, defined as $ g(z) = z \left(1 + \exp(-z) \right)^{-1} $ for $z\in \mathbb{R}$, which is a generalization of the widely used rectified linear unit (ReLU; \citealp{NairHinton2010}). It has been shown to perform well in comparison with ReLU for a range of problems \citep{ramach2018}. In each hidden layer, $ \kappa^{(\ell)} $ is a $ m_\ell $-dimensional vector of intercepts that is added row-wise, and $ \Gamma^{(\ell)} $ is a $ m_{\ell-1} \times m_\ell $ coefficients matrix. We follow the convention to set $ m_0 \equiv K $. The hyperparameter $ h\geq 1 $ denotes the number of hidden layers (the network depth), and $ m_\ell $ denotes the number of columns in $ Z^{(\ell)} $ (the size of the $ \ell $th hidden layer). The choice of network architecture, as characterized by its depth and layer sizes, is a standard model selection problem subject to a bias-variance tradeoff. Model selection is discussed in Section \ref{sec: Model selection (methodology)}.
To compute predictions $\hat y = f(X; \hat \theta_2) + \hat u $ based on the NNARMA model in \ref{eq: NNARMA vector}, neural network predictions of the regression component (without the pre-whitening mechanism) $f(X; \hat \theta_2) $ are combined with ARMA predictions of the disturbance component $\hat u$. The latter is conveniently obtained from the Durbin--Levinson algorithm applied to the neural network residuals.
\subsection{Parameter estimation}\label{sec: estimation}
To jointly estimate the set of ARMA parameters $ \theta_1= \{ \phi_1,\ldots,\phi_p, \omega_1,\ldots,\omega_q \} $ and the set of neural network parameters $ \theta_2 = \{\beta, \kappa^{(1)},\ldots,\kappa^{(h)}, \Gamma^{(1)},\ldots, \Gamma^{(h)} \} $, we minimize the residual sum of squares from the pre-whitened model in \eqref{eq: C formulation}, that is
\begin{equation}\label{eq: loss function}
\min_{\theta_1} \min_{\theta_2} \lVert y^\ast - C(\theta_1)f(X; \theta_2) \rVert_2^2 ,
\end{equation}
subject to two sets of constraints within $\theta _1$ given by
\begin{equation} \label{eq: stationarity condition}
1 - \phi_1 z - \ldots - \phi_p z^p \neq 0, \qquad 1 + \omega_1 z + \ldots + \omega_q z^q \neq 0, \qquad z \in \mathbb{C}, \qquad \lvert z \rvert \leq 1,
\end{equation}
which are imposed to enforce stability, for a stationary process $u_t$, and invertibility, for unique identification of $\omega_1,\ldots,\omega _q$, of the ARMA process in \eqref{eq: model ARMA}, respectively \citep{hamilton1994series}.
We solve the joint minimization problem by profiling out $ \theta_2 $ from the sum of squares. For a given $ \theta_1 $, an estimate of $ \theta_2 $ is obtained by numerically solving the inner minimization problem in \eqref{eq: loss function} using standard techniques from the neural networks literature,
\begin{equation}\label{eq: theta_2}
\hat{\theta}_2(\theta_1) = \arg \underset{\theta_2}{\min} \lVert y^\ast - C(\theta_1)f(X; \theta_2) \rVert_2^2.
\end{equation}
We rely on the popular Adam variant of gradient descent \citep{kingma2014adam} to obtain estimates from \eqref{eq: theta_2}. We use all data available in the estimation sample when evaluating gradients, sometimes referred to as batch learning. An estimate of $ \theta_1 $ is obtained by solving the outer minimization problem in \eqref{eq: loss function}, for the given choice of $ \theta_2 =\hat{\theta}_2(\theta_1)$ and subject to the constraints in \eqref{eq: stationarity condition}, we have
\begin{equation}\label{eq: theta_1}
\hat{\theta}_1 = \arg \underset{\theta_1}{\min} \lVert y^\ast - C(\theta_1)f(X; \hat{\theta}_2(\theta_1)) \rVert_2^2 .
\end{equation}
We use Jones' reparametrization of the parameters in $\theta _1$ to impose the stability and invertibility constraints in \eqref{eq: stationarity condition}, allowing \eqref{eq: theta_1} to be solved via unconstrained optimization; see Appendix \ref{sec: Jones} for details. We therefore rely simply on traditional (gradient free) numerical optimization techniques of \cite{Powell1964} to obtain \eqref{eq: theta_1} subject to the constraints implied by \eqref{eq: stationarity condition}. Finally, the estimate of $ \sigma^2 $ is obtained from
\begin{equation}
\hat \sigma^2 = \frac{\lVert y^\ast - C(\hat \theta_1)f(X; \hat \theta_2(\hat \theta_1)) \rVert_2^2 }{T}.
\end{equation}
We emphasize that the estimation procedure above solves the minimization problem in \eqref{eq: loss function} with respect to $\theta_1$ and $\theta_2$ subject to constraints \eqref{eq: stationarity condition}. It is \emph{not} an iterative step-by-step method such as the Cochrane--Orcutt procedure, where the estimate of $\theta_1$ is fixed when minimizing over $\theta_2$, after which the estimate of $\theta_2$ is fixed when minimizing over $\theta_1$.
The Adam algorithm is used to solve \eqref{eq: theta_2} and the Powell algorithm is used to solve \eqref{eq: theta_1} and they are both local optimization routines. The minimization problem in \eqref{eq: theta_2} is solved with respect to $ \theta_2 $ for different candidate choices of $ \theta_1 $. Each time, we initialize the slope coefficients in the neural network ($ \beta, \Gamma^{(1)},\ldots,\Gamma^{(h)} $) randomly from a Gaussian distribution, as suggested by \cite{He_2015}. Intercepts ($ \kappa^{(1)},\ldots,\kappa^{(h)} $) are initialized at zero, following standard practice in the neural networks literature \citep{Goodfellow-et-al-2016}. Initial values for $\theta_1$ are obtained by estimating the ARMA parameters as in \eqref{eq: theta_1} based on residuals from an initial neural network regression as in \eqref{eq: theta_2}, assuming serially uncorrelated disturbances for $u_t$, such that $\theta _1 \in \emptyset$ and $C(\theta_1)=I_T$.
\subsection{Early stopping}\label{sec: Early stopping}
We employ early stopping \citep{Prechelt2012} to mitigate overfitting and reduce the need for extensive model selection. In practice, early stopping reduces sensitivity to the choice of network architecture; see Section \ref{sec: Model selection (methodology)}.
The observation sample is split into an estimation and a validation sample. Let $T$ and $T_\text{val}$ denote the lengths of the estimation and the validation sample, with $T_{\text{tot}} = T + T_\text{val}$ denoting the length of the joint estimation and validation sample. The validation sample consists of observations $(X_t,y_t)$ for $t=T+1,\ldots,T_{\text{tot}}$, and follows the estimation sample in time to avoid information spillover due to serial correlation (look-ahead bias). We define
\begin{align*}
y_{1:T_{\text{tot}}} := (y_1,\ldots,y_{T_{\text{tot}}})', \qquad X_{1:T_{\text{tot}}} := [x_1,\ldots,x_{T_{\text{tot}}}]', \\
f(X_{1:T_{\text{tot}}}) := (f(x_1),\ldots,f(x_{T_{\text{tot}}}))', \qquad u_{1:T_{\text{tot}}} := (u_1,\ldots,u_{T_{\text{tot}}})'.
\end{align*}
Let $ C(\hat \theta_1)_{1:T_{\text{tot}}} $ denote the Cholesky factor of $ \text{Var}(u_{1:T_{\text{tot}}}) $. Denote by $ C(\hat \theta_1)_{T+1:T_{\text{tot}}} $ the submatrix formed by the last $T_\text{val}$ rows of $ C(\hat \theta_1)_{1:T_{\text{tot}}} $. After each iteration of the Adam algorithm used to estimate the parameters of the neural network $ \theta_2 $, early stopping monitors the loss on the validation sample based on candidate choices $\hat \theta_1, \hat \theta_2$,
\begin{equation}\label{eq: val loss}
L_{\text{val}} = \lVert C(\hat \theta_1)_{T+1:T_{\text{tot}}} ( y_{1:T_{\text{tot}}} - f(X_{1:T_{\text{tot}}}; \hat \theta_2) ) \rVert_2^2.
\end{equation}
Pre-multiplication by $ C(\hat \theta_1)_{T+1:T_{\text{tot}}} $ instead of $ C(\hat \theta_1)_{1:T_{\text{tot}}} $ in \eqref{eq: val loss} restricts attention to the validation sample. The Adam algorithm is stopped when no significant decrease in the validation loss $L_{\text{val}} $ is observed over 100 consecutive iterations, using a tolerance level of $ 10^{-6} $. We restore the estimate of $ \theta_2 $ associated with the lowest validation loss across iterations (not necessarily the last iteration).
\subsection{Model selection}\label{sec: Model selection (methodology)}
Our proposed NNARMA model requires choices for the number of hidden layers in the neural network $h$, the size of each layer $m_1,\ldots,m_h$, and the ARMA orders $p$ and $q$. Model selection can be performed based on in-sample criteria or out-of-sample predictive performances using time series cross-validation. However, repeated re-estimation across candidate network architectures and ARMA specifications quickly becomes computationally prohibitive. We therefore adopt a practical strategy that avoids such repeated estimation and is motivated by simulation evidence. Specifically, we fix a sufficiently flexible network architecture combined with early stopping and select the ARMA orders $p$ and $q$ by minimizing the Bayesian information criterion (BIC) applied to residuals from a preliminary neural network regression that assumes serially uncorrelated errors.
Simulation evidence demonstrates largely stable validation loss across network architectures, with mostly very narrow networks tending to produce higher losses (Appendix \ref{sec:Model selection MC Appendix}). With early stopping, performance is therefore largely insensitive to the precise choice of architecture once the network is sufficiently flexible. Further, it is shown that the ARMA specification selected by our practical strategy generally converges to the true specification as the sample size increases.
\FloatBarrier
\section{Monte Carlo experiments}\label{sec: Monte Carlo (main text)}
\FloatBarrier
In this section, we summarize the results from a comprehensive suite of Monte Carlo experiments. We examine the finite-sample properties of the estimation procedure in Section \ref{sec: estimation}, with early stopping as described in Section \ref{sec: Early stopping}, and evaluate the out-of-sample predictive performance of the NNARMA model relative to existing feedforward neural networks with lagged variables. We consider both correct specification of the NNARMA model and misspecification arising from the ARMA component as well as from omitted lagged variables (dynamic misspecification). Additional details on the Monte Carlo experiments, including the data-generating processes and supplementary results, are provided in Appendix \ref{sec: Monte Carlo (appendix)}.
\subsection{Monte Carlo setup and benchmark models}\label{sec: MC setup Main text}
For all experiments except those involving dynamic misspecification (Section \ref{sec:dynammic misspecification maint text}), we generate a univariate response variable $y_t$ as the sum of a single regression function $f(x_{1t},x_{2t})$ (a hump-shaped or a sinusoidal regression function, respectively) and a noise term. The regression input variables $x_{1t}$ and $x_{2t}$ are generated by stationary $\text{AR}(1)$ zero-mean processes, with autoregressive coefficients $0.8$ and $0.7$, respectively. The noise term is generated by a stationary $\text{ARMA}(1,2)$ zero-mean process \eqref{eq: model ARMA} with autoregressive coefficient $\phi _1 = 0.9$ and moving average coefficients $\omega _1 = -0.5$ and $\omega _2 = 0.2$.
We split each Monte Carlo sample into three consecutive parts: (1) an estimation sample (first $ 60\% $), (2) a validation sample used for early stopping (next $ 20\% $), and (3) a test sample reserved for out-of-sample evaluation (final $ 20 \% $). All models are estimated once, and we compute out-of-sample predictions conditional on the regressors (and lagged response variables) in the test sample. This setup mimics the empirical setting in Section \ref{sec: CFC application}.
For convenience, we denote by $\text{NNARMA}(p,q)$ the model proposed in Section \ref{sec: methodology} with $\text{ARMA}(p,q)$ disturbances. We denote by $\text{NN}(\ell_y,\ell_x)$ the feedforward neural network benchmark models with lagged variables, where $\ell_y,\ell_x \in \mathbb{N} $ denote the number of lagged dependent variables and regressors (including contemporaneous regressors), respectively. Under this notation, $\text{NNARMA}(0,0)$ and $\text{NN}(0,0)$ are equivalent. We consider 8 benchmark models: $\text{NN}(0,0)$, $\text{NN}(1,0)$, $\text{NN}(0,1)$, $\text{NN}(1,1)$, $\text{NN}(2,2)$, $\text{NN}(3,3)$, $\text{NN}(5,5)$, $\text{NN}(10,10)$.
We fix the network architecture to two hidden layers with $ 32 $ units in the first and $ 16 $ in the second, both for the NNARMA model and the benchmark models. This architecture is not intended to be optimal, but rather serves as a representative, sufficiently flexible choice for illustrating model performance without extensive tuning; see the discussion in Section \ref{sec: Model selection (methodology)}. We use the same network architecture in the NNARMA model and the benchmark models to ensure that differences in predictive performance reflect differences in model specification rather than architecture tuning. We denote the signal-to-noise ratio by $r$ and the sample size by $T$. Our results below are aggregated for both regression functions and all combinations of $ r \in \{0.05, 1.0\} $ and $ T \in \{250, 1000, 5000\} $ using 100 Monte Carlo replications. Appendix \ref{sec: Monte Carlo (appendix)} extends the analysis to $ r \in \{0.05, 0.1, 0.2, 1.0\} $ and $ T \in \{250, 500, 1000, 5000\} $ for a more comprehensive evaluation. These additional results do not alter the conclusions of this section.
\begin{figure}[t!]
\centering
\begin{subfigure}[c]{0.47\textwidth}
\centering
\includegraphics[width=\textwidth]{"new_mean_hump_shape_main_0_poster"}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.47\textwidth}
\centering
\includegraphics[width=\textwidth]{"new_mean_sine_waves_main_0_poster"}
\end{subfigure} \\ \vspace*{-0.45cm}
\begin{subfigure}[c]{0.47\textwidth}
\centering
\includegraphics[width=\textwidth]{"new_mean_hump_shape_main_1_poster"}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.47\textwidth}
\centering
\includegraphics[width=\textwidth]{"new_mean_sine_waves_main_1_poster"}
\end{subfigure} \\ \vspace*{-0.45cm}
\begin{subfigure}[c]{0.47\textwidth}
\centering
\includegraphics[width=\textwidth]{"new_mean_hump_shape_main_2_poster"}
\subcaption{Hump-shaped regression function}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.47\textwidth}
\centering
\includegraphics[width=\textwidth]{"new_mean_sine_waves_main_2_poster"}
\subcaption{Sinusoidal regression function}
\end{subfigure}
\caption{Average estimate of the hump-shaped (panel a) and the sinusoidal (panel b) regression functions. Subplots are arranged by $r$ (signal-to-noise ratio; columns) and $T$ (sample size; rows). The true function is shown transparently in the background.}
\label{fig: estimated reg func hump sine}
\end{figure}
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{"arma_hump_sine_same_lines_main"}
\caption{Sampling distribution of the ARMA coefficient estimates for the hump-shaped (solid) and sinusoidal (dashed) regression functions. Subplots are arranged by $r$ (columns) and $T$ (rows). Vertical dashed lines indicate the true coefficient values.}
\label{fig: ARMA coeffs main}
\end{figure}
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{"new_relative_mses_test_hump_sine_box_main_NNARMA_wp"}
\caption{Out-of-sample results for NN and NNARMA (correctly specified) models with Box plots of relative MSEs normalized against NNARMA, for the hump-shaped (in blue) and the sinusoidal (in red) regression functions, where values above one indicate that NNARMA has lowest error. Subplots are arranged by $r$ (columns) and $T$ (rows).}
\label{fig: OOS main}
\end{figure}
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{"new_relative_mses_arma_test_hump_sine_box_main_wp"}
\caption{Out-of-sample results for NNARMA misspecified models ($\text{NNARMA}(1,2)$ is correctly specified) models with Box plots of relative MSEs normalized against $\text{NNARMA}(1,2)$, for the hump-shaped (in blue) and the sinusoidal (in red) regression functions, where values above one indicate that NNARMA has lowest error. Subplots are arranged by $r$ (columns) and $T$ (rows).}
\label{fig: OOS main ARMA mis}
\end{figure}
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{"new_relative_mses_test_hump_sine_box_main_NN22_wp"}
\caption{Out-of-sample results for a selection of NN models ($\text{NN}(2,2)$ is correctly specified model) with Box plots of relative MSEs normalized against $\text{NNARMA}$, for the hump-shaped (in blue) and the sinusoidal (in red) regression functions, where values above one indicate that NNARMA has lowest error. Subplots are arranged by $r$ (columns) and $T$ (rows). The estimation procedure of the NNARMA model includes the optimal selection of $p$ and $q$.}
\label{fig: OOS main dyn mis}
\end{figure}
\subsection{Finite sample properties}\label{sec: Finite sample properties main text}
We start by examining the finite sample properties of the proposed estimation procedure, including early stopping, under the assumption that the NNARMA model is correctly specified.
Figure \ref{fig: estimated reg func hump sine} displays the NNARMA average estimate of the hump-shaped and sinusoidal regression function. For both regression functions, accuracy of the estimated surface is increasing with $r$ and $T$ as expected. For high $r$ (strong signal) and large $T$ (large sample), both functions are estimated with a high degree of accuracy. The hump-shaped function is estimated with reasonable accuracy even for low $r$ (noisy signal) and small $T$ (small sample), especially in the interior of the input region. All estimated surfaces exhibit reduced accuracy near the boundaries of the input region, a pattern that is particularly pronounced for the sinusoidal function and is well known in nonparametric settings \citep{Cattaneo2020} and neural network applications \citep{Jensen2023}. At low $r$, the sinusoidal function is estimated fairly well provided the sample size is large enough. It generally requires a larger sample size than for the hump-shaped function.
Corresponding to the surface estimates in Figure \ref{fig: estimated reg func hump sine}, Figure \ref{fig: ARMA coeffs main} displays the sampling distribution of the ARMA coefficient estimates.\footnote{Throughout, sampling distributions are estimated by kernel density estimation using a Gaussian kernel and the bandwidth suggested by \cite{Silverman}.} The behavior of the coefficient estimates closely mirrors that of the surface estimates: estimation accuracy is decreasing with $r$ and increasing with $T$. For low $r$ and large $T$, all coefficients are accurately estimated for both regression functions. For the hump-shaped function, estimation accuracy is high even for high $r$ and small $T$. For the sinusoidal function, the AR coefficient is accurately estimated for high $r$ and even for small $T$, although a larger sample size than for the hump-shaped function is required for the finite sample distribution of the MA coefficients to peak at their true values.
\subsection{Out-of-sample prediction performance}\label{sec: Out-of-sample main text}
We next extend the in-sample analysis using the same simulated data, and evaluate out-of-sample prediction accuracy of the correctly specified NNARMA model relative to the benchmark models with lagged variables.
Figure \ref{fig: OOS main} presents Box plots of relative out-of-sample mean squared error (MSE) for all benchmark models, normalized against the NNARMA$(1,2)$ model. To avoid distorting Figure \ref{fig: OOS main}, we omit the results for $\text{NN}(0,0)$ and $\text{NN}(0,1)$ because their MSEs are substantially higher, and also for $\text{NN}(10,10)$ its performance is always worse than specifications with fewer lags. Since almost all relative MSEs are above one, we can conclude that the NNARMA model almost always leads to more accurate predictions than the benchmarks, especially for the sinusoidal function. The improvements in accuracy of NNARMA relative to the benchmarks are increasing with $r$ and $T$.
\subsection{ARMA misspecification}\label{sec: ARMA misspecification main text}
We examine the effect of NNARMA misspecification on the finite-sample properties and out-of-sample predictive performances of the NNARMA models using the same simulated data set as above. Specifically, we consider all $16$ $\text{NNARMA}(p,q)$ model specifications with $p,q\in\{0,1,2,3\}$. Here, we focus on the out-of-sample performance. Appendix \ref{sec: MC Misspecification} presents a detailed discussion of both finite-sample properties and out-of-sample performance.
Figure \ref{fig: OOS main ARMA mis} presents Box plots of the relative out-of-sample MSE for all NNARMA misspecified models, except $\text{NNARMA}(0,0)$, normalized against the correct specification $\text{NNARMA}(1,2)$. The $\text{NNARMA}(0,0)$ specification yields substantially higher MSE values, especially for the sinusoidal function, and is therefore omitted to avoid distorting the graphs. Prediction accuracy is lower than that of the correct specification when the ARMA structure is substantially underspecified. In particular, the specifications $\text{NNARMA}(0,j)$, for $j=1,2,3$, produce less accurate predictions than $\text{NNARMA(1,2)}$. The predictions of $\text{NNARMA}(1,0)$ are also less accurate for the smaller samples, but are as accurate as those of $\text{NNARMA(1,2)}$ for $T=5000$, especially for the hump-shape regression function. The accuracy of $\text{NNARMA(1,1)}$ and all specifications where the ARMA structure is overspecified fluctuates around the same level as $\text{NNARMA(1,2)}$, for all combinations of $T$ and $r$.
\subsection{Dynamic misspecification}\label{sec:dynammic misspecification maint text}
In our final simulation experiment, we analyze the out-of-sample predictive performance of the NNARMA model when it is misspecified and the benchmark model $\text{NN}(2,2)$ is correctly specified. For this experiment, we generate a new data set for $y_t$ from the following data generation process,
\[
y_t = f(0.4x_{t}+0.3x_{t-1}+0.2x_{t-2}) + 0.5 y_{t-1} + 0.4 y_{t-2} + e_t, \qquad e_t \sim NID (0, 1/r),
\]
with $2\times 1$ regressor vector $x_t=(x_{1t},x_{2t})'$, for $t=1,\ldots,T$. We consider the same hump-shaped and sinusoidal regression functions for $f(\cdot)$ and the same values for $T$ and $r$. The simulated data set for $x_t$ remains the same as above. To analyze the effect of omitting lagged variables from the $\text{NN}(2,2)$ specification, we compare the out-of-sample MSEs to the model specifications $\text{NN}(1,0)$, $\text{NN}(1,1)$, and NNARMA. In the latter case, the NNARMA orders $p$ and $q$ are selected using our practical strategy discussed in Section \ref{sec: Model selection (methodology)}. For completeness, we also present the results for two overspecified models, $\text{NN}(3,3)$ and $\text{NN}(5,5)$. A detailed discussion of this simulation exercise is provided in Appendix \ref{sec: MC Misspecification Dynamics} together with results from a similar experiment where model $\text{NN}(0,2)$ is the correct specification.
The out-of-sample prediction accuracy for a selection of NN models is presented in Figure \ref{fig: OOS main dyn mis}. The MSEs are normalized against the NNARMA model optimized for each simulated series. The correct specification $\text{NN}(2,2) $ is more accurate than the optimized NNARMA model for high $r$ and large $T$, especially for the sinusoidal function with $r=1$ and $T=5000$. In this scenario, the correct specification $\text{NN}(2,2) $ is considerably more accurate than any other benchmark. For all other combinations of $T$ and $r$, the accuracy of NNARMA is comparable to that of the benchmarks, including the correctly specified $\text{NN}(2,2) $, with various instances of a higher accuracy for the NNARMA model. Appendix \ref{sec: MC Misspecification Dynamics} reports the model rankings by MSEs for each simulated series. For $r \in \{0.05, 0.1, 0.2\}$, the NNARMA model and the correct $\text{NN}(2,2)$ model show broadly similar frequencies of being the most accurate, especially for the sinusoidal function.
\FloatBarrier
\section{Cloud cover prediction}\label{sec: CFC application}
\FloatBarrier
Clouds play a crucial role in the climate system \citep{Hughes1984}. They are central to the global water cycle \citep{Bengtsson2010}, particularly through preciptation formation and distribution \citep{Stephens2005}. They also act as energy gatekeepers for the climate system by reflecting incoming solar radiation (cooling) and blocking outgoing terrestrial radiation (warming; \citealp{StephensETAl2012,LiuEtAl2023}). Future changes in cloud fractional cover (the proportion of the sky covered by clouds) are expected to affect global warming and, in turn, both climate and society \citep{Svennevik2024}. However, projections of future cloud cover are subject to substantial uncertainty \citep{Zelinka2020}, largely because the accurate representation of clouds remains a major bottleneck in numerical climate models, limiting the reliability of climate projections \citep{StevensBony2013,BonyEtAl2015,GrundnerEtAl2025,ThuyEtAl2025}. A more accurate representation of clouds in numerical climate models is therefore required. One solution is a statistical downscaling procedure based on regression as discussed by \citep{Svennevik2024}. We propose NNARMA as a candidate regression model, and show that it is useful for representing the nonlinear statistical relations between cloud cover and other climate variables and for improving cloud cover predictions.
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{"Topography"}
\caption{Topography map. Figure \ref{fig: TS plot 1} shows time series plots for the example locations labeled 1--5.}
\label{fig: topography map}
\end{figure}
\subsection{Data}\label{sec: empirical data}
We consider the European Cloud Cover data set from \cite{Svennevik2024}, available through the Open Science Framework \citep{DataOSF}. The data set comprises satellite observations of cloud fractional cover together with reanalysis fields for $(i)$ air temperature in Kelvin (K), $(ii)$ specific humidity expressed in kilograms of water vapor per kilogram of air (kg/kg), $(iii)$ relative humidity, measured by the amount of water vapor in the air as a percentage of the amount the air could potentially hold at that temperature and pressure, and $(iv)$ surface pressure in Pascals (Pa). Cloud fractional cover is the response variable in our analysis, and is sourced from the METeosat Second Generation (MSG) cloud mask from the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT; \citealp{Schmetz2002}). Sensor failures introduce missing values that are random in time but perfectly synchronous across all spatial locations. The remaining four variables serve as regressors and are sourced from the 5th Generation Reanalysis data (ERA5) from the European Centre for Medium-Range Weather Forecasts \citep{Hersbach2020}. They are retrieved from the surface or the closest pressure level (1000hPa). No missing values occur for the four regressor variables.
Cloud formation occurs when air becomes saturated with water vapor. Saturation can arise either from decreases in temperature (cooling) or from increases in water vapor. These mechanisms are reflected in the four regressors: $(i)$ temperature determines the saturation level; $(ii)$ specific humidity measures the amount of water vapor present; $(iii)$ relative humidity defines proximity to saturation; $(iv)$ surface pressure captures additional dynamical effects through its link to vertical air movement, with low surface pressure typically associated with rising air that cools and promotes condensation. Relative humidity is derived as a nonlinear transformation of temperature and specific humidity, and therefore may imply potential redundancy in the flexible neural network models. Since relative humidity encodes key thermodynamic relationships governing saturation, this variable helps in representing the resulting nonlinearities in the data \citep{Romps2014}. At the same time, relative humidity reflects only proximity to saturation, and cloud formation is influenced by thermodynamic structure in addition to relative humidity \citep{Slingo1987}.
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{"same_model_diff_locs_1_no_trunc"}
\caption{NNARMA predictions over the tests sample (red) together with observed cloud fractional cover (black) for the locations indicated in Figure \ref{fig: topography map}.}
\label{fig: TS plot 1}
\end{figure}
The geographical domain of the data set is between 30\degree \,and 50\degree \,latitude (degrees north) and -15\degree \,and 25\degree \,longitude (degrees east), covering Southern Europe and Northern Africa as shown in the topography map in Figure \ref{fig: topography map}. The data are provided at a 0.25\degree horizontal resolution with hourly sampling from April 2004 to December 2018. We average the data to daily frequency and sample these at a 1.0\degree horizontal resolution. The temporal aggregation eliminates most missing values, as days with at least one valid observation yield a complete daily record. In total, we consider $861$ geographical locations. At each location, the effective sample size is $T=5370$ observations. A total of $11$ observations are omitted due to missing values, and $7$ initial observations are reserved as initial conditions for benchmark models with lagged dependent variables, see the discussion below. The initial observations are also omitted in models without lagged dependent variables to ensure a common sample across specifications.
Appendix \ref{sec: Appendix supplementary figures}, Figure \ref{fig: Sumary plot appendix} summarizes the spatial distribution of the sample mean and standard deviation of cloud fractional cover and the regressors, computed over time at each location. Mean cloud fractional cover is higher over the eastern North Atlantic Ocean, particularly between 40\degree \, and 50\degree \, latitude, than over land and the Mediterranean Sea. Variability follows the opposite pattern. Over the ocean, a steady supply of water vapor and limited temperature variation maintain consistently high and stable relative humidity and cloud cover. Over land, a limited supply of water vapor reduces mean cloud cover, while strong day--night temperature fluctuations increase the variability of relative humidity and cloud cover. The Mediterranean Sea exhibits large variability in specific humidity due to alternating dry continental and moist marine air masses, with cloud cover showing intermediate mean and variability.
\begin{figure}[t!]
\centering
\begin{subfigure}[c]{\textwidth}
\centering
\includegraphics[width=\textwidth]{pq_bar_0_1_He}
\subcaption{Empirical distribution of ARMA specifications.}
\end{subfigure} \\ \vspace*{0.2cm}
\begin{subfigure}[c]{\textwidth}
\centering
\includegraphics[width=\textwidth]{pq_map3_0_1_He}
\subcaption{Spatial distribution of model classes.}
\end{subfigure}
\caption{Empirical distribution of ARMA specifications across the 861 geographical locations (panel a) and spatial distribution of model classes (AR, MA, or ARMA; panel b) selected by NNARMA and a linear model with ARMA errors.}
\label{fig: ARMA selections across locations}
\end{figure}
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{"main_same_loc_diff_models_loc_1_trunc"}
\caption{Predictions over the tests sample (red) together observed cloud fractional cover (black) for the location labeled 1 in Figure \ref{fig: topography map}.}
\label{fig: TS plot many 1}
\end{figure}
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{"linear_MSEDM_He_trunc"}
\caption{Relative MSEs over the test sample normalized to NNARMA (left), where values above one (red) indicate lower error for NNARMA. $p$-values from a two-sided Diebold-Mariano test of equal predictive accuracy to NNARMA (right), where values below $0.05$ (red) indicate significant differences in predictive accuracy at that level.}
\label{fig: Linear MSEDM}
\end{figure}
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{"LSTM_MSEDM_He_trunc"}
\caption{Relative MSEs over the test sample normalized to NNARMA (left), where values above one (red) indicate lower error for NNARMA. The color bar limits differ from Figure \ref{fig: Linear MSEDM}. $p$-values from a two-sided Diebold-Mariano test of equal predictive accuracy to NNARMA (right), where values below $0.05$ (red) indicate significant differences in predictive accuracy at that level.}
\label{fig: LSTM MSEDM}
\end{figure}
\subsection{Empirical design and benchmark methods}
For each of the $861$ geographical locations, we separately predict cloud fractional cover from the four regressors which are treated as given and known. We compare the prediction accuracy of our NNARMA model to that of the following benchmarks: a linear regression model with and without a lagged dependent variable $y_{t-1}$, as adopted in \cite{Svennevik2024}, a linear regression model with ARMA disturbances, and an LSTM neural network with and without a lagged dependent variable $y_{t-1}$. The recurrent structure of the LSTM network incorporates all prior input information when predicting cloud cover. Appendix \ref{sec: Appendix expanded set of benchmarks} presents the key figures from this section for an expanded set of benchmarks, which includes the same models but with seven lagged dependent variables instead of one. It is found that including more lags has a negligible effect on the results presented below. In all specifications, an explicit intercept term is included to align with the linear benchmark models.
For estimating the regression function and the parameters, including those for the ARMA part, we consider the same two-layer network architecture (with $ 32 $ and $ 16 $ units) as adopted in the Monte Carlo study, both for the NNARMA model and the LSTM network. We select the ARMA specifications for the NNARMA model and the linear regression with ARMA disturbances using the practical strategy discussed in Section \ref{sec: Model selection (methodology)}. We consider 36 $\text{ARMA}(p,q)$ specifications, for $p,q \in \{0,1,2,3,4,5\}$, as candidates.
We split the daily sample into three consecutive parts: an estimation sample (2004--2012), a validation sample used for early stopping (2013--2015), and a test sample reserved for out-of-sample evaluation (2016--2018). Only the NNARMA model and the LSTM network rely on early stopping. The linear benchmarks do not, and for these models we collapse the period from 2004--2015 into an extended estimation sample. All models are estimated once. No re-estimation is performed, and predictions are conditional on the regressors (and possibly lagged dependent variables) in the test sample. This design is motivated by \cite{Svennevik2024} and mimics the Monte Carlo design in Section \ref{sec: Monte Carlo (main text)}. We measure prediction accuracy by out-of-sample mean squared errors and use the \cite{Diebold1995} test to determine statistically significant differences across models.\footnote{Diebold--Mariano tests are computed using a Bartlett lag window and a truncation lag of $7$. Similar results are obtained with truncation lags in $\{1,2,3,4,5,6,14,21\}$ and are available upon request.}
The response variable $y_t$ lies in the unit interval $ [0,1]$. To avoid introducing additional model complexity, we follow \cite{Svennevik2024} and do not explicitly enforce this constraint through the loss function nor through structural constraints such as bounded link functions. All models are estimated with regular squared error loss. All regressors are standardized to have zero mean and unit variance. Ex post, predictions below 0 are set to 0 and predictions above 1 are set to 1. Appendix \ref{sec: Appendix supplementary figures}, Figure \ref{fig: Truncations box plot} shows box plots, across geographical locations, of the percentage of unrestricted out-of-sample predictions falling below 0, above 1, and outside the admissible range $[0,1]$. The median percentage of unrestricted predictions outside the admissible range is 1.6\% for NNARMA and 1.8\% or below for the remaining models. Thus, violations of the range constraint are rare in practice.
\subsection{Results}
Figure \ref{fig: TS plot 1} presents NNARMA out-of-sample predictions for the five illustrative locations indicated in Figure \ref{fig: topography map} together with the observed cloud fractional cover time series. These five locations span a representative range of land and ocean environments and altitudes that are present among the total $861$ locations in our study. Consistent with the data discussion above, these time series exhibit distinct temporal dynamics. The NNARMA model appears to capture much of this variation across locations.
The variation in temporal dynamics across locations is reflected in the heterogeneity of ARMA specifications selected for the NNARMA model and its linear counterpart, see Figure \ref{fig: ARMA selections across locations}, panel a.
Across locations, the most frequent selections are $\text{ARMA}(2,1)$ and $\text{ARMA}(2,2)$ for the NNARMA model and the linear model, respectively. The extreme candidate specifications $\text{ARMA}(0,0)$ and $\text{ARMA}(5,5)$ are never selected.
Figure \ref{fig: TS plot many 1} illustrates the NNARMA predictions in comparison with those of the benchmarks models for location 1. Most visibly, the LSTM network with no lagged dependent variable captures the least amount of variation in the observed series. The accuracy of the LSTM network increases substantially with the inclusion of a lagged dependent variable $y_{t-1}$. In contrast, the inclusion of $y_{t-1}$ appears to have a negligible effect on the accuracy of the linear model. Similar patterns are observed across the other locations; see Appendix \ref{sec: Appendix supplementary figures}, Figure \ref{fig: TS plot many 2} for location 2.
Figures \ref{fig: Linear MSEDM} and \ref{fig: LSTM MSEDM} summarize prediction accuracy across the entire geographical domain of the Mediterranean, by reporting relative MSEs for the benchmark models normalized against NNARMA, together with the $p$-values from two-sided Diebold-Mariano tests. Cloud formation is inherently threshold-based and nonlinear with respect to thermodynamic variables (e.g., \citealp{SundqvistEtAl1989, XuRandall1996, Romps2014}). Consistent with this nonlinear behavior, the NNARMA model provides significantly more accurate predictions than the linear model, with or without a lagged dependent variable, across most locations in the geographical domain. The use of ARMA disturbances in place of lagged dependent variables significantly improves the accuracy of the linear model; see also Appendix \ref{sec: Appendix supplementary figures}, Figure \ref{fig: lin vs lin}, which contrasts the linear specifications. This pattern is consistent with the ARMA selections in Figure \ref{fig: ARMA selections across locations}, panel b. Across most of the geographical locations, the selected ARMA specifications include moving-average dynamics, which are often important for accurate predictions but difficult to capture using a finite number of lagged dependent variables.
The NNARMA model often performs significantly better than the linear model with ARMA errors in mountain regions (compare Figures \ref{fig: Linear MSEDM} and \ref{fig: topography map}), characterized by low surface pressure and large variability in relative humidity (compare Figure \ref{fig: Linear MSEDM} and Appendix \ref{sec: Appendix supplementary figures}, Figure \ref{fig: Sumary plot appendix}). This pattern is consistent with how mountains modify airflow and temperature through terrain-induced vertical motion and associated thermodynamic changes (e.g., \citealp{SMITH1979,Houze2012}), which can amplify nonlinear effects, including threshold behavior, regime switching, and interactions. Mountain terrain forces air to rise and sink over slopes, with rapid cooling during ascent and rapid warming during descent. These processes generate large variability in relative humidity and repeatedly push air above and below the saturation threshold, creating abrupt, value-dependent changes in cloud-response sensitivity. Related to this, mountain regions exhibits strong regime switching because the terrain makes vertical wind motion highly sensitive to wind direction and atmospheric stability. Small changes in regional wind patterns or atmospheric stability can trigger a switch between cloud-forming and cloud-suppressing states, which may be better captured by the flexibility of the NNARMA model. In addition, mountains increase the importance of interaction effects. One example is the coupling of temperature and specific humidity. As air is forced to rise or sink over slopes, cooling or warming occurs without an immediate change in specific humidity, bringing the air closer to or farther away from saturation. As a result, temperature directly modulates the effect of specific humidity on cloud formation. Vertical motion also promotes mixing between air masses with different temperature and specific humidity. The NNARMA model naturally captures these interaction effects. In contrast, oceans typically experience smoother and more coherent atmospheric variability, leading to more gradual cloud responses and allowing the linear models to perform comparatively well.
Throughout most of the geographical domain, the NNARMA model achieves significantly higher predictive accuracy than the LSTM network without a lagged dependent variable. Including a lagged dependent variable in the LSTM network improves performance and leads to statistically indistinguishable accuracy relative to NNARMA across mainland Europe and large parts of the Mediterranean Sea, with each model occasionally outperforming the other. Across most of the remaining domain, including large portions of North Africa and the eastern North Atlantic Ocean, the NNARMA model continues to significantly outperform the LSTM network. In conclusion, incorporating ARMA structures can improve predictive accuracy in both linear and neural network models in this application of cloud cover prediction.
\FloatBarrier
\section{Conclusion}\label{sec: Conclusion}
In this paper, we have introduced the NNARMA model, which integrates autoregressive moving average (ARMA) error structures into feedforward neural networks for flexible nonlinear regression and prediction with serially correlated disturbances. We have developed a method for the joint estimation of the neural network and ARMA coefficients. Simulation experiments have demonstrated accurate estimation of unknown regression functions and the underlying dynamics in the disturbances across a range of time series lengths and signal-to-noise ratios. An empirical study of cloud cover has shown improved prediction accuracy of NNARMA relative to existing methods, including LSTM networks, with particularly pronounced gains in mountain regions compared to linear models with ARMA errors.
Our findings underscore the complementary roles of flexible nonlinear regression and structured stochastic modeling of the disturbances within the proposed framework. Explicitly representing temporal dependence through ARMA disturbances may provide advantages relative to the inclusion of lagged variables in the regression component in settings characterized by structured temporal dynamics. The simulations further indicate that ARMA disturbances can partially mitigate the impact of omitted lags in the regression equation, suggesting that the disturbance-based representation of temporal dependence can accommodate certain forms of dynamic misspecification. The neural network regression component enables the model to accommodate complex relationships that may not be well captured by linear specifications. The pronounced differences observed in mountain regions provide a concrete illustration of this mechanism. Together, these elements illustrate the potential of integrating structured stochastic modeling into flexible machine learning frameworks such as neural networks. From a practical perspective, the NNARMA model shows robustness to the choices of ARMA specification and network architecture, with noticeable deterioration only under substantial underspecification.
While the proposed NNARMA model emphasizes flexible nonlinear regression and structured stochastic modeling of the disturbances, several aspects remain open for further development. In particular, a rigorous theoretical analysis of the statistical properties of the proposed estimation procedure under serial dependence has not yet been established, formal uncertainty quantification is not yet available, and interpretability of the estimated nonlinear relationships may warrant further investigation. These considerations reflect broader methodological challenges in reconciling flexible neural network methods with formal statistical inference, and with particular attention to dynamic stochastic systems. Further, as alternatives to ARMA structures for the disturbances, stochastic processes implied by hidden Markov models and spatial-temporal models can, for example, be considered. Continued efforts to strengthen the theoretical foundation and inferential tools for such advanced methods may further enhance their applicability in time series analysis and forecasting.
\clearpage
\bibliographystyle{chicago}
\bibliography{References}
\newpage