Abstract

Recently, data with complex characteristics such as epilepsy electroencephalography (EEG) time series has emerged. Epilepsy EEG data has special characteristics including nonlinearity, nonnormality, and nonperiodicity. Therefore, it is important to find a suitable forecasting method that covers these special characteristics. In this paper, we propose a coercively adjusted autoregression (CA-AR) method that forecasts future values from a multivariable epilepsy EEG time series. We use the technique of random coefficients, which forcefully adjusts the coefficients with and 1. The fractal dimension is used to determine the order of the CA-AR model. We applied the CA-AR method reflecting special characteristics of data to forecast the future value of epilepsy EEG data. Experimental results show that when compared to previous methods, the proposed method can forecast faster and accurately.

1. Introduction

Forecasting time series data predicts future values by discovering a set of rules or identifying patterns from past data. Linear regression models for forecasting time series such as autoregressive (AR), moving average (MA), and autoregressive moving average (ARMA) are widely used [1]. However, these methods have difficulty obtaining accurate forecasts when the time series data has nonlinear characteristics that constantly change. Thus, soft computing techniques, such as fuzzy logic and neural networks, have been developed to resolve the problems of linear approach considering nonlinear properties and uncertainty of time series data [2, 3].

EEG time series signals obtained from a brain have irregular and complex wave structures. They also include a large amount of noise. Epilepsy EEG data is a representative example of a complex time series. Epilepsy is a disease defined by abnormal electrical activity in the brain that is central to the diagnosis of epilepsy. Epilepsy EEG signals display changes over time through constant interaction with external factors [4]. Noise is included within the complexity of epilepsy EEG data during measurements. Epilepsy EEG data is difficult to forecast because it has special characteristics, such as nonlinearity, abnormalities, and noise. Therefore, it is important to select an appropriate forecasting method because these characteristics affect the forecasting accuracy.

In recent years, studies have been conducted to automatically detect and predict epilepsy seizures using EEG data. Univariate, bivariate, and multivariate algorithms were proposed to solve the problem of seizure detection and prediction based on the EEG analysis of single or multiple electrodes [57]. Rabbi et al. applied nonlinear dynamics based on unvaried characteristic measurements to extract a correlation dimension from the intracranial EEG recordings and designed a fuzzy rule-based system for seizure prediction [8]. Iasemidis et al. proposed an adaptive seizure prediction algorithm (ASPA) based on the convergence of the short-term maximum Lyapunov exponents (STLmax) among critical electrodes in the preseizure phase [9]. Liu et al. introduced a seizure prediction approach using particle filtering [10]. Also, Shahidi Zandi et al. proposed a method to predict seizures by analyzing the entropy level corresponding to zero-crossing intervals in scalp EEG and its derivatives [11]. Many researchers also used autoregressive and spectral analysis for forecasting by extracting seizure precursors from the EEG data [12, 13]. However, these researchers used approaches that were based on the linear data or unvaried characteristics of data.

We can give the following problem definition.

Problem Definition. Given a time series data which includes special characteristics such as nonlinearity, nonnormality, and nonperiodicity, a forecasting model attempts to forecast the values over some future time period. More formally, given a time series of Epilepsy EEG , we forecast the value of . is the value of the time series at time , and is the forecasting length.

In this paper, we propose an adaptive forecasting algorithm that adjusts its coefficients of the autoregressive (AR) model forcedly. To forecast the future values of epilepsy EEG data including special characteristics, we use the random coefficients with −1 and 1 and the fractal dimension which the order of the CA-AR model determines. We conduct experiments with sets of EEG time series to evaluate the suitability of our forecasting approach. The experimental results demonstrate that the proposed method provides better forecasting performance than previous methods. The proposed algorithm provides the following benefits: seizure forecasting and warning to patients about seizures and actively probing the characteristics of seizure onset.

The remainder of the paper is organized as follows. In Section 2, we describe the proposed method for forecasting. The experiment results are compared with the other methods in Section 3. In Section 4, we discuss the related work on the prediction and analysis of seizures. Finally, we present our conclusions in Section 5.

2. Materials and Methods

An autoregressive model is a simple model to estimate the future value of a series using previous input values. The AR model represents a stochastic process using a general form of order as shown in (1): where is the observation at time , is the autoregressive coefficient of th order, and is white noise normally distributed with mean zero and variance at time .

