Earthquake forecasting (prediction) is known to be a challenging research program (problem) for which no single prediction method can claim to be the best and is hard to obtain the high and consistent prediction accuracy in worldwide. At large, earthquake data when viewed as a time series over a long time, exhibits a complex pattern that is composed of a mix of statistical features. This is due to the fact that the magnitude variation is high; hence the white noise component is dominating. A single prediction algorithm often does not yield an optimal forecast by analyzing over a long series that is composed of a large variety of features. In the past many attempts by using traditional time series forecasting algorithms have been made, and the prediction models are usually based on the input of univariate time series. In this thesis, a new analytic framework is proposed that allows these mixed features from the time series to be automatically extracted by a computer program, and fed into a decision tree classifier for choosing an appropriate method for the current forecasting task. The motivation behind this concept of a new analytic framework is to let the data decide which prediction algorithm should be adopted, adaptively across different periods of the time series. With a case study of global maximal earthquake prediction by different time series, the performance results obtained by traditional univariate time series forecasting methods, GMDH 1 alone, and extended GMDH with residual-feedback are generated and compared. The experiment results show that GMDH with residual-feedback can yield the lowest error of all relatively. It is forecasted that in merely one year from now, an earthquake of magnitude 10.2136 in Richter scales will strike the world, at the MAPE (mean absolute percentage error) of 2.211%. This case study is a novel approach that converts univariate time series to multivariate data is proposed, namely residual-feedback, for providing relevant multivariate inputs to GMDH or polynomial neural networks. This is important because the strength of GMDH just like any neural network is on predicting outcomes from multivariate data. GMDH is a well-known ensemble type of prediction method that is capable of modelling highly non-linear relations, and achieving an optimal accuracy by inducing through all possible structures of polynomial forecasting (prediction) models. The contributions of this thesis are: (1) a framework of automatic forecasting which is very suitable for real-time earthquake monitor is proposed, (2) an investigation on how different features of the data series are coupled with different prediction algorithms for the best possible accuracy, and (3) the case study of global maximal earthquake event prediction by different time series period by GMDH methods is significant and will potentially make a high impact to the development of forecasting techniques. In lieu of earthquake prediction, the same residual-feedback GMDH model could be used to forecast (predict) any kind of time series speculation with very low error measures.