Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
68,211 characters · 20 sections · 70 citation commands
Neural networks for nonlinear regression with serially correlated disturbances: Evidence from cloud cover
}
\affil[a]{Aarhus University, Denmark}
\affil[b]{Vrije Universiteit Amsterdam, The Netherlands}
{\let\relax}
Keywords: Neural networks; Time series; Nonlinear regressions; Serially correlated disturbances; Autoregressive moving average model; Durbin-Levinson; Cloud cover prediction. \\
Machine learning methods, including neural networks, are increasingly used in data science, econometrics and statistics 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 boxjenkins2015 for a textbook treatment. We focus on feedforward neural networks for a number of reasons. They are universal approximators HORNIK1989, Cybenko1989, LESHNO1993 and are known to perform well in practice on a number of regression and predictions tasks Kelly2020,Jensen2023. Theoretical justification is starting to be developed Farrell2021. 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 Goodfellow-et-al-2016. Bach2017 and Bauer2019 show that feedforward neural networks solve the curse of dimensionality that affects classical nonparametric methods 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 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 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) discusses the literature on regression with serially correlated disturbances and our contribution. In Section (ref), we present our methodology for neural network regression with ARMA errors, including estimation and model selection. Section (ref) 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). In Section (ref), we present our empirical study of cloud cover prediction. Section (ref) briefly summarizes and concludes.
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
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 (ref); its inclusion is straightforward. In the linear case, re-centering $y_t$ and each element in $x_t$ provides an equivalent representation.
It is common practice for the statistical treatment of model (ref) 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 (ref) with
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 (ref) 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$ 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
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\}$ hamilton1994series. The least squares method applied to the regression model (ref) with ARMA disturbances from (ref) is originally treated by Pierce1971. The exact maximum likelihood approach of 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.
The nonlinear regression model (ref) with a known parametric functional form for $f(\cdot)$ and serially uncorrelated disturbances $\{u_t\}$ is treated in textbooks such as Gallant1987 and BatesWatts1988, which establish the asymptotic theory of the nonlinear least squares (NLS) estimator. When disturbances follow an ARMA process (ref), 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 harvey1990, hamilton1994series. These approaches build on the pre-whitening ideas in 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 Nadaraya1964, Watson1964 and local polynomial estimators, together with spline methods, have been shown to be consistent under mixing or other weak‑dependence conditions Robinson1983, Robinson1988, HardleTsybakov1997; see also PaganUllah1999 and 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 Xiao2003, Su2006, LIU2010.
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 WhiteDomowitz1984 established consistency of feedforward networks as flexible nonlinear regressors, while the universal approximation results of HORNIK1989 provided the theoretical foundation for their use in regression; see also White1989. Extensions to time series settings have been considered, for example, by Chen1999, who establish results under mixing conditions.
In the context of time series forecasting, Zhang2003 considers hybrid neural network methods for ARMA models, while Terasvirta2006 show how neural network regressions can be embedded in variable selection methods. Finally, the research of sun2021 is most closely related to our study. It considers a neural network regression model with serially correlated disturbances. The key differences are that 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 (ref) with output variable $y_t$ and input vector $x_t$, assume ARMA disturbances $u_t$ as in (ref), 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, 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 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 (ref) 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 HendryNielsen2007. Within nonlinear and nonparametric settings, this strategy is standard in machine learning methods 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 rumelhart1986,Hochreiter1997. The LSTM variant has been shown to achieve state-of-the-art performance for a wide range of problems Jozefowicz2015, Schmidhuber2017 and thus serves as a natural benchmark in our empirical analysis.
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
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 (ref) 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). The (inverse) Cholesky factor $C(\theta_1)$ is used for pre-whitening the data. Pre-multiplying both sides of (ref) by $ C(\theta_1) $ yields a regression model with serially uncorrelated disturbances with mean zero and variance $\sigma^2$, we obtain
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 (ref) 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 (ref) through the scheme
where $\beta$ remains a regression coefficient vector and the $h\in\mathbb{N}$ hidden layers jointly construct the recursive system of matrix equations
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 (ref) 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 (ref) 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 (ref). 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; NairHinton2010). It has been shown to perform well in comparison with ReLU for a range of problems 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).
To compute predictions $\hat y = f(X; \hat \theta_2) + \hat u $ based on the NNARMA model in (ref), 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.
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 (ref), that is
subject to two sets of constraints within $\theta _1$ given by
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 (ref), respectively 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 (ref) using standard techniques from the neural networks literature,
We rely on the popular Adam variant of gradient descent kingma2014adam to obtain estimates from (ref). 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 (ref), for the given choice of $ \theta_2 =\hat{\theta}_2(\theta_1)$ and subject to the constraints in (ref), we have
We use Jones' reparametrization of the parameters in $\theta _1$ to impose the stability and invertibility constraints in (ref), allowing (ref) to be solved via unconstrained optimization; see Appendix (ref) for details. We therefore rely simply on traditional (gradient free) numerical optimization techniques of Powell1964 to obtain (ref) subject to the constraints implied by (ref). Finally, the estimate of $ \sigma^2 $ is obtained from
We emphasize that the estimation procedure above solves the minimization problem in (ref) with respect to $\theta_1$ and $\theta_2$ subject to constraints (ref). It is 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 (ref) and the Powell algorithm is used to solve (ref) and they are both local optimization routines. The minimization problem in (ref) 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 He_2015. Intercepts ($ \kappa^{(1)},\ldots,\kappa^{(h)} $) are initialized at zero, following standard practice in the neural networks literature Goodfellow-et-al-2016. Initial values for $\theta_1$ are obtained by estimating the ARMA parameters as in (ref) based on residuals from an initial neural network regression as in (ref), assuming serially uncorrelated disturbances for $u_t$, such that $\theta _1 \in \emptyset$ and $C(\theta_1)=I_T$.
We employ early stopping 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).
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
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$,
Pre-multiplication by $ C(\hat \theta_1)_{T+1:T_{\text{tot}}} $ instead of $ C(\hat \theta_1)_{1:T_{\text{tot}}} $ in (ref) 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).
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)). 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
\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), with early stopping as described in Section (ref), 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).
For all experiments except those involving dynamic misspecification (Section (ref)), 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 (ref) 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).
For convenience, we denote by $\text{NNARMA}(p,q)$ the model proposed in Section (ref) 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). 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) 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.
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) 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 Cattaneo2020 and neural network applications 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), Figure (ref) 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 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.
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) 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), 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$.
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) presents a detailed discussion of both finite-sample properties and out-of-sample performance.
Figure (ref) 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$.
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). 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) 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). 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) 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
\FloatBarrier
Clouds play a crucial role in the climate system Hughes1984. They are central to the global water cycle Bengtsson2010, particularly through preciptation formation and distribution Stephens2005. They also act as energy gatekeepers for the climate system by reflecting incoming solar radiation (cooling) and blocking outgoing terrestrial radiation (warming; 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 Svennevik2024. However, projections of future cloud cover are subject to substantial uncertainty Zelinka2020, largely because the accurate representation of clouds remains a major bottleneck in numerical climate models, limiting the reliability of climate projections 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 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.
We consider the European Cloud Cover data set from Svennevik2024, available through the Open Science Framework 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; 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 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 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 Slingo1987.
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). 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), Figure (ref) 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.
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 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) 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). 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 Svennevik2024 and mimics the Monte Carlo design in Section (ref). We measure prediction accuracy by out-of-sample mean squared errors and use the 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 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), Figure (ref) 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.
Figure (ref) presents NNARMA out-of-sample predictions for the five illustrative locations indicated in Figure (ref) 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), 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) 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), Figure (ref) for location 2.
Figures (ref) and (ref) 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., 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), Figure (ref), which contrasts the linear specifications. This pattern is consistent with the ARMA selections in Figure (ref), 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) and (ref)), characterized by low surface pressure and large variability in relative humidity (compare Figure (ref) and Appendix (ref), Figure (ref)). This pattern is consistent with how mountains modify airflow and temperature through terrain-induced vertical motion and associated thermodynamic changes (e.g., 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
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.