System Identification and Data‐Driven Forecasting of AE Index and Prediction Uncertainty Analysis Using a New Cloud‐NARX Model

Severe geomagnetic storms caused by the solar wind disturbances have harmful influences on the operation of modern equipment and systems. The modeling and forecasting of AE index are extremely useful to understand the geomagnetic substorms. This study presents a novel cloud‐nonlinear autoregressive with exogenous input (NARX) model to predict AE index 1 hr ahead. The cloud‐NARX model provides AE index forecasting results, with a correlation coefficient of 0.87 on the data of whole year 2015. The benchmarks on the data of the two interested periods of 17–21 March 2015 and 22–26 June 2015 are presented. The presented model uses uncertainty “cloud” model and cloud transformation to quantify the uncertainty throughout the structure detection, parameter estimation, and model prediction. The new predicted band can be generated to forecast AE index with confidence interval. The proposed method provides a new way to evaluate the model based on uncertainty analysis, revealing the reliability of model, and visualize the bias of model prediction.


Introduction
Many modern technological systems are sensitive to space weather disturbances, such as geomagnetic storms and substorms and ionosphere variability (Knipp, 2012;MacAlester & Murtagh, 2014;Schrijver, 2015). The severe situations of space weather can have harmful effect on power grid, navigation systems, and satellite system. Thus, it is extremely important to forecast the space weather disturbances to avoid damages and losses. In addition to the traditional first principle and statistical approaches for understanding the interactions between the solar wind and magnetosphere (e.g., Ala-Lahti et al., 2018, and the references therein), application of data-based methods and in particular techniques based on machine learning to the prediction of various geomagnetic indices resulted in many innovative forecasting models (e.g., Camporeale et al., 2018;Chandorkar et al., 2017;Wintoft & Cander, 2000).
The AE index, along with the Al and AU indices, was introduced by Davis and Sugiura (1966) as a measurement of global auroral electrojet activity (Mayaud, 1980). Changes in AE are driven by variations in the solar wind convection electric field produced by fluctuations in the solar wind velocity and interplanetary magnetic field (IMF). These two factors govern the efficiency of the coupling between the solar wind and terrestrial magnetosphere with the dominant role being played by a southward directed IMF. In this coupling process, the energy associated with the solar wind flow is converted into magnetic energy, which is transferred into the magnetosphere via reconnection processes on the dayside and is stored in the magnetotail. This energy is eventually released, energizing the plasma sheet, ring current, and ionosphere.
Three classes of interactions have been identified, depending upon the southward turnings of the IMF (see, e.g., Gonzalez et al., 1994). Short-lived southward turnings of the IMF with modest (Bz ≃ 3nT) give rise to minor intensifications of the ring current, yielding substorm events. Repeated southward turnings, referred to as high intensity, long duration, continuous AE activity events arise due to the occurrence of interplanetary Alfven wave train embedded within the solar wind flow (Tsurutani & Gonzalez, 1987). These events result in a continued period of AE activity. Finally, coronal mass ejections (CMEs) or magnetic clouds exhibit extended periods in which a strong Bz is observed. This coupling, between the CME and terrestrial magnetosphere, results in a major intensification of the ring current, and large deviations in both AE and Dst, and is referred to as an intense magnetic storm.
Previous studies of substorm using AE index have provided accumulated evidence that the magnetosphere behaves as a nonlinear dynamic system, and it can be described by a small number of variables (Kamide et al., 1998). There are plenty of studies aiming to forecast AE index from solar wind measurements. Among the many approaches of modeling and forecast, neural network (NN) is a commonly used method. Early in 1997, NN models were constructed to study prediction of the AE index (Takalo & Timonen, 1997). Later, an ANN algorithm based at IMF measured on Lagrangian point L1 and plasma measurements was introduced in 2008 to predict AE index from 5 to 60 min ahead (Pallocchia et al., 2008). The ANN models were further improved to achieve a correlation coefficient of 0.83 for 1-hr-ahead forecast and 0.80 for 3-hr-ahead forecast (Bala & Reiff, 2012). In addition, some other approaches are also applied for the analysis; for example, a correlation analysis with a technique of wavelet decomposition and selective reconstruction was applied to analyze the relationship between AE index and solar wind variables (Guarnieri et al., 2018). The advantage of NNs and its variants is that it can provide an efficient nonlinear representation to generate good model predictions. However, the identification process of NNs often involves a large number of variables, so that the model structure of NNs can be very complicated. From such model structure, it is quite difficult to further understand the nonlinear dynamic of the system, for example, which model term/variables are superior for describing the index and which model terms/variables are redundant. Nevertheless, it is obvious that such a model cannot provide a model structure that is simple and easy for understanding.
Another widely used approach for the modeling and forecast of magnetosphere is the nonlinear autoregressive with exogenous input (NARX) model. The NARX model is developed for the nonlinear system identification and can detect an appropriate model structure by selecting the most important model terms from a dictionary consisting of a huge number of candidate model terms (Billings, 2013). Thus, it is very efficient method for the space weather forecast due to the fact that the magnetosphere is a nonlinear process (Kamide et al., 1998). The NARX model have successfully solved the modeling and prediction of many magnetic indices, for example, the Dst index Boynton et al., 2011; and the Kp index (Ayala Solares et al., 2016). Comparing to the NNs, the NARX method only uses a small number of effective model terms to describe the system, so that the system can be represented a linear-in-the-parameter form, which is parsimonious and transparent. It is achieved by a model selection algorithm called orthogonal forward regression (OFR), which was initially developed as subset selection method for nonlinear modeling problem with unknown structure . In recent years, several variants have been introduced to improve the performance of NARX model and OFR algorithm, for example, the wavelet NARX model (Billings & Wei, 2005;, the iterative search algorithm with mutual information , and the common/robust model structure selection (RMSS) method (Gu & Wei, 2018;Li et al., 2016).
Under the assumption that the identified individual model structure elements can perfectly describe the true system components, most of the models are capable to provide accurate representations of the system. However, in many practical scenarios, this assumption is usually invalid due to the existence of uncertainty. Generally speaking, there are several types of uncertainty that may cause deleterious effect on the system modeling process. First, the uncertainty in data collection, for example, the experimental uncertainty and epistemic uncertainty, may lead to incomplete and inaccurate information in the data set. If the data samples are insufficient or some important variables are missing in the data set, it would be extremely difficult to find a model. Second, model type or model structure uncertainty directly affects the model performance. It is known that models are usually designed to represent some specific system features; that is, there are no single model type or structure that can perfectly describe all the true systems. Thus, it is essential to choose a correct model type to represent the system, and an inappropriate model type can largely reduce the model performance. Third, noise/disturbance is one of the main sources of uncertainty. The noise can be brought to the data through many ways, for example, measurement error from physical equipment, and external disturbances. The existence of noise could lead to biased parameter estimation and incorrected selected model terms. Therefore, quantifying uncertainty is essential for establishing the significance of findings and making predictions with known confidence.
From the literature, it is known that estimating the true uncertainty remains as an exclusive goal for AE index forecast. Under the effect of uncertainty, the identified model usually cannot perfectly represent the system but only approximately describe the system behaviors. In these situations, a single model may not always work well on future new data, as there might be a risk on trusting and relying on a single model for future system behavior forecasting. Thus, the model robustness becomes extremely important. Given the above concerns, this study focuses on the modeling and forecasting of AE index using NARX model structure with new uncertainty analysis method. In this study, a novel cloud-NARX model is presented to predict AE index 1 hr ahead and to quantify the uncertainty in the system modeling process. The reliability of the model can be quantified by the proposed uncertainty analysis, which makes the cloud-NARX model more robust than the conventional NARX model.
In summary, the main contribution of this work lies in the cloud-NARX model for (a) describing model structure and parameter uncertainty using a new uncertainty concept "cloud" model; (b) generating a new predicted band, which provides the confidence interval of predicted AE index; and (c) providing a new way to evaluate the model reliability based on uncertainty analysis.
The remainder of this paper is organized as follows. In section 2, the NARX model and cloud model are briefly reviewed. In section 3, the proposed cloud-NARX model is introduced. The AE and solar wind data are described in section 4. Section 5 presents the identified cloud-NARX model with model evaluation and uncertainty analysis. The study is summarized in section 6.

Construction of the NARX Model
The nonlinear autoregressive moving average with exogenous input (NARMAX) model structure can be described as where y(t) and u(t) are systems output and input signals; e(t) is a noise component that can be assumed to be white Gaussian (but such an assumption is not always necessary) in many applications. n y , n u , and n e are the maximum lags for the system output, input, and noise, respectively. F[ ] is some nonlinear function. A polynomial NARX model can be written as the following linear-in-the-parameters form: where φ m (t) = φ m (ϑ(t)) are the model terms generated from the regressor vector ϑ(t) = [y(t À 1), …, y(t À n y ), u(t À 1), …, u(t À n u )] T (T indicates the transpose of the vector), θ m are the unknown parameters, and M is the number of candidate model terms.
The OFR algorithm is briefly introduced as follows . First, the regression model and prediction error can be written in a compact matrix form: where be the initial dictionary of all the candidate model terms, the objective of OFR algorithm is to select a number of significant model terms to form a subset, which can be described as The output can then be described with the selected terms as follows: At first step of the term selection, the ERR index of each candidate model term of the initial dictionary can be calculated by where i = 1, 2, …, M. The first selected model term is the candidate model term with highest ERR value, as The first selected model term is φ l1 , and its associated orthogonal variable can be defined as q 1 ¼ φ l1 . The selected term φ l1 is then removed from the initial dictionary, and the dictionary D is then reduced to a subdictionary D M À 1 , which consists of M À 1 model candidates. The residual sum of squares can be calculated as At step s (s ≥ 2), the M À s + 1 bases are first transformed into new group of orthogonalized base q i with an orthogonal transformation as below: where q r (r = 1, 2, …, s À 1) are orthogonal vectors, φ j (j = 1, 2, …, M À s + 1) are the basis of unselected model terms of subset D M À s + 1 , and q s ð Þ j j ¼ 1; 2; …; M À s þ 1 ð Þ are the new orthogonalized bases. The rest of the model terms can then be identified step by step using the ERR index of orthogonalized subsets D M À s + 1 : The sth significant model term of the subset is φ ls , and its associated orthogonal variable can be defined as q s ¼ q s ð Þ ls . The residual sum of squares can be updated by : Recursively, the model terms of the subset φ l1 ; …; φ ln È É can be identified step by step, each at one step. By summing (11) for s from 1 to n yields The ‖r n ‖ 2 is called residual sum of squares, or sum squared error of the final model. The mean square error of the model can be calculated as ‖r n ‖ 2 /n. The number of the model terms that should be included in the NARX model can be decided by a modified generalized cross-validation criterion, also known as adjusted predicted sum of squares (APRESS). The APRESS is calculated as ) where p(k) is a penalty function defined in terms of the cost function C(k, α) = k × α with α being an tuning parameter.