An important feature of the AR model is utilizing recent past observations in the process of estimating the current observation at time . That is, the current observation can be estimated by a linear weighted sum of previous observations. The weights denote the auto-regression coefficients. The problem in AR analysis is the assumption that data is stationary and linear, and it must derive the best values for coefficients given a series . Several methods have been used to estimate AR parameters, such as Yule-Walker, least squares, and Burg’s method [1]. It has been shown that for large data samples these estimation techniques should lead to approximately the same parameter values. The Yule-Walker method applied the AR model to the signals by minimizing the forward forecasting error in a least squares sense. The unknown parameter is estimated as follows: Under the assumption that is normally distributed, this is also the maximum likelihood estimate of . The distribution of has been studied extensively. Unfortunately, the exact distribution of is unknown. Asymptotically, if , it has a normal distribution, while if , it is a Cauchy distribution. In addition, if , it is a nonstandard distribution [14]. These distributions can be used to approximate the finite sample distribution of . This suggests that the distribution would not adequately approximate the finite sample distribution, especially near the discontinuity point of , because the exact distribution of is continuous for all values of . It has been found that, unless is close to zero, these distributions do not approximate the distribution of finite samples well. The nonstandard limiting distribution when seems to give a good approximation of the finite sample distribution when is close to 1. However, it is too complex for practical use since an accurate approximation to this nonstandard limiting distribution can be obtained from the asymptotic expansion [15]. In this study, our aim is to find a forecasting method suitable for epilepsy EEG data. More specifically, suppose that is an independent and identically distributed sequence defined by where and . From (3), we have , , and [16].

In particular, if we take , we obtain a special case of an AR process [16]. We study this special case in this paper. As a motivation for the model in (3) for with , consider the order of the standard flexible coefficient AR model given in (1). If , we have (4) leading to the standard random model from (1): If we put , we have Both models (4) and (5) correspond to the standard (critical case) unit root [17]. A model that generalizes (4) and (5) is where the random coefficients take the values 1 or −1. The model in (6) can be viewed as a generalization of the standard random where the successive jumps go up or down according as or −1.

Epilepsy EEG data has special characteristics, such as nonlinearity and abnormal and nonstandard distributions [4]. Therefore, in this study, when EEG data has an abnormal distribution, we forcefully adjust the coefficients of AR. In this paper, the AR model is the basis of our coercively adjusted AR model (CA-AR). This model can be expressed as follows: where is random coefficients with and , with and for all and . has an important role in AR modeling since it determines the order of the coefficients. In the autoregressive model, to determine the order of the AR model is an important issue [18]. The order of an AR model, , must be appropriately selected because it determines the efficiency of the autoregressive model. If is smaller, then the estimation error is higher while calculation speed is faster. On the other hand, if is bigger, there are drawbacks requiring more computation time without any decreases in estimation error. Therefore, in order to resolve these drawbacks, an optimal way to determine the order of the AR model is required.

In this paper, the fractal dimension is used to determine the order of the CA-AR model. To calculate the fractal dimension, we apply the box-counting method [19]. The box-counting method is one of the most common methods to obtain the fractal dimensions using boxes that are big enough to cover the measured signal [20]. In other words, when the length of one side of the square is and the number of square boxes is , the box-counting dimension of is and . It can be expressed as the box-counting dimension of , , and the positive constant, : By taking logs on both sides of (8), we get Fractal dimension is given by (10): where is excluded while the denominator . Also, , and if is a negative number, will be a positive number. If the log-diagram of verses is a straight line, the fractal dimension is the slope of this straight line (as shown in Figure 3).

In this paper, the measured value using (10) is defined as the order of AR, , and it is applied to the CA-AR() model. In other words, to predict epileptic seizures from EEG data, the future value is predicted using , the observed values from the past. This paper will show that special case of epileptic seizure where is predicted by , where the series is AR() using (7). The estimated AR model is then used for prediction by applying the least squares estimation.

3. Results

In this section we present the empirical verification of our data analysis to forecast epilepsy EEG data. EEG datasets are provided in [4], and epilepsy EEG data set composed the five EEG datasets 5 (denoted by A~E). A and B datasets are recorded in the relaxed awake state of healthy volunteers (eye open or closed). C and D are measured during seizure free intervals, and E contained seizure activity. These five EEG datasets contain 100 single channel EEG segments of 23.6 sec duration, and they are sampled at 173.61 Hz. For our experiments, we used only three datasets such as A, C, and E.

3.1. Detection of Special Characteristics

In this paper, we proposed a novel approach to help in the improvement of epileptic seizure forecasting in nonlinear and nonperiodic EEG signals. In this section, we first analyze the characteristics of epilepsy EEG data which show nonlinearity and periodicity by applying cepstrum and lag plots. The cepstrum is employed to extract periodicities or repeated patterns [21]. The cepstrum analysis of a spectrum will have peaks corresponding to the spacing of the harmonics and sidebands. The -axis of the cepstrum shows frequency, and peaks in the cepstrum are related to periodicities. The cepstrum is employed to find the periodicity in Subjects A, C, and E. Figure 1(a) demonstrates the 50th original signal of Subject E recorded during seizure activity and Figure 1(b) displays the measured periodicity using the cepstrum. Figure 1(a) seems to have a periodic wave within the original signal. However, Figure 1(b) does not have any periodicity. In addition, Subjects A and C also do not show the periodicity.

Even though the periodicities in the original signal repeatedly appear as a sinusoidal wave during seizure activity, when we applied the cepstrum to the seizure activity signals, the results differ from the original signal. Seizure activity signals do not have any periodicity. We observed that our experimental results of the seizure activity signals by the cepstrum do not have any periodicity. Therefore, since most conventional forecasting or prediction approaches require periodicity in observed data, these approaches are not appropriate for the nonperiodic seizure activity signals.

In this paper we also applied lag plots to find hidden characteristics in the data. Lag plots are useful in the analysis of cyclical data [22]. A lag plot checks whether a dataset is random or not. In addition, they provide the autocorrelation of the data. Figure 2(a) shows the first raw signal of Subject A’s epilepsy EEG time series. The lag plot of Subject A is shown in Figures 2(b) and 2(c), where the lag and , respectively. Figure 2(b) shows a definite linear structure in the lag plot, which was hidden in Figure 2(a). That is, this lag plot exhibits a linear pattern. If the data is strongly nonrandom, we are able to apply an autoregressive model that might be appropriate for prediction. Figure 2(c) shows the Gaussian distribution of Subject A plotted with a lag of by plotting versus . Figure 2(d) is the first raw signal of Subject C. Figures 2(e) and 2(f) show lag plots for and , respectively. In the case of , linear patterns are shown for both Subject A and Subject C. Figure 2(f) shows similar results to Figure 2(c). However, Figures 2(h) and 2(i) for Subject E differ from Subject A and Subject C. In the case of shown in Figure 2(i), it is a mixture of Gaussian distribution. Generally, when a lag plot has a nonrandom pattern, the data can be predicted using conventional methods. However, the epilepsy EEG dataset shows random patterns that are difficult to be predicted by conventional approaches. Thus, we need a suitable method for forecasting of EEG signals whose characteristics are nonlinearity and nonperiodicity.

3.2. The Order of CA-AR Choice (Box-Counting) for Forecasting

In this section, we present how the order of CA-AR is determined using fractal dimension. We use box-counting analysis which is a common method for fractal dimension estimation. It is also known that it is easy, automatically computable, and applicable to patterns with or without self-similarity [20]. However, this technique, including the processing of data and definition of the range of grid size, requires proper implementation to be effective in practice. In this study, the grid size is changed from 0.1 to 10000 in multiplication of 2 [22]. The slope of the linear part of the plot is the estimated fractal dimension of the epilepsy dataset. In this method each signal is covered by a sequence of grids of ascending sizes. Two values are recorded for each of the grids: the number of square boxes intersected by the signal, , and the side-length of the squares, . The slope of the straight line formed by plotting against indicates the degree of complexity or fractal dimension between 1 and 2 [23]. A signal with a fractal dimension of 1 or 2 is considered as completely differentiable or very rough and irregular, respectively.

We measured the fractal dimension of the 100 single signals from each subject to determine the order of CA-AR using box-counting analysis. Figure 3 illustrates the plot of versus based on the grid sizes. The log-log plots of Figure 3 are used to estimate the fractal dimension that is computed from the slope of the plot. Figures 3(a), 3(b), and 3(c) are the log-log plot for the fractal dimension of the 1st signals of Subject A, Subject C, and Subject E. This graph clearly displays a horizontal straight line when the grid size is small or too big. Deviation from a linear straight line can be expected to lead to underestimation of the fractal dimension value for the skeleton of a signal.