Cloud Model and Cloud Transformation
Cloud model is a cognitive model that provides a way of bidirectional transformation between a qualitative concept "cloud" and the quantitative data "cloud drops" (Wang et al., 2014). The concept cloud is described by three numerical characteristics, namely, ex (expectation), en (entropy), and he (hyper entropy). Similar to normal distribution, ex is the expectation of all the elements in the set, and en is the variance of the distribution. he depicts the degree of departure from normal distribution of cloud model (Wang et al., 2014). Based on the theorem that any distribution can be represented by the sum of several normal distributions, the cloud model can be seen as an extension of normal distribution: When he equals to 0, the cloud model becomes actually a normal distribution. he is often regarded as an extra variable in practical situation, such as psychological quality of an athlete.
The bridges between cloud model and cloud drops is cloud transformation. The commonly used cloud transformation is the generic forward cloud transformation (GFCT) and generic backward cloud transformation (GBCT). The forward transformation is used to generate cloud drops from a known cloud model. The backward transformation is used to identify the cloud model from a sequence of cloud drops. In previous research, an ideal cloud backward transformation is also studied (Zhang et al., 2016). However, it is not feasible in real life as the groups of cloud drops could hardly be obtained in advance. The representation of GFCT and GBCT can be illustrated in Figure 1 and as follows: where the notation cloud (ex, en, he) represents a cloud concept of characteristics modeled from r samples numerical data cloud drops [x 1 , x 2 , …x r ]. The parameter ex, en, and he of the associated cloud model can be calculated from these samples. The GBCT can be described as follows (Wang et al., 2014). First, ex can be calculated by where u and v are the number of resampling folds parameters for GBCT u × v = r. Then calculate the sample variance of each group (i = 1, 2, …, v): Finally, en and he can be calculated as where Þ 2 are the sample mean and variance of γ i .
With an identified cloud model, a series of cloud drops can be generated with GFCT. Let ex, en, and he be the numerical characteristics and n be the number of cloud drops; let a and b be the parameters of number of folds of GFCT (a × b = z). The generated cloud drops x ij with certainty degrees μ(x ij ) (i = 1, 2, …, a, j = 1, 2, …, b) can be generated by several steps: First, generate a series of normally distributed random numbers γ i with expectation en and variance he 2 ; next, for each γ i , generate b normally distributed random numbers x ij with expectation ex and variance γ i 2 and calculate the certainty degree as The generic cloud transformation achieves the transformation between intension and extension of the cloud concept. The advantage of cloud model is that it provides a way to describe a distribution with only three parameters that cannot be characterized by traditional normal distribution. The cloud transformation is better and more powerful than normal distribution in that (i) it includes normal distribution as a special case and (ii) many data in real life do not follow a normal distribution. Based on cloud model and cloud transformation, a cloud-NARX model is proposed. The idea behind the metrics is to use a new uncertainty "concept" (cloud model) to replace the traditional model parameters.
A series of predicted points can be calculated by performing a transformation (generic cloud transformation) to generate a series of model parameters (cloud drops) from the concept. These predicted points form a predicted distribution (surface/band) with confidence intervals, describing the uncertainty and risk brought by model uncertainty. The cloud-NARX model can be described as: where Cloud li ex; en; he ð Þ i ¼ 1; 2; …; n ð Þ are the cloud models, which represent the estimated parameters and the uncertainty of these parameters. The parameters ex, en, and he are the characteristics of each cloud model.