We applied the box-counting method to estimate the fractal dimension of the Phase Space from a signal of each subject. The vector space of the delay coordinate vectors is termed the Phase Space [22]. The observation sequence is represented by the series , which gives the value of the time series at time . That is, we can define . is a real number greater than zero termed the time delay, and is any integer greater than zero. The vector is termed the delay coordinate vector, because its terms are the time-delayed data values from the time series. Given time series and lag , we form all the delay coordinate vectors from . The Phase Space is a -dimensional space.

Figure 4 demonstrates the estimated fractal dimensions. The -axis and -axis denote the Phase Space (time delay space) and slope, respectively. This implies that a lag length of one is sufficient to reconstruct the state space.

Figure 4 shows fractal dimension of few signal, and we confirm some particular results; if the dimension of signal increases, the fractal dimension (slope) increases. As a result of Figure 4, in the case of Subjects A and C, signals mostly have a fractal dimension between 5 and 7, while Subject E exhibits a fractal dimension between 2 and 4 when time delay is 20. That is, Subjects A and C exhibit an average slope between 4 and 5. However, in the case of Subject E, the seizure activity represents an average slope of 2.5. When the time delay dimension increases, the fractal dimension increases. However, the plot of fractal dimension versus lag length shows that fractal dimension does not significantly increase, as lag length is incremented. This experiment shows the determination of the value of parameters using log-log plots for time series prediction. That is, we select the round-up integer of the average slope values of all normal signals as Subjects A and C for the order of CA-AR, because we must predict abnormal behavior that dropped out of the past pattern. We ran the CA-AR model with the round-up integer “5” of the average of the slope of all signals for forecasting. To verify that the selected is the optimal order, we measured the forecasting error of each signal from each subject. That result can be confirmed in Section 3.3.

3.3. Forecast Accuracy

To evaluate the reliability of optimal order for our model, we measured Root Mean Square Error (RMSE) of forecasting from all signals of each subject. An autoregressive model of order implies that the current value of the time series is being predicted based on past th data of the same random variable. Thus, an autoregressive model of order can be expressed using the previous values of the time series. Let , be successive instances of the random variable , measured at regular intervals of time. We applied the standard AR model and CA-AR (coercively adjusted coefficient with or −1) model to forecast the new of the epilepsy EEG data.

We forecasted the signals from 481 to 500 time points by the proposed model and compared the forecasting errors between the optimal order that decided by the average of fractal dimensions and the other orders. For forecasting, we used the past values of time series, and is selected by fractal dimension. If , our model uses from 476 to 480 time points to forecast 481 time point. Table 1 shows RMSE results of several signals that were measured between the original values and the generated values by the model. RMSE is compared when the orders are 3, 5, 10, and 15 for AR and the proposed method. In the case of in the proposed model, RMSE has a higher accuracy than , , and in Subjects A, C, and E. In case of the standard AR, Subjects A, C and E exhibited the lowest RMSE in , similar to the proposed method. As shown in Table 1, the proposed method is exhibits a lower error rate than standard AR, and when is 5, the lowest RMSE exhibited in Subject E among the subjects. Thus, we confirm that optimal order is 5, and it is used as the order of the CA-AR during experiments to verify the efficiency. The order of the proposed method is determined with the round-up integer of the average fractal dimension that it is measured from all normal signals as Subjects A and C.

Figure 5 shows the result of the forecast snapshot of Subjects A, C, and E, using the CA-AR and standard AR models. Plots (a), (c), and (e) show a specific case of a 20-time step prediction of the 20th electrode, and plots (b), (d), and (f) provide the prediction result for the 80th electrode signal. The original signal is shown in Figure 5(a), from the time 481 to 500 in red line with the star point marker. The plus sign marker shows the forecasted signals by the proposed method (green line) and the point marker plots the forecasted signals by standard AR (blue line). These plots confirm that our forecasting method outperforms the conventional AR method.

In this paper, we compared the forecasting results among several existing methods and CA-AR method. Table 2 shows the forecasting results using existing methods of linear and nonlinear prediction (Artificial Neural Networks [3], Fuzzy [2, 24], Nearest Neighbors [25], and the proposed method) of the 7th electrode signal in Subject E. For the experiments, we used a single signal of each subject to measure the forecasting error in each forecasting method. Also, to measure forecasting errors of existing methods, we derived the partial training signals from original signals that are from the starting points to 500 time points or 1000 time points of signal. For example, when the existing method tries to forecast the length of 150 time points, 1 to 500 time points or 1000 time point of the signal are used for training. However, our method still uses only 5 time points to forecast the future values of 150 time points. CA-AR considers from the time 496 to 500 to forecast the future values of the time 501, and it used from the time 497 to 501 to forecast the time 502.