Estimation of the Cloud-NARX Model
The estimation of cloud-NARX model consists of three steps, which are data resampling, submodel identification, and cloud parameter estimation. The general process of estimating the cloud-NARX model is shown in Figure 2.
First, the original data set can be regrouped to form K subdata sets through some resampling methods, for example, random sampling or bootstrap (Wei & Billings, 2009). Assume that the input and output sequence for the kth data set is u k ð Þ t ð Þ È É Nk t¼1 and y k ð Þ t ð Þ È É Nk t¼1 , respectively, for k = 1, 2, …, K. The model terms φ of the kth data set can be generated from the associated regressor vector relating to the kth data set [y (k) (t À 1), …, y (k) (t À n y ), u (k) (t À 1), …, u (k) (t À n u )] T . The submodel for the kth subdata set can be written in the compact matrix form Second, for each subdata set, the OFR algorithm can be applied to select a number of significant model terms to establish an individual NARX model. A common model structure can be formed by model terms selected in and important for all the subdata sets. In addition, a RMSS method is developed as an alternative method, for small-size data modeling problem (Gu & Wei, 2018 where i = 1, 2, …, n; in this way, the cloud-NARX model can be identified.
It is known that when the model structure is perfect, and the data are not corrupted with noises, any of the subdata sets will lead to exact the same model. However, the model structures of the submodels might be different when there is model uncertainty brought by the noises/disturbances/insufficient information. In these situations, any single model might be unreliable. By removing or adding some data points in the K subdata sets, the uncertainty can be quantified by the submodels with different structures and parameters.

Model Predicted Band and Averaged Prediction
With the identified cloud model of each parameter, K ' groups of cloud drops are generated using cloud forward transformation, as follows: is the generated parameters for the model term θ li with k ' = 1, 2, …, K ' . A number of K ' predicted time series of output y can then be calculated as where k ' is the index of predicted time series (k ' = 1, 2, …, K ' ). The K ' model prediction can then form a predicted band with density. The upper and lower boundaries of the predicted band can be determined as The averaged model prediction can also be calculated as

Model Performance Evaluation
To evaluate the averaged prediction of the model, the correlation coefficient (ρ), prediction efficiency (PE), and normalized root-mean square error (NRMSE) are calculated. The PE is defined as where σ 2 observed is the variance of the observed AE values and σ 2 error is the variance of the error between the predicted AE values and observed AE values. The accuracy of the predicted band can be defined as where N t is the total number of observed data points in test data set and N t ' is number of the observed data points within the predicted band.

10.1029/2018JA025957
Journal of Geophysical Research: Space Physics

Data
A full description of the solar wind variables and the magnetic indices is given in Table 1. The AE index is one of the most widely used indices for researchers in geomagnetism, aeronomy, and solar-terrestrial physics, to understand the geomagnetic activity. The AE index is the maximum deviation of the horizontal components of geomagnetic field variations from a set of globally distributed ground-based magnetometers located in and near the auroral zone in the Northern Hemisphere (Guarnieri et al., 2018). It increases when a substorm event is happening and represents the overall disturbance in both eastward and westward ionospheric electrojets located at around 100-km altitude (Davis & Sugiura, 1966).
The AE index and solar wind variables used in this study were all sampled hourly. The AE index and solar wind variables are used as the output and input of the systems modeling, respectively. The amplitude of the solar wind velocity is around 250-1,000 km/s, which is much larger than those of the other input signals. To avoid producing extreme parameter estimations, the solar wind speed/velocity variable is first normalized by V → V ' /1, 000, where V ' is the original signal and V is the normalized signal. Two derived variables, ffiffiffi p p and VBst = V × B T sin 6 (θ/2) , which are effective in describing the magnetic indices, are also used as input variables for the system model.

The Cloud-NARX Model
The AE and solar wind data from 2011 to 2013 (around 26,000 sampled data points) were used for training the model, and the data of 2015 (around 9,000 sampled data points) were used for model validation.  In order to determine the maximum time lags for both the input and output variables, we have carried out premodeling experiments and simulations; the results suggest that the maximum time lags of the input and output were chosen to be nu = 2 and ny = 2. The initial full model was chosen to be a polynomial form with nonlinear degree of 2. The input-output data points of training data set were first resampled 100 times with replacement, to form 100 subdata sets. For each subdata set, a NARX model with six model terms is identified. For convenience of description, these single NARX models are referred to as "individual NARX models." Thus, there are a total number of 100 different individual NARX models, and each has its own parameters. A total of 12 different model terms are selected during the 100 runs, and these terms are used for cloud-NARX model construction. The cloud parameters of each of these selected model terms are shown in Table 2.
It is noteworthy that the cloud-NARX model consists of 12 model terms, rather than six terms; this is because that each individual NARX model has its own structure. There are some common terms that are included in nearly all the individual NARX models, for example, VBst(t À 02) and y(t À 01). Also, some terms, for example, VBst(t À 02) × y(t À 01), are selected and included in a relatively small number of times out of the 100 runs. These rarely selected model terms are usually ignored in conventional NARX model because they make small contributions to the whole data set. However, in some of the subdata sets, they might be effective in some rare situations, for example, the peak times of the AE index. Figure 4 shows the normal cloud membership functions of the 12 selected model terms. The estimated parameters of the some model terms are normally distributed, for example, Bst(t À 02), V(t À 01) * Bst(t À 01), ð Þ , and V(t À 02) × p(t À 01). The distributions of the parameters Note. Bst = B T sin 6 (θ/2) (nT; Boynton et al., 2011).

10.1029/2018JA025957
Journal of Geophysical Research: Space Physics of some other model terms (e.g., VBs(t À 1), y(t À 1), VBst(t À 1) × y(t À 1)) are beyond normal distributions. Note that the normal distributions are not always sufficient to describe the distribution of the estimated parameters of these model terms due to the existence of uncertainty, which do not necessarily follow a normal distribution law. The three characteristics ex, en, and he are used to analyze the uncertainty of each model term. As discussed earlier, ex is the mean of estimated parameters of each model term, which is consistent with the conventional model parameter; en is the variance of the parameter estimation; and he is the hype-parameter to describe the degree of departure of the distribution to normal distribution. The values of en of some model terms (e.g., y(t À 1)) are quite small, which indicates that the parameters of these model terms in the individual models are very close. In other words, the contributions of these model terms are consistent in each individual model. On the contrary, the values of en of some model terms (e.g., VBs(t À 1)) are quite large, which means that uncertainty of the estimated parameters of these model terms are strong. In

10.1029/2018JA025957
Journal of Geophysical Research: Space Physics the disturbed periods, the contributions of these model terms are different in each individual model, and the width of the predicted band increases due to the prediction uncertainty. The cloud parameter he describes how much the distribution is beyond normal distribution. If the value of he is much smaller than that of en, it means the estimated parameters of the model term are normally distributed. With the hyper cloud parameter he, the cloud model can better describe the estimated model parameters that are not normally distributed.

One-Hour-Ahead Prediction of AE Cloud-NARX Model
As mentioned earlier, the cloud-NARX model is built on hourly sampled data, so the model can be directly used to generate 1-hr-ahead (i.e., one-step-ahead) predictions of AE index. With the cloud parameters, the generic cloud forward transformation was applied to generate 100 sets of model parameters (that is called cloud drops in the transformation) for all the selected terms. A total number of 100 time series of the AE index prediction were calculated. The average prediction and predicted band are presented in Figure 5. The predicted band À0.0795 0.2808 0.4866 VBs(t À 01) × VBst (t À 01) À5.6495 0.4665 0.7133 VBst (t À 02) × y(t À 01) À0.0195 0.0029 0.0008 Note. NARX = nonlinear autoregressive with exogenous input; ex = expectation; en = entropy; ℎe = hyper entropy. Bst = B T sin 6 (θ/2) (nT), VBst = V × Bst .

10.1029/2018JA025957
Journal of Geophysical Research: Space Physics is the quantification of uncertainty throughout the structure detection, parameter estimation, and model prediction. If the model structure is perfect and the parameters are estimated unbiased, the predicted band will be narrow. Otherwise, if there are strong uncertainty in the data itself or the model structure and parameter, the uncertainty will be propagated to model prediction, and the width of predicted band will be increased.
From Figure 5, the predicted band is very wide over the period from 17 Mar 2015 to 22 Jun 2015. This can be explained or understood as follows. First, from the input signals shown in Figure 3, we know that there were interplanetary disturbances over the two periods. It is known that in general most storms last quite a short period in the long-term evolution of the process. As a consequence, most of the training data were sampled at "quiet" times, and only a very small fraction of the training data is for the storm period. This results in that the training data are severely "imbalanced" (Ayala Solares et al., 2016). Therefore, while a single model may well capture the features and dynamics of the system at quiet times, it may not sufficiently capture the system dynamics at the severe situation times. That is why the prediction band is so wide for these stormy periods. Second, the wide prediction band over the period of 17 March 2015 and 22 June 2015 implies that no single model would produce reliable prediction of AE over stormy periods, no matter what/which method is used to build the model. That is why we propose to carry out uncertainty analysis in this study.
Note that the predicted band in Figure 5 provides only rough quantification of the uncertainty. In many situations, the detailed information of the predicted AE index at a specific time point is often needed. Figures 6 and 7 are the predicted bands with density over an 8-hr period on 17 March 2015 and 23 June 2015, respectively. These figures show the probability of the predicted AE index being in each interval. As shown in the two figures, the interval of the predicted band for each time point is divided into 100 bins. The histogram shows the probability (frequency) of a single predicted AE value occurs in each bin. The boundaries of the predicted band are also displayed with the histogram, to visualize the prediction uncertainty and make it easier to understand. In addition, it is straightforward to compare the observation (green line) and averaged prediction (blue line) of AE index in the figure. The overall accuracy of the predicted band on the test data set is 65%. The accuracy of high AE period (AE > 1000) is 70%.
The only way to reduce the width of the predicted band is to find a model structure that can better describe the true system. However, it is very hard, if not impossible, to obtain a perfect model structure for real-world system identification data modeling problem in the presence of strong uncertainty. Nevertheless, it should be noted that the performance of the model given by Table 2 outperforms previous models, for example,

10.1029/2018JA025957
Journal of Geophysical Research: Space Physics the NN model (Bala & Reiff, 2012; as shown in Table 3). Therefore, a wide predicted band might indicate that a severe situation (interplanetary disturbances) is likely to happen. The property of the predicted band could potentially be used to forecast the arrival of the interplanetary disturbances.

Performance and Advantage of the Cloud-NARX Model
The performance of the averaged prediction of cloud-NARX model is comparable to the best NARX model with very similar structure but fixed model parameters, as shown in Table 3. Figure 8 presents the scatter plot of the averaged prediction and observation. The correlation coefficient, PE, and NRMSE of the averaged prediction are 0.872, 75.97%, and 0.0589 (for data of year 2015), respectively, which are consistent with the best NARX model. The NARX model outperform the NN model for 1-hr-ahead prediction, as the previous NN model achieves the correlation of 0.83 (Bala & Reiff, 2012). More importantly, the cloud-NARX provides a transparent and parsimonious representation. As shown in Table 2, the NARX model reveals which of the variables/model terms are significant and which are not, for example, the model terms V(t À 02) × Bst(t À 01) indicates that the dayside reconnection 20-40 min prior (Balikhin et al., 2010) is an important component of the system, and the model terms y(t À 1) suggests that the autoregressive term has a significant effect on the AE index. On the contrary, the NN models are usually very complicated, and the training process involves a huge number of model terms and takes a lot of time.
The cloud-NARX model holds all the good properties of conventional NARX model and possesses an extra advantage. It provides a tool for understanding and analyzing uncertainty in the model structure and forecasting. For example, the uncertainty band in Figure 5   variables, e.g., Bst/ VBst). As discussed earlier, this property could potentially be used to forecast the arrival of a solar wind storm. Note that the model also works well and even better in nondisturbed periods. This is because that the model was trained on the data set where most of the data were sampled at nondisturbed period. Therefore, the system behaviors in nondisturbed periods are well captured by the identified model. A comparison between the observed and predicted AE index in two selected nondisturbed periods (23 April to 5 May and 19 October to 1 November) is given in Figure 9. According to the figure, the predicted band is narrow, which means that the uncertainty of the model is not strong. From these results, the cloud-NARX model also generates good prediction results in the nondisturbed times.
The model prediction of the cloud-NARX model and the conventional NARX model are consistent in nondisturbed periods. In some disturbed periods, the prediction performance of the cloud-NARX model is better than that of the NARX model. In addition, it is easy to generate long-term prediction using the cloud-NARX model. For example, the 3-hr-ahead AE index forecast can be  achieved by generating three-step-ahead model predicted output with the cloud-NARX model. The correlation coefficient, prediction efficiency, and NRMSE of the three-step-ahead model predicted output of the cloud-NARX model are 0.8167, 0.6667, and 0.0694, respectively. It is reasonable that the performance of 3-hr-ahead prediction is lower than that of the 1-hr-ahead prediction. It is because that at each step of the multiple-step-ahead prediction, the predicted AE index at previous step is used as the model input (as autoregressive variable). Thus, the prediction uncertainty of long-term prediction is increased due to the propagation of the error.