As a result of Table 2, the proposed method shows lower error rates than existing forecasting methods. In case of ANN, RMSE is gradually increased when the forecasting time points increase. However, our proposed method is unaffected by the length of the forecasting time points because it used only the order of AR. Besides, even though the length of forecasting time points increased, the error rate of the forecasting result did not show much change in our model. In case of the Fuzzy-based method, it performs the batch process using all of the previous values to estimate the new value. Thus, we quit measuring the forecasting error because it requires very long execution time.

Table 3 measures the forecasting time using the existing methods as in Table 2. For example, if the training time duration of ANN (Artificial Neural Networks) is 500 and the forecasting time duration is 150, then 0.5304 represents the execution time to forecast from the time 501 to 650. Fuzzy-based method [2, 24] was excluded from the 1500-forecasting time point test, because it already exhibited the highest forecasting time compared to the other methods in the 150-time point forecasting test. The proposed method achieves much faster forecasting time than ANN (Artificial Neural Networks) and NN (Nearest Neighbors). When the training time point or forecasting time point increases, the existing methods incrementally increase the forecasting time. However, since the proposed CA-AR method uses only a few observed values, it maintains a steady time. In the case of epilepsy seizure, we assume that it should be able to inform a patient a few minutes or several hours before the beginning of a seizure.

We evaluated forecasting error with each signal in each subject, and Table 4 shows the results. This experiment is done to measure the future values of the length of 500 time points using the observed signal from the time 1 to 500. That is, single signal has 4096 time points, and the existing method used from the time 1 to 500 of single signal for learning. Models forecast from the time 501 to 1000 through learning. In this paper, we measured RMSE between the original signals and the forecasted signals. As a result of Table 4, Subjects C and E indicated the lowest average forecasting error rate when NN algorithm was used. However, NN method required much learning time compared to our method. That is, the values of 1 to 500 in NN will be the training period for the forecasted value of 501. For the forecast of 502, the training period is 1 : 501. On the other hand, our method provides the fast computation time because its training period is 407 to 501 for the forecast of 502. Besides, the difference of error rate between NN and CA-AR is very small.

Our method needs only the past values of that are determined by fractal dimension to forecast the future values. Therefore, it guarantees the fastest computation time for learning to forecast the epileptic seizure. We measured the execution time of each model during forecasting the time length of 1000, 2000, and 3000. We can confirm the result of several signals from Table 5. In this experiment, we used each signal of Subject E and training signal used from the time 1 to 500 or 1000 time points of each signal in the existing method. As a result of Table 5, the lowest RMSE of forecasting is appeared in NN method except when CA-AR forecasted the length of the 2000 times using the length of the 1000 times as training signals. However, we need to look at the execution time of the whole forecasting method. NN method provided good forecasting results. However, it needs longer execution time. In addition, when the length of forecasting time increases, the execution time also increases.

The accuracy of time series forecasting is a very important factor to many decision processes, and hence the research for improving the effectiveness of forecasting models has lasted. Both the neural network and the AR model capture all of the patterns in the data [26]. Our method also can capture the patterns of data because it was based on the AR model. In epilepsy EEG data, the amplitude between normal and seizure signal presents a great difference. If a pattern of the generated signals by the proposed model deviates from the past pattern, CA-AR can regard these as the epileptic seizure. However, the proposed model does not separately determine or measure the sliding time window length to detect the change of pattern in this paper. The goal of this study is to provide the fast run time and high forecasting accuracy in time series data that has special characteristics such as the nonperiodicity and non-linearity. Through our experiments results, we can guarantee the fast execution time and accuracy between original signal and generated signal from our model.

4. Discussion

Epilepsy is a common neurological disorder in which some nerve cells spasmodically incur excessive electricity for a short time. Seizure predictions are mostly handled by statistical analysis methods from the EEG recordings of brain activity. The forecasting of epilepsy seizures can be used as a warning about seizures occurring on certain time scales by estimating the change in brain waves. That is, the forecasting of seizures alerts patients before an epilepsy seizure occurs. As a result, they could avoid potentially dangerous situations such as brain damage or injury during seizures.