10.1029/2018JA025957
Journal of Geophysical Research: Space Physics

Conclusion
In this paper, a new cloud-NARX model was applied to the modeling and forecasting of AE index. Good forecasting results were obtained for 1-hr-ahead AE index prediction. The correlation coefficient between averaged prediction and observation is 0.87 and prediction efficiency of 0.81 when benchmarked for the period of 17-21 March 2015 and 22-26 June 2015, which is nearly identical to that produced by the best NARX model. The cloud-NARX model outperforms the previous models, for example, NNs. More importantly, the cloud-NARX model is capable to quantify the uncertainty of model structure, model parameter, and model prediction. The advantages of this new model can be summarized as follows. First, the model structure of cloud-NARX model is more robust than that of the conventional NARX model, as the model terms of cloud-NARX model are selected from resampled subdata sets. Second, the estimated parameters (ex, en, and he) of cloud-NARX model can provide more comprehensive information on the model parameter uncertainty. Third, based on cloud forward transformation, the cloud-NARX model can generate the predicted band, which clearly indicates the confidence interval of each predicted AE index. It is extremely important for space weather forecast, because when model becomes unreliable under some severe situations, the biased prediction could cause damages and losses. With the predicted band, the bias of model prediction can be identified, and the reliability of model can be evaluated.
One of the limitations of the paper is that there is still room for improvement for the accuracy of predicted band. Because for magnetic storm period and quiet times, the uncertainty is at different levels. Thus, it is highly desirable to further improve the cloud-NARX model to make it more adaptive for the high and low AE data.