In recent years, much research has looked into the prediction of epilepsy seizures using EEG data. Mormann et al. [27] analyzed bivariate EEG signals for seizure prediction. They analyzed the synchrony of EEG data using mean phase coherence (MPC) and maximum linear cross-correlation between EEG signal pairs. Schelter et al. [28] used MPC and obtained a proportion of seizures that were correctly predicted. Chávez et al. [29] analyzed the focal epilepsy EEG data for seizure prediction using non-linear regression analysis and phase synchrony. Winterhalder et al. [30] suggested the “seizure prediction characteristic” based on clinical and statistical considerations and compared to the performance of seizure prediction methods using concepts of linear and nonlinear time series analysis. This work indicates the uncertainty of predictions made by the use of the seizure occurrence period (SOP), in which the seizure is expected. However, it can be expressed when the independent variables are continuous. Moreover, these methods assume that the data has normality and independence.

Li and Yao [31] proposed prediction methods based on the wavelet transform and fuzzy similarity measurements of EEG data. This method is divided into two steps: to calculate the entropy of the EEG data and to calculate similarity between variables. Li and Ouyang [32] proposed the dynamical similarity measure based on a similarity index to predict epileptic seizures using EEG data. Gigola et al. [33] analyzed the time domain of different types of epilepsy to predict epileptic seizures using wavelet analysis based on the evolution of accumulated energy. Maiwald et al. [34] evaluated three nonlinear methods for seizure prediction: dynamical similarity index, correlation dimension, and accumulated energy. These methods can extract robust features from EEG data. However, part of these methods is performed based on the window unit and it is not sufficient for clinical applications.

Several techniques have been proposed to analyze characteristics of seizures via various methods. Liu et al. [35] measured the fractal dimension of the human cerebellum in magnetic resonance images (MRI) of 24 healthy young subjects using the box-counting method. Esteller et al. [36] determined the fractal dimension in the cortex electroencephalogram (IEEG, ECoG), using the Katz algorithm. Their results show that an electrographic seizure in the Electrocorticography (ECoG) occurs when there is an increase of complexity. Sackellares et al. [37] found that temporal lobe epilepsy is characterized by episodic paroxysmal electrical discharges (ictal activity). These discharges consist of organized synchronous activity of mesial temporal neurons, particularly those of the hippocampus. However, proper interpretation of such analyses has not been thoroughly addressed.

In this paper, we proposed a new CA-AR forecasting method based on the AR model that can forecast the seizure of complex epilepsy EEG data by applying the property of nonstandard distribution from [14]. The CA-AR model is suited to time series data with special characteristics, such as abnormality, noise, nonlinearity, and nonperiodicity.

5. Conclusions

Epilepsy may be caused by a number of unrelated conditions, including damage resulting from high fever, stroke, toxicity, or electrolyte imbalances. An algorithm capable of effective real-time epileptic seizure prediction will allow the patient to take appropriate precautions minimizing the risk of a seizure attack or injuries resulting from such an attack. Conventional methods for forecasting or prediction of data require periodicity in the observed data. However, when we applied the cepstrum, seizure activity signals did not exhibit periodicity. In addition, we could distinguish whether the epilepsy EEG data is random or nonrandom using the lag plot. If the lag plot has a nonrandom pattern, it can be used for prediction by conventional approaches. However, our data appears to have a random distribution.

This study proposed the random coefficients appropriate for random distribution data. Further, we used the log-log plots (box-counting) using the concept of fractal dimensions to forecast epilepsy EEG data to estimate the vital forecasting optimal order in our CA-AR model. Our experimental results demonstrate that CA-AR (coercively adjusted autoregressive) is the most suitable forecasting method for nonperiodic data. It does not require complex calculations and conducts fast forecasting compared to other methods. In addition, our method generates future values more accurately or similar than other methods. The experiments on epilepsy EEG data show that our method is not only fast and scalable but also accurate in achieving low prediction errors.

Future research could focus on extending CA-AR to perform forecasting on a multiple, coevolving time series which includes linear or non-linear correlations and periodicity or nonperiodicity. A more ambitious direction would be to automatically readjust the parameter and coefficient equations.

Acknowledgments

This research was supported by the MKE (the Ministry of Knowledge Economy), Korea, under the ITRC (Information Technology Research Center) support program (NIPA-2013-H0301-13-3005) supervised by the NIPA (National IT Industry Promotion Agency). It was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2012-0007810).