# Integrated Feature Selection of ARIMA with Computational Intelligence Approaches for Food Crop Price Prediction

Authors: Yuehjen E. Shao; Jun-Ting Dai
Journal: Complexity
 (2018)

---

## Abstract

Because of global climate change, lack of arable land, and rapid population growth, the supplies of three major food crops (i.e., rice, wheat, and corn) have been gradually decreasing worldwide. The rapid increase in demand for food has contributed to a continuous rise in food prices, which directly threatens the lives of over 800 million people around the world who are reported to be chronically undernourished. Consequently, food crop price prediction has attracted considerable attention in recent years. Recent integrated forecasting models have developed various feature selection methods (FSMs) to capture fewer, but more important, explanatory variables. However, one major problem is that the future values of these important explanatory variables are not available. Thus, predictions based on these variables are not actually possible. Because an autoregressive integrated moving average (ARIMA) can extract important self-predictor variables with future values that can be calculated, this study incorporates an ARIMA as the FSM for computational intelligence (CI) models to predict three major food crop (i.e., rice, wheat, and corn) prices. Other than the ARIMA, the components of the proposed integrated forecasting models include artificial neural networks (ANNs), support vector regression (SVR), and multivariate adaptive regression splines (MARS). The predictive accuracies of ARIMA, ANN, SVR, MARS, and the proposed integrated model are compared and discussed. Experimental results reveal that the proposed integrated model achieves superior forecasting performance for predicting food crop prices.


---

## Body

## 1. Introduction

### 1.1. Background

Everyone needs food. However, not everyone has enough food to survive. Food crops are a primary source of human food, with rice, wheat, and corn being the most widely consumed sources of grains around the world. In this study, food crops refer to plants, which provide food for human consumption, cultivated by agriculturists. The demand and consumption for food crops will rapidly increase in the future, propelled by a 2.3 billion person increase in global population and greater per capita incomes anticipated through the midcentury [1, 2]. According to the Food and Agriculture Organization of the United Nations (FAO) [3], the demand for food is expected to grow substantially by 2050. A major factor for this increase is world population growth. Today, world population has reached 7.6 billion, and we may reach 9.7 billion by 2050. This growth, along with rising incomes in developing countries, is driving up global food demand. Consequently, humanity may directly face one difficulty, which is “will enough food crops be produced at affordable prices or will rising prices drive more of the humanity into hunger?”The three most important food crops, rice, wheat, and corn, directly provide more than 50% of all calories consumed by the world population. While the area harvested for wheat each year is 214 million ha, the areas harvested for rice and corn are 154 million ha and 140 million ha, respectively [4]. Additionally, the rapid rise in food prices has been a burden on the poor in developing countries, who spend roughly half their household income on food. The issue of being able to affordably purchase the aforementioned food crops plays a very important role in human life. Therefore, from the human point of view, the prediction of prices for rice, wheat, and corn has become a significant research topic.


### 1.2. Related Work

Owing to various agricultural and environmental factors, meteorological factors, and biophysical factors, it was indicated in [5] that the exports of rice were very difficult to predict. They employed the autoregressive integrated moving average (ARIMA) and artificial neural network (ANN) methodologies to predict Thailand’s rice exports. The wheat price in the Chinese market was predicted by using the ARIMA, ANN, and linear combination models [6]. The overall results showed that the ANN technique was the best prediction model. In Africa, around 40% of the rice consumed is imported [7]. This high dependence on rice imports indicates that Africa is highly exposed to international rice market shocks with sometimes grave consequences for its food security and political stability. In addition, models were constructed to predict the quarterly prices of two types of food crops, barley and wheat [8]. In [9], it was reported that the world food prices would rise around 32% by 2022.They showed the dynamic relationship between acreage response, food crop yield, and price volatility by developing an optimization model [10]. In [11], it was stated that most countries would like to predict their annual food requirements in order to provide food security for the people. They proposed an artificial intelligent support vector regression (SVR) model to predict the output energy in rice production.Because food crop prices are seasonally affected [5, 6], this study uses the ARIMA to predict the prices of rice, wheat, and corn. However, the assumptions inherent to the linear form of ARIMA may encounter problems in adopting nonlinear relationships for practical data [12]. Computational intelligence (CI) techniques have been widely used in many forecasting applications because of data-driven features and fewer a priori assumptions. Accordingly, in addition to ARIMA modeling, this study uses CI schemes, including ANN, SVR, and multivariate adaptive regression splines (MARS), for predicting the prices of the three food crops because they allow nonlinearity modeling and provide good forecasting characteristics.However, CI modeling may face difficulties in its training process for designing an optimal topology owing to the use of a high number of input variables. Therefore, feature selection methods (FSMs) have been incorporated in order to reduce the number of explanatory or predictor variables [13, 14]. In this study, feature selection refers to the process of identifying a subset of relevant explanatory or predictor variables for use during forecasting model construction. This subset of variables contains fewer but more important input variables that aid in predicting the outcome.Feature selection techniques have become a research hot spot in many forecasting applications. For example, in order to accurately predict wind speed, an SVR forecasting model with a feature selection procedure called phase space reconstruction has been proposed [15]. Additionally, a set of general ARIMA models was used for performance comparison. A hybrid forecasting model with a feature selection technique based on mutual information, extreme learning machines, and bootstrap techniques was proposed to predict day-ahead electricity prices [16]. The authors of [17] proposed a hybrid filter-wrapper approach to predict electrical load and the price of electricity. The feature selection method used in their approach identifies a minimum subset of the most informative features by considering the relevancy, redundancy, and interactions of candidate inputs. In [18], the authors proposed a two-step approach that selects a set of candidate features based on data characteristics. A hybrid ANN-based model was then used to predict the day-ahead price of electricity. A hybrid forecasting methodology was designed to predict electrical load in [19]. The study used a feature selection method that was developed based on entropy and several CI techniques. The prediction performance of the proposed model was compared to various ARIMA models. A regression model was proposed to predict online sales in [20]. Additionally, a feature selection methodology based on a multiobjective evolutionary algorithm was combined with the forecasting model. A combined ANN and Kalman filter approach was proposed to predict wind speed in [21]. An ARIMA feature selection method was used to determine the input structure for their hybrid model.When reviewing related research, we noticed that, although many studies have focused on using FSMs to obtain accurate forecasts, little attention has been devoted to the use of FSMs for food crop price prediction. Additionally, we encounter practical problems when using FSMs to predict food crop prices. As mentioned earlier, one of the major purposes for using an FSM is to select a subset of explanatory variables for further use in forecasting. However, the future values of the explanatory variables that are selected by the FSM are unknown. Relatively little research has attempted to address this problem. As a consequence, predictions cannot be computed even though the most important explanatory variables have been identified by effective FSMs. We propose an integrated forecasting model to remedy the problems caused by unknown explanatory variables. The proposed integrated model employs an ARIMA as an FSM to capture important self-predictor variables, which then serve as input variables for the ANN, SVR, and MARS models. Because the future values of a self-predictor variable can be computed based on its previous and current values, food crop price predictions can be obtained.In this study, we propose four single-stage models and three integrated models for predicting food crop prices. A real monthly dataset was obtained, which contains the prices of rice, wheat, and corn from January 1990 to September 2015. This real dataset makes it possible to compare predictions of food crop prices using single-stage models and integrated models. The remainder of the study is organized as follows: Section2 presents the modeling methodologies of different forecasting schemes. In Section 3, we describe the design structure for our single- and two-stage integrated models. Real data on the prices of three food crops are collected to verify the single-stage models and integrated models. Section 4 contains forecasting comparisons for four single-stage models and three integrated models. Some practical implications and other concerns are addressed in Section 5. The final section summarizes our research findings and contains our conclusions.


## 1.1. Background

Everyone needs food. However, not everyone has enough food to survive. Food crops are a primary source of human food, with rice, wheat, and corn being the most widely consumed sources of grains around the world. In this study, food crops refer to plants, which provide food for human consumption, cultivated by agriculturists. The demand and consumption for food crops will rapidly increase in the future, propelled by a 2.3 billion person increase in global population and greater per capita incomes anticipated through the midcentury [1, 2]. According to the Food and Agriculture Organization of the United Nations (FAO) [3], the demand for food is expected to grow substantially by 2050. A major factor for this increase is world population growth. Today, world population has reached 7.6 billion, and we may reach 9.7 billion by 2050. This growth, along with rising incomes in developing countries, is driving up global food demand. Consequently, humanity may directly face one difficulty, which is “will enough food crops be produced at affordable prices or will rising prices drive more of the humanity into hunger?”The three most important food crops, rice, wheat, and corn, directly provide more than 50% of all calories consumed by the world population. While the area harvested for wheat each year is 214 million ha, the areas harvested for rice and corn are 154 million ha and 140 million ha, respectively [4]. Additionally, the rapid rise in food prices has been a burden on the poor in developing countries, who spend roughly half their household income on food. The issue of being able to affordably purchase the aforementioned food crops plays a very important role in human life. Therefore, from the human point of view, the prediction of prices for rice, wheat, and corn has become a significant research topic.


## 1.2. Related Work

Owing to various agricultural and environmental factors, meteorological factors, and biophysical factors, it was indicated in [5] that the exports of rice were very difficult to predict. They employed the autoregressive integrated moving average (ARIMA) and artificial neural network (ANN) methodologies to predict Thailand’s rice exports. The wheat price in the Chinese market was predicted by using the ARIMA, ANN, and linear combination models [6]. The overall results showed that the ANN technique was the best prediction model. In Africa, around 40% of the rice consumed is imported [7]. This high dependence on rice imports indicates that Africa is highly exposed to international rice market shocks with sometimes grave consequences for its food security and political stability. In addition, models were constructed to predict the quarterly prices of two types of food crops, barley and wheat [8]. In [9], it was reported that the world food prices would rise around 32% by 2022.They showed the dynamic relationship between acreage response, food crop yield, and price volatility by developing an optimization model [10]. In [11], it was stated that most countries would like to predict their annual food requirements in order to provide food security for the people. They proposed an artificial intelligent support vector regression (SVR) model to predict the output energy in rice production.Because food crop prices are seasonally affected [5, 6], this study uses the ARIMA to predict the prices of rice, wheat, and corn. However, the assumptions inherent to the linear form of ARIMA may encounter problems in adopting nonlinear relationships for practical data [12]. Computational intelligence (CI) techniques have been widely used in many forecasting applications because of data-driven features and fewer a priori assumptions. Accordingly, in addition to ARIMA modeling, this study uses CI schemes, including ANN, SVR, and multivariate adaptive regression splines (MARS), for predicting the prices of the three food crops because they allow nonlinearity modeling and provide good forecasting characteristics.However, CI modeling may face difficulties in its training process for designing an optimal topology owing to the use of a high number of input variables. Therefore, feature selection methods (FSMs) have been incorporated in order to reduce the number of explanatory or predictor variables [13, 14]. In this study, feature selection refers to the process of identifying a subset of relevant explanatory or predictor variables for use during forecasting model construction. This subset of variables contains fewer but more important input variables that aid in predicting the outcome.Feature selection techniques have become a research hot spot in many forecasting applications. For example, in order to accurately predict wind speed, an SVR forecasting model with a feature selection procedure called phase space reconstruction has been proposed [15]. Additionally, a set of general ARIMA models was used for performance comparison. A hybrid forecasting model with a feature selection technique based on mutual information, extreme learning machines, and bootstrap techniques was proposed to predict day-ahead electricity prices [16]. The authors of [17] proposed a hybrid filter-wrapper approach to predict electrical load and the price of electricity. The feature selection method used in their approach identifies a minimum subset of the most informative features by considering the relevancy, redundancy, and interactions of candidate inputs. In [18], the authors proposed a two-step approach that selects a set of candidate features based on data characteristics. A hybrid ANN-based model was then used to predict the day-ahead price of electricity. A hybrid forecasting methodology was designed to predict electrical load in [19]. The study used a feature selection method that was developed based on entropy and several CI techniques. The prediction performance of the proposed model was compared to various ARIMA models. A regression model was proposed to predict online sales in [20]. Additionally, a feature selection methodology based on a multiobjective evolutionary algorithm was combined with the forecasting model. A combined ANN and Kalman filter approach was proposed to predict wind speed in [21]. An ARIMA feature selection method was used to determine the input structure for their hybrid model.When reviewing related research, we noticed that, although many studies have focused on using FSMs to obtain accurate forecasts, little attention has been devoted to the use of FSMs for food crop price prediction. Additionally, we encounter practical problems when using FSMs to predict food crop prices. As mentioned earlier, one of the major purposes for using an FSM is to select a subset of explanatory variables for further use in forecasting. However, the future values of the explanatory variables that are selected by the FSM are unknown. Relatively little research has attempted to address this problem. As a consequence, predictions cannot be computed even though the most important explanatory variables have been identified by effective FSMs. We propose an integrated forecasting model to remedy the problems caused by unknown explanatory variables. The proposed integrated model employs an ARIMA as an FSM to capture important self-predictor variables, which then serve as input variables for the ANN, SVR, and MARS models. Because the future values of a self-predictor variable can be computed based on its previous and current values, food crop price predictions can be obtained.In this study, we propose four single-stage models and three integrated models for predicting food crop prices. A real monthly dataset was obtained, which contains the prices of rice, wheat, and corn from January 1990 to September 2015. This real dataset makes it possible to compare predictions of food crop prices using single-stage models and integrated models. The remainder of the study is organized as follows: Section2 presents the modeling methodologies of different forecasting schemes. In Section 3, we describe the design structure for our single- and two-stage integrated models. Real data on the prices of three food crops are collected to verify the single-stage models and integrated models. Section 4 contains forecasting comparisons for four single-stage models and three integrated models. Some practical implications and other concerns are addressed in Section 5. The final section summarizes our research findings and contains our conclusions.


## 2. Forecasting Methodologies

### 2.1. Autoregressive Integrated Moving Average

Because seasonal effects are possibly involved in food crop price forecasting, seasonal ARIMA modeling should be developed [22]. A seasonal ARIMA model can be described as follows:
(1)ϕpBΦPBL1−Bd1−BLDZt=δ+θqBΘQBLat,whereB is the backward shift operator, defined as Bjyt=yt−j;L is the number of seasons in a year, and L=12 for monthly data;d is the values of nonseasonal difference transformations;D is the values of seasonal difference transformations;Zt is the working series, which are stationary after fitting a suitable difference transformation from original time series Yt;δ is an unknown constant;at is the white noise at time t, which is independent and identical (iid) with normal distribution;p is the order of nonseasonal autoregressive (AR) models;P is the order of seasonal AR models;q is the order of nonseasonal moving average (MA) models;Q is the order of seasonal MA models;ϕpB is a polynomial function for a nonseasonal AR model, defined as ϕpB=1−ϕ1B−ϕ2B2−ΛΛ−ϕpBp;ΦPBL is a polynomial function for a seasonal AR model, defined as ΦPBL=1−Φ1,LBL−Φ2,LB2L−⋯−ΦP,LBPL;θqB is a polynomial function for a nonseasonal MA model, defined as θqB=1−θ1B−θ2B2−ΛΛ−θqBq;ΘQBL is a polynomial function for a seasonal MA model, defined as ΘQBL=1−Θ1,LBL−Θ2,LB2L−⋯−ΘP,LBQL.The original nonstationary time series should be transformed into a stationary working series through differencing. Typically, the appropriate transformations can be performed by four combinations ofd and D; that is, d,D=0,0, d,D=1,0, d,D=0,1, and d,D=1,1, respectively. Once the stationary working series has been attained, we can observe the behavior of the sample autocorrelation function (ACF) and sample partial autocorrelation function (PACF) to determine the values of p, P, q, and Q for the seasonal ARIMA models.They prefer using the maximum likelihood (ML) technique to obtain estimates for model parameters [22]. For the ML technique, the likelihood function is maximized via nonlinear least squares using Marquardt’s method. Because the ML approach is more computationally expensive than conditional least squares (CLS) estimates, most computer packages employ CLS as a default approach. LS refers to the parameter estimate associated with the smallest sum of squared errors (SSE). For the CLS approach, we assume that the number of past unobserved errors is zero. The data series yt can be represented as follows:
(2)yt=∑i=1∞πiyt−i+at.Theπ weights are calculated as follows:
(3)ϕBθB=1−∑i=1∞πiBi.The CLS approach should produce parameter estimates that can minimize the following:(4)SSE=∑t=1nat∧2=∑t=1nyt−∑i=1∞πi∧yt−12,where the unobserved past values of yt are set to 0 and πi is computed from the estimates of the ARIMA model parameters ϕ and θ at each iteration.After estimating the model parameters, the diagnostic checks for testing the lack of fit of the seasonal ARIMA model should be conducted. A logical way to check the adequacy of a seasonal ARIMA model is to analyze the residuals obtained from the underlying model. The Ljung-Box test was developed to examine whether the firstK sample autocorrelations of residuals indicate adequacy of the model [23]. The null hypothesis for this test is that the first K autocorrelations are jointly zero; that is,
(5)H0:ρ1=ρ2=⋯=ρk=0.The Ljung-Box statistics are described as(6)Q=ndnd+2∑l=1kn′−l−1rl2a^,nd=n−d,whereQ is the Ljung-Box statistics, and the asymptotic distribution of Q follows a chi-square distribution with k−np degrees of freedom;n is the number of observations;d is the degree of nonseasonal differencing;np is the number of parameters other than δ that must be estimated in the ARIMA model under consideration;rl2a^ is the square of the sample autocorrelation of the residuals at lag l.The Ljung-Box test rejects the null hypothesis, which indicates that the underlying model has significant lack of fit if(7)Q>χα2k−np,where χα2k−np is the chi-square distribution table value with k−np degrees of freedom and significance level α.


### 2.2. Artificial Neural Networks

ANN is composed of a large number of highly interconnected processing elements, which are referred to as nodes or neurons, working in parallel to solve a particular problem. While the leftmost layer of the network is referred to as the input layer, the rightmost layer is called the output layer. The middle layer is called the hidden layer. In the ANN structure, each layer consists of a number of nodes connected by links. The ANN contains nodes that are connected to themselves, enabling a node to influence other nodes and itself. ANN includes input variables and output variables. In addition, a certain number of nodes are contained in the hidden layer, and the hidden nodes are nonlinear functions of the input variables. The output variable is a function of the nodes in the hidden layer.ANN modeling is briefly described as follows. For an ANN model, the relationship between inputsI1,I2,…,Ia and output (O) can be represented as
(8)O=α0+∑j=1bαjgδ0j+∑i=1aδijI+ε,where αjj=0,1,2,…,b and δiji=0,1,2,…,a;j=0,1,2,…,b are the model connection weights; a is the number of input nodes; b is the number of hidden nodes; and ε is the error term. In the hidden layer, the transfer function is often used with a logistic function.
(9)gz=11+exp−z.Consequently, the ANN transports nonlinear functional mapping from the inputsI1,I2,…,Ia to the output O; namely,
(10)O=fI1,I2,…,Ia,w+ε,where w is a vector of model parameters and f is a function determined by the ANN structure and connection weights.For the ANN structure, the nodes in the input layers receive input signals from an external source and the nodes in the output layers generate the target output signals. The output of each neuron in the input layer is the same as the input to that neuron. The hidden layers adjust the weights of those inputs until the neural network’s error is minimized. For ANN processing, backpropagation is one method for computing the error contribution of each neuron after a batch of data is processed. This method can be used to adjust the weight of each neuron, thereby completing the learning process for that case. For each neuronj in the hidden layer and each neuron k in the output layer, the net inputs are given by
(11)netj=∑iwji×oi,netk=∑jwkj×oj,where ij is a neuron in the previous layer, oioj is the output of node ij, and wjiwkj is the connection weight from neuron ij to neuron jk. The neuron outputs are given by
(12)oi=neti,oj=11+e−netj+θj=fjnetj,θj,ok=11+e−netk+θk=fknetk,θk,where netjnetk is the input signal from the external source to node jk in the input layer and θjθk is a bias.The generalized delta rule is the conventional technique used to derive the connection weights in the network. First, a set of random numbers is assigned to the connection weights. Then for a presentation of a patternp with a target output vector tp=tp1,tp2,…,tpMT, the sum of squared errors to be minimized is given by
(13)Ep=12∑j=1Mtpj−opj2,where M is the number of output nodes. By minimizing the error Ep using the gradient descent technique, the connection weights can be updated by applying the following equations:
(14)Δwjip=ηδpjopj+αΔwjip−1,where for output nodes,
(15)δpj=tpj−opjopj1−opj,and for hidden nodes,
(16)δpj=opj1−opj∑kδpk×wkj.Note that the learning rate affects the network’s generalization ability and learning speed to a great extent.


### 2.3. Support Vector Regression

SVR is an adaptation of the support vector machine (SVM), one of the most powerful techniques in CI [24, 25]. The basic concept of SVR is to find a model function fx to represent the relationship between the features and target. The modeling of SVR can be described as follows. Suppose
(17)y=fx=w⋅Φx+b,where x and y represent the model input and output vectors, w is the weight vector, b is a constant, Φx denotes a mapping function in the feature space, and w⋅Φx describes the dot production in the feature space F.Typical regression modeling estimates the coefficients by minimizing the square error, which can be considered as an empirical risk based on loss function. Theε-insensitivity loss function was introduced [25] and can be described as follows:
(18)Lεfx,y=x−y−ε,iffx−y≥ε,0,otherwise,where y is the model output and ε is the region of ε-insensitivity. When the predicted value falls into the band area, the loss is zero. However, when the predicted value falls outside the band area, the loss is defined as the difference between the predicted value and margin.When both empirical risk and structure risk are considered, the SVR can be designed to minimize the following quadratic programming problem [25]:
(19)min12w2+C∑i=1nξi+ξi∗,subjecttoyi−w⋅Φxi−b≤ε+ξi,w⋅Φxi+b−yi≤ε+ξi∗,ξi,ξi∗≥0,fori=1,…,n,where i=1,…,n is the number of training observations; ξi,ξi∗ is the empirical risk; 1/2w2 is the structure risk preventing overlearning and lack of applied universality; and C is a modifying coefficient representing the trade-off between empirical and structure risks. With an appropriate modifying coefficient C, band area width ε, and kernel function, the optimum value of each parameter can be solved by the Lagrange procedure. In addition, the general form of the SVR-based regression function can be expressed as follows [25]:
(20)fx,w=fx,r,r∗=∑i=1Nri−ri∗Kx,xi+b,where ri and ri∗ are the Lagrangian multipliers and satisfy the equality riri∗=0. In addition, Kx,xi is the kernel function. The values of the kernel are equal to the inner product of two vectors, xi and xj, in the feature space Φxi and Φxj; that is, Kxi,xj=ΦxiΦxj. Because the radial basis function (RBF) is the most commonly used kernel function [26], we employ it for the experimental study. The RBF is written as follows:
(21)Kxi,xj=exp−xi−xj22σ2,where σ is the width of the RBF.


### 2.4. Multivariate Adaptive Regression Splines

MARS was developed for solving regression-type problems [27]. It is a nonparametric regression procedure that makes no assumption about the functional relationship between the response and explanatory variables. MARS modeling is based on a divide-and-conquer strategy where training datasets are partitioned into separate regions, each of which is assigned to its own regression equation. Consequently, MARS is appropriate and effective for problems with more than two input variables. Particularly, MARS can select important explanatory variables and relationships for complex data structures that often hide in higher dimensional data series.A MARS model can be described as follows [27]:
(22)fx=β0+∑m=1Mβm∏j=1JmSjmxνj,m−ljm,where β0 and βm are the parameters, M is the number of basis functions (BFs), Jm is the number of knots, Sjm takes on values of either 1 or −1 and indicates the right or left sense of the associated step function, νj,m is the label of the independent variable, and ljm is the knot location. The optimal MARS model is determined in a two-step procedure. In the first step, a model is grown by adding BFs until an overly large model is obtained. The BFs are then deleted in the order of least contribution to most contribution by using the generalized cross-validation (GCV) criterion in the second step. The measure of variable importance is provided by observing the decrease in the computed GCV values when a variable is removed from the model. The GCV is described as follows:
(23)GCVM=1/N∑i=1Nyi−fMxi21−CM/N2,where N is the number of observations and CM is the cost penalty measures of a model containing M BF.


## 2.1. Autoregressive Integrated Moving Average

Because seasonal effects are possibly involved in food crop price forecasting, seasonal ARIMA modeling should be developed [22]. A seasonal ARIMA model can be described as follows:
(1)ϕpBΦPBL1−Bd1−BLDZt=δ+θqBΘQBLat,whereB is the backward shift operator, defined as Bjyt=yt−j;L is the number of seasons in a year, and L=12 for monthly data;d is the values of nonseasonal difference transformations;D is the values of seasonal difference transformations;Zt is the working series, which are stationary after fitting a suitable difference transformation from original time series Yt;δ is an unknown constant;at is the white noise at time t, which is independent and identical (iid) with normal distribution;p is the order of nonseasonal autoregressive (AR) models;P is the order of seasonal AR models;q is the order of nonseasonal moving average (MA) models;Q is the order of seasonal MA models;ϕpB is a polynomial function for a nonseasonal AR model, defined as ϕpB=1−ϕ1B−ϕ2B2−ΛΛ−ϕpBp;ΦPBL is a polynomial function for a seasonal AR model, defined as ΦPBL=1−Φ1,LBL−Φ2,LB2L−⋯−ΦP,LBPL;θqB is a polynomial function for a nonseasonal MA model, defined as θqB=1−θ1B−θ2B2−ΛΛ−θqBq;ΘQBL is a polynomial function for a seasonal MA model, defined as ΘQBL=1−Θ1,LBL−Θ2,LB2L−⋯−ΘP,LBQL.The original nonstationary time series should be transformed into a stationary working series through differencing. Typically, the appropriate transformations can be performed by four combinations ofd and D; that is, d,D=0,0, d,D=1,0, d,D=0,1, and d,D=1,1, respectively. Once the stationary working series has been attained, we can observe the behavior of the sample autocorrelation function (ACF) and sample partial autocorrelation function (PACF) to determine the values of p, P, q, and Q for the seasonal ARIMA models.They prefer using the maximum likelihood (ML) technique to obtain estimates for model parameters [22]. For the ML technique, the likelihood function is maximized via nonlinear least squares using Marquardt’s method. Because the ML approach is more computationally expensive than conditional least squares (CLS) estimates, most computer packages employ CLS as a default approach. LS refers to the parameter estimate associated with the smallest sum of squared errors (SSE). For the CLS approach, we assume that the number of past unobserved errors is zero. The data series yt can be represented as follows:
(2)yt=∑i=1∞πiyt−i+at.Theπ weights are calculated as follows:
(3)ϕBθB=1−∑i=1∞πiBi.The CLS approach should produce parameter estimates that can minimize the following:(4)SSE=∑t=1nat∧2=∑t=1nyt−∑i=1∞πi∧yt−12,where the unobserved past values of yt are set to 0 and πi is computed from the estimates of the ARIMA model parameters ϕ and θ at each iteration.After estimating the model parameters, the diagnostic checks for testing the lack of fit of the seasonal ARIMA model should be conducted. A logical way to check the adequacy of a seasonal ARIMA model is to analyze the residuals obtained from the underlying model. The Ljung-Box test was developed to examine whether the firstK sample autocorrelations of residuals indicate adequacy of the model [23]. The null hypothesis for this test is that the first K autocorrelations are jointly zero; that is,
(5)H0:ρ1=ρ2=⋯=ρk=0.The Ljung-Box statistics are described as(6)Q=ndnd+2∑l=1kn′−l−1rl2a^,nd=n−d,whereQ is the Ljung-Box statistics, and the asymptotic distribution of Q follows a chi-square distribution with k−np degrees of freedom;n is the number of observations;d is the degree of nonseasonal differencing;np is the number of parameters other than δ that must be estimated in the ARIMA model under consideration;rl2a^ is the square of the sample autocorrelation of the residuals at lag l.The Ljung-Box test rejects the null hypothesis, which indicates that the underlying model has significant lack of fit if(7)Q>χα2k−np,where χα2k−np is the chi-square distribution table value with k−np degrees of freedom and significance level α.


## 2.2. Artificial Neural Networks

ANN is composed of a large number of highly interconnected processing elements, which are referred to as nodes or neurons, working in parallel to solve a particular problem. While the leftmost layer of the network is referred to as the input layer, the rightmost layer is called the output layer. The middle layer is called the hidden layer. In the ANN structure, each layer consists of a number of nodes connected by links. The ANN contains nodes that are connected to themselves, enabling a node to influence other nodes and itself. ANN includes input variables and output variables. In addition, a certain number of nodes are contained in the hidden layer, and the hidden nodes are nonlinear functions of the input variables. The output variable is a function of the nodes in the hidden layer.ANN modeling is briefly described as follows. For an ANN model, the relationship between inputsI1,I2,…,Ia and output (O) can be represented as
(8)O=α0+∑j=1bαjgδ0j+∑i=1aδijI+ε,where αjj=0,1,2,…,b and δiji=0,1,2,…,a;j=0,1,2,…,b are the model connection weights; a is the number of input nodes; b is the number of hidden nodes; and ε is the error term. In the hidden layer, the transfer function is often used with a logistic function.
(9)gz=11+exp−z.Consequently, the ANN transports nonlinear functional mapping from the inputsI1,I2,…,Ia to the output O; namely,
(10)O=fI1,I2,…,Ia,w+ε,where w is a vector of model parameters and f is a function determined by the ANN structure and connection weights.For the ANN structure, the nodes in the input layers receive input signals from an external source and the nodes in the output layers generate the target output signals. The output of each neuron in the input layer is the same as the input to that neuron. The hidden layers adjust the weights of those inputs until the neural network’s error is minimized. For ANN processing, backpropagation is one method for computing the error contribution of each neuron after a batch of data is processed. This method can be used to adjust the weight of each neuron, thereby completing the learning process for that case. For each neuronj in the hidden layer and each neuron k in the output layer, the net inputs are given by
(11)netj=∑iwji×oi,netk=∑jwkj×oj,where ij is a neuron in the previous layer, oioj is the output of node ij, and wjiwkj is the connection weight from neuron ij to neuron jk. The neuron outputs are given by
(12)oi=neti,oj=11+e−netj+θj=fjnetj,θj,ok=11+e−netk+θk=fknetk,θk,where netjnetk is the input signal from the external source to node jk in the input layer and θjθk is a bias.The generalized delta rule is the conventional technique used to derive the connection weights in the network. First, a set of random numbers is assigned to the connection weights. Then for a presentation of a patternp with a target output vector tp=tp1,tp2,…,tpMT, the sum of squared errors to be minimized is given by
(13)Ep=12∑j=1Mtpj−opj2,where M is the number of output nodes. By minimizing the error Ep using the gradient descent technique, the connection weights can be updated by applying the following equations:
(14)Δwjip=ηδpjopj+αΔwjip−1,where for output nodes,
(15)δpj=tpj−opjopj1−opj,and for hidden nodes,
(16)δpj=opj1−opj∑kδpk×wkj.Note that the learning rate affects the network’s generalization ability and learning speed to a great extent.


## 2.3. Support Vector Regression

SVR is an adaptation of the support vector machine (SVM), one of the most powerful techniques in CI [24, 25]. The basic concept of SVR is to find a model function fx to represent the relationship between the features and target. The modeling of SVR can be described as follows. Suppose
(17)y=fx=w⋅Φx+b,where x and y represent the model input and output vectors, w is the weight vector, b is a constant, Φx denotes a mapping function in the feature space, and w⋅Φx describes the dot production in the feature space F.Typical regression modeling estimates the coefficients by minimizing the square error, which can be considered as an empirical risk based on loss function. Theε-insensitivity loss function was introduced [25] and can be described as follows:
(18)Lεfx,y=x−y−ε,iffx−y≥ε,0,otherwise,where y is the model output and ε is the region of ε-insensitivity. When the predicted value falls into the band area, the loss is zero. However, when the predicted value falls outside the band area, the loss is defined as the difference between the predicted value and margin.When both empirical risk and structure risk are considered, the SVR can be designed to minimize the following quadratic programming problem [25]:
(19)min12w2+C∑i=1nξi+ξi∗,subjecttoyi−w⋅Φxi−b≤ε+ξi,w⋅Φxi+b−yi≤ε+ξi∗,ξi,ξi∗≥0,fori=1,…,n,where i=1,…,n is the number of training observations; ξi,ξi∗ is the empirical risk; 1/2w2 is the structure risk preventing overlearning and lack of applied universality; and C is a modifying coefficient representing the trade-off between empirical and structure risks. With an appropriate modifying coefficient C, band area width ε, and kernel function, the optimum value of each parameter can be solved by the Lagrange procedure. In addition, the general form of the SVR-based regression function can be expressed as follows [25]:
(20)fx,w=fx,r,r∗=∑i=1Nri−ri∗Kx,xi+b,where ri and ri∗ are the Lagrangian multipliers and satisfy the equality riri∗=0. In addition, Kx,xi is the kernel function. The values of the kernel are equal to the inner product of two vectors, xi and xj, in the feature space Φxi and Φxj; that is, Kxi,xj=ΦxiΦxj. Because the radial basis function (RBF) is the most commonly used kernel function [26], we employ it for the experimental study. The RBF is written as follows:
(21)Kxi,xj=exp−xi−xj22σ2,where σ is the width of the RBF.


## 2.4. Multivariate Adaptive Regression Splines

MARS was developed for solving regression-type problems [27]. It is a nonparametric regression procedure that makes no assumption about the functional relationship between the response and explanatory variables. MARS modeling is based on a divide-and-conquer strategy where training datasets are partitioned into separate regions, each of which is assigned to its own regression equation. Consequently, MARS is appropriate and effective for problems with more than two input variables. Particularly, MARS can select important explanatory variables and relationships for complex data structures that often hide in higher dimensional data series.A MARS model can be described as follows [27]:
(22)fx=β0+∑m=1Mβm∏j=1JmSjmxνj,m−ljm,where β0 and βm are the parameters, M is the number of basis functions (BFs), Jm is the number of knots, Sjm takes on values of either 1 or −1 and indicates the right or left sense of the associated step function, νj,m is the label of the independent variable, and ljm is the knot location. The optimal MARS model is determined in a two-step procedure. In the first step, a model is grown by adding BFs until an overly large model is obtained. The BFs are then deleted in the order of least contribution to most contribution by using the generalized cross-validation (GCV) criterion in the second step. The measure of variable importance is provided by observing the decrease in the computed GCV values when a variable is removed from the model. The GCV is described as follows:
(23)GCVM=1/N∑i=1Nyi−fMxi21−CM/N2,where N is the number of observations and CM is the cost penalty measures of a model containing M BF.


## 3. Food Crop Price Forecasting

### 3.1. Datasets and Forecasting Criteria

In this study, the price data of the three most important food crops were collected for the period of January 1990 to September 2015 from the web sites of IndexMundi [28–31]. The datasets consist of 309 records for each food crop price. Among them, the first 297 data records were used to develop the different forecasting models, while the remaining 12 data records were used to perform model validation. This study has presented an experiment based on a larger population of data and has foreseen obtaining accurate forecasts of food crop price, with the forecast horizon for the out-of-sample forecasting experiment of 12. Good prediction of prices may assist in planning agricultural supply, but the investigated data does not contain any information on population.The forecasting capability of the models was compared using three forecasting accuracy criteria, including mean absolute percentage error (MAPE), root mean square error (RMSE), and mean absolute error (MAE). These forecasting measurements are expressed as follows:
(24)MAPE=1n∑t=1netYt×100%,MAE=1n∑t=1net,RMSE=1n∑t=1net2,where et is the value of the residual at time t and Yt is the observation at time t.Obviously, it can be noted that the lower theMAPE, MAE, and RMSE values, the closer the forecasted values to the actual values.


### 3.2. Forecast Modeling of Rice Prices

Figure1 shows the original time plot for the rice price data series. This study uses an SAS package to run the ARIMA modeling. Table 1 shows the parameter estimates, and all the estimated parameters are significantly different from zero. In Table 1, the notations “AR1,1,” “AR1,2,” “AR1,3,” and “AR1,4” correspond to the parameters ϕ1, ϕ2, ϕ7, and ϕ13 of the AR model. The ARIMA model presented in Table 1 was a subset model. Since the parameters, ϕ1, ϕ2, ϕ7, and ϕ13, are significantly different from zero, this study refers to the model in Table 1 as an AR (1, 2, 7, and 13) model.Figure 1
Time plot for rice price (unit: US dollars per metric ton).Table 1
Parameter estimates for rice prices.
The ARIMA procedureConditional least squares estimation
Parameter
Estimate
Standard error
t value
Approx Pr > |t|
Lag
AR1,1
0.46514
0.05586
8.33
<0.0001
1
AR1,2
−0.23278
0.05559
−4.19
<0.0001
2
AR1,3
−0.17055
0.05183
−3.29
0.0011
7
AR1,4
−0.11264
0.05170
−2.18
0.0301
13In addition, Table2 shows the Ljung-Box test results. As mentioned earlier, the Ljung-Box statistic can be used to check the fit of an ARIMA model. A simple method to derive the hypothesis testing result is described as follows. We should reject the model under consideration by setting the type I error equal to α if and only if the p value is less than α. Generally, α is set equal to 0.05. The SAS package computes the Ljung-Box statistics Q and their associated p values for values of K equal to 6, 12, 18, 24, 30, and 36. In Table 2, the first column contains the values of K, the second column contains the values of Q, and the fourth column contains the associated p values. By studying Table 2, we observe that all the associated p values are greater than α. Accordingly, a conclusion can be made concerning the fit of the underlying ARIMA model. Thus, the model expressed in (25) is suitable for modeling rice prices.
(25)Zt=0.46514Zt−1−0.502Zt−2−0.17055Zt−7−0.11264Zt−13+at,where
(26)Zt=yt−yt−1.Table 2
Ljung-Box test results for ARIMA modeling of rice prices.
Autocorrelation check of residuals
To Lag
Chi-square
DF
Pr > ChiSq
Autocorrelations
6
4.61
2
0.0995
−0.002
0.007
−0.024
−0.073
0.006
−0.096
12
6.34
8
0.6097
0.009
0.007
0.016
0.027
−0.031
0.059
18
13.10
14
0.5190
−0.023
0.028
0.028
−0.09
0.061
−0.087
24
14.40
20
0.8096
0.059
0.012
0.008
−0.006
−0.012
−0.015
30
19.29
26
0.8242
−0.102
−0.045
−0.025
−0.002
0.026
0.036
36
20.10
32
0.9494
0.018
−0.022
−0.002
−0.024
−0.003
0.032For ANN designs, there is no fixed mode to decide the number of hidden nodes. Too few hidden nodes confine the network generalization capability, whereas too many hidden nodes may lead to overtraining difficulties. Therefore, in this study, we consider the hidden nodes to be set from (2n−2) to (2n+2) if n≦8, where n represents the number of input variables. When n>8, we consider the hidden nodes to be set from (n−2) to (n+2). In this study, we denote the term {ni-nh-no} for ANN parameter settings, where ni is the number of neurons in the input layer, nh stands for the number of neurons in the hidden layer, and no represents the number of neurons in the output layer, respectively. Additionally, this study employs the learning rate for all ANN models at the default value (i.e., 0.01) to ensure consistency [32]. The network topology with the minimum MAPE is considered as the optimal network topology, due to the fact that the MAPE is one of the most important performance measurements of forecasting capability.In this study, because there are no explanatory variables available, the inputs of ANN modeling for food crop price forecasting becomes unfeasible. Consequently, this study uses self-predictor variables to serve as the inputs for ANN. In here, the self-predictor variable is defined as the lags or past values of a variable. This study proposes two input designs for the modeling of CI techniques. Because seasonal effects may influence food crop price forecasts, the first design rationally selects the preceding 12 observations to forecast food crop prices at timet. Accordingly, the first design used 12 self-predictor variables (i.e., yt−1,yt−2,yt−3,…, and yt−12) to serve as input variables and employed a single variable (i.e., yt) to serve as the output for ANN modeling. That is, the first design considered 12 input nodes and one output node for the ANN structures. This design is denoted as ANN1.The second input design for CI technique modeling is performed using FSM. For this proposed modeling, this study incorporated FSMs with ANN, SVR, and MARS to develop forecasting models. The proposed input design used ARIMA as an FSM to extract important self-predictor variables, which serve as the inputs for the ANN, SVR, and MARS models. Because the ARIMA model for rice price forecasting is obtained in (25), it indicates that Zt−1,Zt−2,Zt−7, and Zt−13 will influence Zt. Consequently, this study selected Zt−1,Zt−2,Zt−7, and Zt−13 to serve as self-predictor, or input, variables for the proposed integrated modeling on the rice price data series. In this study, this second input design is denoted as ARIMA-ANN.To forecast the rice prices, the two ANN models contain twelve and four input nodes for the first and second designs, respectively. The hidden nodes were selected as 10, 11, 12, 13, and 14 for the first design and 6, 7, 8, 9, and 10 for the second design. After performing the ANN modeling, the ANN1 design showed that the {12-14-1} structure provided the best results and minimum testing MAPE for rice prices. For the ARIMA-ANN design, the best structure was {4-7-1}. Three forecasting measurements for different settings of the ANN topologies for the two designs are shown in Table 3.Table 3
Forecasting measurements for two ANN designs (rice price forecasting).
ANN1 topology
MAE
RMSE
MAPE (%)
ARIMA-ANN topology
MAE
RMSE
MAPE (%)
{12-10-1}
15.438
16.648
3.939
{4-6-1}
10.384
12.053
2.654
{12-11-1}
14.737
16.022
3.765
{4-7-1}
10.365
12.037
2.650
{12-12-1}
15.001
16.262
3.825
{4-8-1}
10.374
12.044
2.652
{12-13-1}
15.182
16.499
3.869
{4-9-1}
10.379
12.062
2.653
{12-14-1}
13.915
15.320
3.557
{4-10-1}
10.389
12.039
2.655For modeling SVR on rice prices, this study adopted the input structure as ANN modeling. SVR modeling employed two designs for the input variables. The first design used 12 variables,yt−1,yt−2,yt−3,…, and yt−12 as the input variables and yt as the output variable. The second design used four self-predictor variables, Zt−1,Zt−2,Zt−7, and Zt−13 to serve as the input variables and Zt as the output variable. The first and second SVR designs are denoted as SVR1 and ARIMA-SVR, respectively. Because the parameter settings of C and ε often affect the performance of SVR modeling, an analytic parameter selection approach and grid search are used in this study [26]. Accordingly, we denote the term {ε, C}SVR as the best SVR parameter settings. After performing the SVR modeling, the SVR1 design reported that the parameter settings of {2−4, 211}SVR provided MAPE = 6.912% for rice price forecasting. For the second design, ARIMA-SVR, the parameter settings of {2−5, 29}SVR provided MAPE = 2.057%.For MARS modeling, this study also adopted the same input structure of ANN or SVR modeling. The first and second MARS designs are denoted as MARS1 and ARIMA-MARS, respectively. Table 4 lists the variable selection results and BFs after modeling rice price data using MARS1. In Table 4, the first column contains the variables that should be included in the model, the second column contains the relative importance of the variables that are listed in the first column, the third column contains the various BFs, and the final column contains the estimated coefficient values for the BFs that are listed in the third column. Based on the results in Table 4, it can be seen that three input variables (i.e., yt−12,yt−11, and yt−10) played important roles in building the MARS forecasting models. Their relative importance values (%) are 100, 13.1, and 12.1. The construction of the BFs can be expressed as follows. Consider the variable yt−12 as an example. We have two BFs (i.e., BF1 and BF2) to be considered. The observed values of BF1 are determined by the values larger than 0 or yt−12−367. The corresponding coefficient is estimated to be 0.985. The observed values of BF2 are determined by the values larger than 0 or 552.09−yt−12. The corresponding coefficient is estimated to be −1.258. Accordingly, the MARS forecasting model for rice prices can be expressed as follows:
(27)Y=542.484+0.985BF1−1.258BF2+1.539BF3−1.880BF4+0.524BF5−2.410BF6+1.783BF7+0.454BF8−0.234BF9+0.346BF10.Table 4
The results of MARS1 for rice price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
542.484
yt−12
100.0
BF1 = Max (0, yt−12 − 367)
0.985
BF2 = Max (0, 552.09 − yt−12)
−1.258
yt−11
13.1
BF3 = Max (0, yt−11 − 356)
1.539
BF4 = Max (0, yt−11 − 472.48)
−1.880
BF5 = Max (0, 585.95 − yt−11)
0.524
yt−10
12.1
BF6 = Max (0, yt−10 − 363)
−2.410
BF7 = Max (0, yt−10 − 472.48)
1.783
BF8 = Max (0, yt−10 − 543.14)
0.454
BF9 = Max (0, 623 − yt−10)
−0.234
BF10 = Max (0, yt−10 − 623)
0.346For ARIMA-MARS modeling, the variable selection results and BFs are summarized in Table5. In addition, it can be seen that three input variables (i.e., Zt−13,Zt−7, and Zt−1) played important roles in building the ARIMA-MARS forecasting models. The MARS forecasting model for rice prices can be described as follows:
(28)Y=−13.404+0.407BF1+23.543BF2−73.468BF3+49.716BF4+0.479BF5,where
(29)Zt=yt−yt−1.Table 5
The results of ARIMA-MARS for rice price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
−13.404
Zt−13
100.0
BF1 = Max (0, Zt−13 + 35)
0.407
Zt−7
45.4
BF2 = Max (0, Zt−7 − 0)
23.543
BF3 = Max (0, Zt−7 + 1.30)
−73.468
BF4 = Max (0, Zt−7 + 1.87)
49.716
Zt−1
28.5
BF5 = Max (0, −39.85 − Zt−1)
0.479


### 3.3. Forecast Modeling of Wheat Prices

Figure2 displays the original time plot for wheat prices. After modeling the wheat price data series with the ARIMA procedure, the parameter estimates are shown in Table 6. In Table 6, the notations “MA1,1,” “MA1,2,” “MA1,3,” and “MA1,4” correspond to the parameters θ1, θ7, θ12, and θ13 of the MA model. Table 7 displays the Ljung-Box test results, and a conclusion is made on the appropriateness of the underlying ARIMA model. Accordingly, (30) is suitable for modeling wheat prices.
(30)Zt=at−0.2402at−1+0.14629at−7+0.11171at−12+0.13721at−13,where
(31)Zt=yt−yt−1.Figure 2
Time plot for wheat price (unit: US dollars per metric ton).Table 6
Parameter estimates for wheat prices.
The ARIMA procedureConditional least squares estimation
Parameter
Estimate
Standard error
t value
Approx Pr > |t|
Lag
MA1,1
−0.2402
0.05601
−4.29
<0.0001
1
MA1,2
0.14629
0.05754
2.54
0.0115
7
MA1,3
0.11171
0.0599
1.86
0.0632
12
MA1,4
0.13721
0.06018
2.28
0.0233
13Table 7
Ljung-Box test results for ARIMA modeling of wheat prices.
Autocorrelation check of residuals
To Lag
Chi-square
DF
Pr > ChiSq
Autocorrelations
6
4.4
2
0.111
0.019
−0.045
0.024
−0.036
0.069
0.074
12
11.07
8
0.1977
−0.03
−0.126
0.008
−0.055
0.042
−0.001
18
14.41
14
0.4197
−0.035
−0.084
0.021
−0.038
0.021
0.006
24
19.99
20
0.4589
0.067
0.081
−0.033
−0.03
0.012
−0.064
30
30.52
26
0.2465
−0.007
−0.041
−0.11
−0.132
−0.009
0.026
36
33.96
32
0.3733
−0.025
0.025
−0.071
0.011
0.054
0.029In addition,at−j stands for the white noise at time t−j.For ANN modeling on wheat prices, this study also used two ANN designs, same as the structures for modeling rice prices. The first design used 12 self-predictor variables (i.e.,yt−1,yt−2,yt−3,…, and yt−12) to serve as input variables and a single variable (i.e., yt) to serve as the output variable. The second design considers the ARIMA as an FSM to extract important self-predictor variables, which serve as the inputs for the ANN models. The ARIMA model for wheat price forecasting is obtained in (30) and indicates that at−1,at−7,at−12, and at−13 will influence Zt.Becauseat−1=Zt−1−Zt−1∧, where Zt−1∧ is the forecast at time t−1, we observe that Zt−1 will influence Zt. Therefore, Zt−1 should be selected as an input variable as long as the component at−1 is in the MA model. Following the same logic, because the model (e.i., (30)) contains at−1,at−7,at−12, and at−13, we conclude that Zt−1,Zt−7,Zt−12, and Zt−13 will influence Zt. Accordingly, we have selected Zt−1,Zt−7,Zt−12, and Zt−13 to serve as input variables for the second proposed ANN modeling design for the wheat price data series.After performing ANN modeling on wheat price data, Table8 shows the corresponding forecast validity measures for different settings of ANN topologies for the two designs. As shown in Table 8, we observe that the ANN1 design with the {12-11-1} structure has the smallest MAPE. For the second, or ARIMA-ANN design, the {4-8-1} structure was associated with the smallest MAPE.Table 8
Forecast validity measures for two ANN designs (rice price forecasting).
ANN1 topology
MAE
RMSE
MAPE (%)
ARIMA-ANN topology
MAE
RMSE
MAPE (%)
{12-10-1}
12.121
13.521
5.448
{4-6-1}
11.833
14.336
5.374
{12-11-1}
11.752
13.230
5.285
{4-7-1}
11.865
14.353
5.388
{12-12-1}
12.084
13.515
5.423
{4-8-1}
8.984
10.864
4.077
{12-13-1}
11.984
13.429
5.385
{4-9-1}
8.984
10.883
4.078
{12-14-1}
12.120
13.516
5.457
{4-10-1}
11.845
14.346
5.379For modeling SVR on wheat prices, this study adopted the input structure as ANN modeling. The first design employed 12 variables,yt−1,yt−2,yt−3,…, and yt−12 as the input variables and yt as the output variable. The second design employed four self-predictor variables, Zt−1,Zt−7,Zt−12, and Zt−13 to serve as the input variables and Zt as the output variable. After performing SVR modeling on the wheat price data series, we obtained the parameter settings of {2−8, 210}SVR and associated MAPE = 4.661% for the first design. In addition, we had parameter settings of {2−7, 2−1}SVR and associated MAPE = 4.347% for the second design.For MARS modeling with the first design, Table9 presents the variable selection results and BFs. As shown in Table 9, we notice that three input variables (i.e., yt−12,yt−6, and yt−11) played important roles in building the MARS forecasting models. In addition, the MARS forecasting model for wheat prices can be expressed as follows:
(32)Y=360.973+0.765BF1−1.272BF2−0.185BF3+0.298BF4−0.218BF5.Table 9
The results of MARS1 for wheat price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
360.973
yt−12
100.0
BF1 = Max (0, yt−12 − 321.18)
0.765
BF2 = Max (0, 321.81 − yt−12)
−1.272
yt−6
7.2
BF3 = Max (0, yt−6 − 178.32)
−0.185
BF4 = Max (0, 178.32 − yt−6)
0.298
yt−11
3.2
BF5 = Max (0, yt−11 − 269.03)
−0.218For MARS modeling with the second design, the variable selection results and BFs are summarized in Table10. Furthermore, the MARS forecasting model for wheat prices is expressed as follows:
(33)Y=−13.665−1.105BF1+8.333BF2−13.064BF3+4.955BF4−0.155BF5.Table 10
The results of ARIMA-MARS for wheat price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
−13.665
Zt−7
100.0
BF1 = Max (0, Zt−7 − 26.58)
−1.105
Zt−13
63.3
BF2 = Max (0, Zt−13 + 26.91)
8.333
BF3 = Max (0, Zt−13 + 22.55)
−13.064
BF4 = Max (0, Zt−13 + 17.56)
4.955
Zt−1
14.4
BF5 = Max (0, Zt−1 + 10.65)
−0.155


### 3.4. Forecast Modeling of Corn Prices

Figure3 presents the original time plot for corn prices. Table 11 presents the parameter estimates after performing ARIMA modeling. In Table 11, the notations “AR1,1,” “AR1,2,” “AR1,3,” “AR1,4,” “AR1,5,” and “AR1,6” correspond to as parameters ϕ1, ϕ3, ϕ6, ϕ12, ϕ14, and ϕ24 of the AR model. In Table 11, we still contain the parameters of “AR1,2” and “AR1,4” in the model, although the absolute t values are less than or equal to 2 (i.e., the typical type I error is chosen as 0.05). The main reason is that those two parameters are important, and the Ljung-Box statistics indicate that the model is not appropriate if those two parameters are not involved in the underlying ARIMA model. Additionally, Table 12 demonstrates the Ljung-Box test results and a conclusion is made on the appropriateness of the ARIMA model. Thus, (34) is a suitable ARIMA model for modeling the corn price data series.
(34)Zt=0.17626Zt−1+0.09344Zt−3−0.145Zt−6−0.0925Zt−12−0.1476Zt−14−0.2149Zt−24+at,where
(35)Zt=yt−yt−1.Figure 3
Time plot for corn price (unit: US dollars per metric ton).Table 11
Parameter estimates for corn prices.
The ARIMA procedureConditional least squares estimation
Parameter
Estimate
Standard error
t value
Approx Pr > |t|
Lag
AR1,1
0.17626
0.05909
2.98
0.0031
1
AR1,2
0.09344
0.0651
1.44
0.1524
3
AR1,3
−0.145
0.06527
−2.22
0.0272
6
AR1,4
−0.0925
0.06619
−1.4
0.1636
12
AR1,5
−0.1476
0.06719
−2.2
0.029
14
AR1,6
−0.2149
0.07024
−3.06
0.0024
24Table 12
Ljung-Box test results for ARIMA modeling of corn prices.
Autocorrelation check of residuals
To Lag
Chi-square
DF
Pr > ChiSq
Autocorrelations
6
.
0
.
−0.011
0.082
0.026
0.002
−0.063
−0.027
12
6.71
6
0.3481
0.013
−0.032
0.001
0.012
0.091
−0.014
18
14.3
12
0.2822
−0.069
−0.036
0.042
0.113
0.038
−0.046
24
16.43
18
0.5624
0.014
−0.033
−0.015
0.061
−0.034
−0.015
30
29.11
24
0.2159
0.035
−0.073
0.103
−0.079
−0.025
−0.12
36
30.86
30
0.4225
−0.036
0
0.06
0.005
0.015
−0.001For ANN modeling on corn prices, while the first design used 12 input variables, the second design used six input variables extracted by using the FSM of ARIMA modeling. After performing ANN modeling, Table13 presents the corresponding forecast validity measures for different ANN topology settings for the two designs. From Table 13, we observe that the ANN1 design with the {12-14-1} structure has the smallest MAPE. For the second design, the {6-13-1} structure has the smallest MAPE.Table 13
Forecast validity measures for two ANN designs (corn price forecasting).
ANN1 topology
MAE
RMSE
MAPE (%)
ARIMA-ANN topology
MAE
RMSE
MAPE (%)
{12-10-1}
7.536
9.698
4.387
{6-10-1}
7.304
10.773
4.267
{12-11-1}
7.464
9.635
4.343
{6-11-1}
7.307
10.853
4.269
{12-12-1}
7.444
9.638
4.331
{6-12-1}
7.291
10.805
4.259
{12-13-1}
7.487
9.633
4.357
{6-13-1}
7.286
10.771
4.256
{12-14-1}
7.433
9.589
4.326
{6-14-1}
7.290
10.809
4.259In the modeling of SVR for the corn price data series, this study adopted the same input structure as ANN modeling. The first design employed 12 variables (i.e.,yt−1,yt−2,yt−3,…, and yt−12) as the input variables and yt as the output variable. The second design employed six self-predictor variables (i.e., Zt−1,Zt−3,Zt−6,Zt−12,Zt−14, and Zt−24), which were selected by using the FSM of ARIMA modeling to serve as the input variables and Zt as the output variable. After performing SVR modeling, we obtained the parameter settings of {2−9, 27}SVR and associated MAPE = 4.060% for the first design. In addition, we obtained parameter settings of {2−11, 211}SVR and associated MAPE = 4.215% for the second design.For MARS modeling on corn price data with the first design, Table14 illustrates the variable selection results and BFs. As shown in Table 14, we notice that two variables (i.e., yt−12 and yt−11) play important roles in building the MARS forecasting models. Consequently, the MARS forecasting model for corn prices can be described as follows:
(36)Y=262.148−1.354BF1+0.969BF2+0.342BF3.Table 14
The results of MARS1 for corn price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
262.148
yt−12
100.0
BF1 = Max (0, 268.79 − yt−12)
−1.354
BF2 = Max (0, yt−12 − 268.79)
0.969
yt−11
7.2
BF3 = Max (0, 295.29 − yt−11)
0.342For MARS modeling with the second design, Table15 lists the variable selection results and BFs. The MARS forecasting model for the corn prices is expressed as follows:
(37)Y=3.787−0.277BF1−0.745BF2−0.466BF3+0.175BF4.Table 15
The results of ARIMA-MARS for corn price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
3.787
Zt−24
100.0
BF1 = Max (0, 15.08 − Zt−24)
−0.277
Zt−12
74.0
BF2 = Max (0, Zt−12 − 13.9)
−0.745
Zt−1
44.8
BF3 = Max (0, Zt−1 − 11.59)
−0.466
Zt−6
17.7
BF4 = Max (0, Zt−6 + 7.36)
0.175


## 3.1. Datasets and Forecasting Criteria

In this study, the price data of the three most important food crops were collected for the period of January 1990 to September 2015 from the web sites of IndexMundi [28–31]. The datasets consist of 309 records for each food crop price. Among them, the first 297 data records were used to develop the different forecasting models, while the remaining 12 data records were used to perform model validation. This study has presented an experiment based on a larger population of data and has foreseen obtaining accurate forecasts of food crop price, with the forecast horizon for the out-of-sample forecasting experiment of 12. Good prediction of prices may assist in planning agricultural supply, but the investigated data does not contain any information on population.The forecasting capability of the models was compared using three forecasting accuracy criteria, including mean absolute percentage error (MAPE), root mean square error (RMSE), and mean absolute error (MAE). These forecasting measurements are expressed as follows:
(24)MAPE=1n∑t=1netYt×100%,MAE=1n∑t=1net,RMSE=1n∑t=1net2,where et is the value of the residual at time t and Yt is the observation at time t.Obviously, it can be noted that the lower theMAPE, MAE, and RMSE values, the closer the forecasted values to the actual values.


## 3.2. Forecast Modeling of Rice Prices

Figure1 shows the original time plot for the rice price data series. This study uses an SAS package to run the ARIMA modeling. Table 1 shows the parameter estimates, and all the estimated parameters are significantly different from zero. In Table 1, the notations “AR1,1,” “AR1,2,” “AR1,3,” and “AR1,4” correspond to the parameters ϕ1, ϕ2, ϕ7, and ϕ13 of the AR model. The ARIMA model presented in Table 1 was a subset model. Since the parameters, ϕ1, ϕ2, ϕ7, and ϕ13, are significantly different from zero, this study refers to the model in Table 1 as an AR (1, 2, 7, and 13) model.Figure 1
Time plot for rice price (unit: US dollars per metric ton).Table 1
Parameter estimates for rice prices.
The ARIMA procedureConditional least squares estimation
Parameter
Estimate
Standard error
t value
Approx Pr > |t|
Lag
AR1,1
0.46514
0.05586
8.33
<0.0001
1
AR1,2
−0.23278
0.05559
−4.19
<0.0001
2
AR1,3
−0.17055
0.05183
−3.29
0.0011
7
AR1,4
−0.11264
0.05170
−2.18
0.0301
13In addition, Table2 shows the Ljung-Box test results. As mentioned earlier, the Ljung-Box statistic can be used to check the fit of an ARIMA model. A simple method to derive the hypothesis testing result is described as follows. We should reject the model under consideration by setting the type I error equal to α if and only if the p value is less than α. Generally, α is set equal to 0.05. The SAS package computes the Ljung-Box statistics Q and their associated p values for values of K equal to 6, 12, 18, 24, 30, and 36. In Table 2, the first column contains the values of K, the second column contains the values of Q, and the fourth column contains the associated p values. By studying Table 2, we observe that all the associated p values are greater than α. Accordingly, a conclusion can be made concerning the fit of the underlying ARIMA model. Thus, the model expressed in (25) is suitable for modeling rice prices.
(25)Zt=0.46514Zt−1−0.502Zt−2−0.17055Zt−7−0.11264Zt−13+at,where
(26)Zt=yt−yt−1.Table 2
Ljung-Box test results for ARIMA modeling of rice prices.
Autocorrelation check of residuals
To Lag
Chi-square
DF
Pr > ChiSq
Autocorrelations
6
4.61
2
0.0995
−0.002
0.007
−0.024
−0.073
0.006
−0.096
12
6.34
8
0.6097
0.009
0.007
0.016
0.027
−0.031
0.059
18
13.10
14
0.5190
−0.023
0.028
0.028
−0.09
0.061
−0.087
24
14.40
20
0.8096
0.059
0.012
0.008
−0.006
−0.012
−0.015
30
19.29
26
0.8242
−0.102
−0.045
−0.025
−0.002
0.026
0.036
36
20.10
32
0.9494
0.018
−0.022
−0.002
−0.024
−0.003
0.032For ANN designs, there is no fixed mode to decide the number of hidden nodes. Too few hidden nodes confine the network generalization capability, whereas too many hidden nodes may lead to overtraining difficulties. Therefore, in this study, we consider the hidden nodes to be set from (2n−2) to (2n+2) if n≦8, where n represents the number of input variables. When n>8, we consider the hidden nodes to be set from (n−2) to (n+2). In this study, we denote the term {ni-nh-no} for ANN parameter settings, where ni is the number of neurons in the input layer, nh stands for the number of neurons in the hidden layer, and no represents the number of neurons in the output layer, respectively. Additionally, this study employs the learning rate for all ANN models at the default value (i.e., 0.01) to ensure consistency [32]. The network topology with the minimum MAPE is considered as the optimal network topology, due to the fact that the MAPE is one of the most important performance measurements of forecasting capability.In this study, because there are no explanatory variables available, the inputs of ANN modeling for food crop price forecasting becomes unfeasible. Consequently, this study uses self-predictor variables to serve as the inputs for ANN. In here, the self-predictor variable is defined as the lags or past values of a variable. This study proposes two input designs for the modeling of CI techniques. Because seasonal effects may influence food crop price forecasts, the first design rationally selects the preceding 12 observations to forecast food crop prices at timet. Accordingly, the first design used 12 self-predictor variables (i.e., yt−1,yt−2,yt−3,…, and yt−12) to serve as input variables and employed a single variable (i.e., yt) to serve as the output for ANN modeling. That is, the first design considered 12 input nodes and one output node for the ANN structures. This design is denoted as ANN1.The second input design for CI technique modeling is performed using FSM. For this proposed modeling, this study incorporated FSMs with ANN, SVR, and MARS to develop forecasting models. The proposed input design used ARIMA as an FSM to extract important self-predictor variables, which serve as the inputs for the ANN, SVR, and MARS models. Because the ARIMA model for rice price forecasting is obtained in (25), it indicates that Zt−1,Zt−2,Zt−7, and Zt−13 will influence Zt. Consequently, this study selected Zt−1,Zt−2,Zt−7, and Zt−13 to serve as self-predictor, or input, variables for the proposed integrated modeling on the rice price data series. In this study, this second input design is denoted as ARIMA-ANN.To forecast the rice prices, the two ANN models contain twelve and four input nodes for the first and second designs, respectively. The hidden nodes were selected as 10, 11, 12, 13, and 14 for the first design and 6, 7, 8, 9, and 10 for the second design. After performing the ANN modeling, the ANN1 design showed that the {12-14-1} structure provided the best results and minimum testing MAPE for rice prices. For the ARIMA-ANN design, the best structure was {4-7-1}. Three forecasting measurements for different settings of the ANN topologies for the two designs are shown in Table 3.Table 3
Forecasting measurements for two ANN designs (rice price forecasting).
ANN1 topology
MAE
RMSE
MAPE (%)
ARIMA-ANN topology
MAE
RMSE
MAPE (%)
{12-10-1}
15.438
16.648
3.939
{4-6-1}
10.384
12.053
2.654
{12-11-1}
14.737
16.022
3.765
{4-7-1}
10.365
12.037
2.650
{12-12-1}
15.001
16.262
3.825
{4-8-1}
10.374
12.044
2.652
{12-13-1}
15.182
16.499
3.869
{4-9-1}
10.379
12.062
2.653
{12-14-1}
13.915
15.320
3.557
{4-10-1}
10.389
12.039
2.655For modeling SVR on rice prices, this study adopted the input structure as ANN modeling. SVR modeling employed two designs for the input variables. The first design used 12 variables,yt−1,yt−2,yt−3,…, and yt−12 as the input variables and yt as the output variable. The second design used four self-predictor variables, Zt−1,Zt−2,Zt−7, and Zt−13 to serve as the input variables and Zt as the output variable. The first and second SVR designs are denoted as SVR1 and ARIMA-SVR, respectively. Because the parameter settings of C and ε often affect the performance of SVR modeling, an analytic parameter selection approach and grid search are used in this study [26]. Accordingly, we denote the term {ε, C}SVR as the best SVR parameter settings. After performing the SVR modeling, the SVR1 design reported that the parameter settings of {2−4, 211}SVR provided MAPE = 6.912% for rice price forecasting. For the second design, ARIMA-SVR, the parameter settings of {2−5, 29}SVR provided MAPE = 2.057%.For MARS modeling, this study also adopted the same input structure of ANN or SVR modeling. The first and second MARS designs are denoted as MARS1 and ARIMA-MARS, respectively. Table 4 lists the variable selection results and BFs after modeling rice price data using MARS1. In Table 4, the first column contains the variables that should be included in the model, the second column contains the relative importance of the variables that are listed in the first column, the third column contains the various BFs, and the final column contains the estimated coefficient values for the BFs that are listed in the third column. Based on the results in Table 4, it can be seen that three input variables (i.e., yt−12,yt−11, and yt−10) played important roles in building the MARS forecasting models. Their relative importance values (%) are 100, 13.1, and 12.1. The construction of the BFs can be expressed as follows. Consider the variable yt−12 as an example. We have two BFs (i.e., BF1 and BF2) to be considered. The observed values of BF1 are determined by the values larger than 0 or yt−12−367. The corresponding coefficient is estimated to be 0.985. The observed values of BF2 are determined by the values larger than 0 or 552.09−yt−12. The corresponding coefficient is estimated to be −1.258. Accordingly, the MARS forecasting model for rice prices can be expressed as follows:
(27)Y=542.484+0.985BF1−1.258BF2+1.539BF3−1.880BF4+0.524BF5−2.410BF6+1.783BF7+0.454BF8−0.234BF9+0.346BF10.Table 4
The results of MARS1 for rice price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
542.484
yt−12
100.0
BF1 = Max (0, yt−12 − 367)
0.985
BF2 = Max (0, 552.09 − yt−12)
−1.258
yt−11
13.1
BF3 = Max (0, yt−11 − 356)
1.539
BF4 = Max (0, yt−11 − 472.48)
−1.880
BF5 = Max (0, 585.95 − yt−11)
0.524
yt−10
12.1
BF6 = Max (0, yt−10 − 363)
−2.410
BF7 = Max (0, yt−10 − 472.48)
1.783
BF8 = Max (0, yt−10 − 543.14)
0.454
BF9 = Max (0, 623 − yt−10)
−0.234
BF10 = Max (0, yt−10 − 623)
0.346For ARIMA-MARS modeling, the variable selection results and BFs are summarized in Table5. In addition, it can be seen that three input variables (i.e., Zt−13,Zt−7, and Zt−1) played important roles in building the ARIMA-MARS forecasting models. The MARS forecasting model for rice prices can be described as follows:
(28)Y=−13.404+0.407BF1+23.543BF2−73.468BF3+49.716BF4+0.479BF5,where
(29)Zt=yt−yt−1.Table 5
The results of ARIMA-MARS for rice price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
−13.404
Zt−13
100.0
BF1 = Max (0, Zt−13 + 35)
0.407
Zt−7
45.4
BF2 = Max (0, Zt−7 − 0)
23.543
BF3 = Max (0, Zt−7 + 1.30)
−73.468
BF4 = Max (0, Zt−7 + 1.87)
49.716
Zt−1
28.5
BF5 = Max (0, −39.85 − Zt−1)
0.479


## 3.3. Forecast Modeling of Wheat Prices

Figure2 displays the original time plot for wheat prices. After modeling the wheat price data series with the ARIMA procedure, the parameter estimates are shown in Table 6. In Table 6, the notations “MA1,1,” “MA1,2,” “MA1,3,” and “MA1,4” correspond to the parameters θ1, θ7, θ12, and θ13 of the MA model. Table 7 displays the Ljung-Box test results, and a conclusion is made on the appropriateness of the underlying ARIMA model. Accordingly, (30) is suitable for modeling wheat prices.
(30)Zt=at−0.2402at−1+0.14629at−7+0.11171at−12+0.13721at−13,where
(31)Zt=yt−yt−1.Figure 2
Time plot for wheat price (unit: US dollars per metric ton).Table 6
Parameter estimates for wheat prices.
The ARIMA procedureConditional least squares estimation
Parameter
Estimate
Standard error
t value
Approx Pr > |t|
Lag
MA1,1
−0.2402
0.05601
−4.29
<0.0001
1
MA1,2
0.14629
0.05754
2.54
0.0115
7
MA1,3
0.11171
0.0599
1.86
0.0632
12
MA1,4
0.13721
0.06018
2.28
0.0233
13Table 7
Ljung-Box test results for ARIMA modeling of wheat prices.
Autocorrelation check of residuals
To Lag
Chi-square
DF
Pr > ChiSq
Autocorrelations
6
4.4
2
0.111
0.019
−0.045
0.024
−0.036
0.069
0.074
12
11.07
8
0.1977
−0.03
−0.126
0.008
−0.055
0.042
−0.001
18
14.41
14
0.4197
−0.035
−0.084
0.021
−0.038
0.021
0.006
24
19.99
20
0.4589
0.067
0.081
−0.033
−0.03
0.012
−0.064
30
30.52
26
0.2465
−0.007
−0.041
−0.11
−0.132
−0.009
0.026
36
33.96
32
0.3733
−0.025
0.025
−0.071
0.011
0.054
0.029In addition,at−j stands for the white noise at time t−j.For ANN modeling on wheat prices, this study also used two ANN designs, same as the structures for modeling rice prices. The first design used 12 self-predictor variables (i.e.,yt−1,yt−2,yt−3,…, and yt−12) to serve as input variables and a single variable (i.e., yt) to serve as the output variable. The second design considers the ARIMA as an FSM to extract important self-predictor variables, which serve as the inputs for the ANN models. The ARIMA model for wheat price forecasting is obtained in (30) and indicates that at−1,at−7,at−12, and at−13 will influence Zt.Becauseat−1=Zt−1−Zt−1∧, where Zt−1∧ is the forecast at time t−1, we observe that Zt−1 will influence Zt. Therefore, Zt−1 should be selected as an input variable as long as the component at−1 is in the MA model. Following the same logic, because the model (e.i., (30)) contains at−1,at−7,at−12, and at−13, we conclude that Zt−1,Zt−7,Zt−12, and Zt−13 will influence Zt. Accordingly, we have selected Zt−1,Zt−7,Zt−12, and Zt−13 to serve as input variables for the second proposed ANN modeling design for the wheat price data series.After performing ANN modeling on wheat price data, Table8 shows the corresponding forecast validity measures for different settings of ANN topologies for the two designs. As shown in Table 8, we observe that the ANN1 design with the {12-11-1} structure has the smallest MAPE. For the second, or ARIMA-ANN design, the {4-8-1} structure was associated with the smallest MAPE.Table 8
Forecast validity measures for two ANN designs (rice price forecasting).
ANN1 topology
MAE
RMSE
MAPE (%)
ARIMA-ANN topology
MAE
RMSE
MAPE (%)
{12-10-1}
12.121
13.521
5.448
{4-6-1}
11.833
14.336
5.374
{12-11-1}
11.752
13.230
5.285
{4-7-1}
11.865
14.353
5.388
{12-12-1}
12.084
13.515
5.423
{4-8-1}
8.984
10.864
4.077
{12-13-1}
11.984
13.429
5.385
{4-9-1}
8.984
10.883
4.078
{12-14-1}
12.120
13.516
5.457
{4-10-1}
11.845
14.346
5.379For modeling SVR on wheat prices, this study adopted the input structure as ANN modeling. The first design employed 12 variables,yt−1,yt−2,yt−3,…, and yt−12 as the input variables and yt as the output variable. The second design employed four self-predictor variables, Zt−1,Zt−7,Zt−12, and Zt−13 to serve as the input variables and Zt as the output variable. After performing SVR modeling on the wheat price data series, we obtained the parameter settings of {2−8, 210}SVR and associated MAPE = 4.661% for the first design. In addition, we had parameter settings of {2−7, 2−1}SVR and associated MAPE = 4.347% for the second design.For MARS modeling with the first design, Table9 presents the variable selection results and BFs. As shown in Table 9, we notice that three input variables (i.e., yt−12,yt−6, and yt−11) played important roles in building the MARS forecasting models. In addition, the MARS forecasting model for wheat prices can be expressed as follows:
(32)Y=360.973+0.765BF1−1.272BF2−0.185BF3+0.298BF4−0.218BF5.Table 9
The results of MARS1 for wheat price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
360.973
yt−12
100.0
BF1 = Max (0, yt−12 − 321.18)
0.765
BF2 = Max (0, 321.81 − yt−12)
−1.272
yt−6
7.2
BF3 = Max (0, yt−6 − 178.32)
−0.185
BF4 = Max (0, 178.32 − yt−6)
0.298
yt−11
3.2
BF5 = Max (0, yt−11 − 269.03)
−0.218For MARS modeling with the second design, the variable selection results and BFs are summarized in Table10. Furthermore, the MARS forecasting model for wheat prices is expressed as follows:
(33)Y=−13.665−1.105BF1+8.333BF2−13.064BF3+4.955BF4−0.155BF5.Table 10
The results of ARIMA-MARS for wheat price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
−13.665
Zt−7
100.0
BF1 = Max (0, Zt−7 − 26.58)
−1.105
Zt−13
63.3
BF2 = Max (0, Zt−13 + 26.91)
8.333
BF3 = Max (0, Zt−13 + 22.55)
−13.064
BF4 = Max (0, Zt−13 + 17.56)
4.955
Zt−1
14.4
BF5 = Max (0, Zt−1 + 10.65)
−0.155


## 3.4. Forecast Modeling of Corn Prices

Figure3 presents the original time plot for corn prices. Table 11 presents the parameter estimates after performing ARIMA modeling. In Table 11, the notations “AR1,1,” “AR1,2,” “AR1,3,” “AR1,4,” “AR1,5,” and “AR1,6” correspond to as parameters ϕ1, ϕ3, ϕ6, ϕ12, ϕ14, and ϕ24 of the AR model. In Table 11, we still contain the parameters of “AR1,2” and “AR1,4” in the model, although the absolute t values are less than or equal to 2 (i.e., the typical type I error is chosen as 0.05). The main reason is that those two parameters are important, and the Ljung-Box statistics indicate that the model is not appropriate if those two parameters are not involved in the underlying ARIMA model. Additionally, Table 12 demonstrates the Ljung-Box test results and a conclusion is made on the appropriateness of the ARIMA model. Thus, (34) is a suitable ARIMA model for modeling the corn price data series.
(34)Zt=0.17626Zt−1+0.09344Zt−3−0.145Zt−6−0.0925Zt−12−0.1476Zt−14−0.2149Zt−24+at,where
(35)Zt=yt−yt−1.Figure 3
Time plot for corn price (unit: US dollars per metric ton).Table 11
Parameter estimates for corn prices.
The ARIMA procedureConditional least squares estimation
Parameter
Estimate
Standard error
t value
Approx Pr > |t|
Lag
AR1,1
0.17626
0.05909
2.98
0.0031
1
AR1,2
0.09344
0.0651
1.44
0.1524
3
AR1,3
−0.145
0.06527
−2.22
0.0272
6
AR1,4
−0.0925
0.06619
−1.4
0.1636
12
AR1,5
−0.1476
0.06719
−2.2
0.029
14
AR1,6
−0.2149
0.07024
−3.06
0.0024
24Table 12
Ljung-Box test results for ARIMA modeling of corn prices.
Autocorrelation check of residuals
To Lag
Chi-square
DF
Pr > ChiSq
Autocorrelations
6
.
0
.
−0.011
0.082
0.026
0.002
−0.063
−0.027
12
6.71
6
0.3481
0.013
−0.032
0.001
0.012
0.091
−0.014
18
14.3
12
0.2822
−0.069
−0.036
0.042
0.113
0.038
−0.046
24
16.43
18
0.5624
0.014
−0.033
−0.015
0.061
−0.034
−0.015
30
29.11
24
0.2159
0.035
−0.073
0.103
−0.079
−0.025
−0.12
36
30.86
30
0.4225
−0.036
0
0.06
0.005
0.015
−0.001For ANN modeling on corn prices, while the first design used 12 input variables, the second design used six input variables extracted by using the FSM of ARIMA modeling. After performing ANN modeling, Table13 presents the corresponding forecast validity measures for different ANN topology settings for the two designs. From Table 13, we observe that the ANN1 design with the {12-14-1} structure has the smallest MAPE. For the second design, the {6-13-1} structure has the smallest MAPE.Table 13
Forecast validity measures for two ANN designs (corn price forecasting).
ANN1 topology
MAE
RMSE
MAPE (%)
ARIMA-ANN topology
MAE
RMSE
MAPE (%)
{12-10-1}
7.536
9.698
4.387
{6-10-1}
7.304
10.773
4.267
{12-11-1}
7.464
9.635
4.343
{6-11-1}
7.307
10.853
4.269
{12-12-1}
7.444
9.638
4.331
{6-12-1}
7.291
10.805
4.259
{12-13-1}
7.487
9.633
4.357
{6-13-1}
7.286
10.771
4.256
{12-14-1}
7.433
9.589
4.326
{6-14-1}
7.290
10.809
4.259In the modeling of SVR for the corn price data series, this study adopted the same input structure as ANN modeling. The first design employed 12 variables (i.e.,yt−1,yt−2,yt−3,…, and yt−12) as the input variables and yt as the output variable. The second design employed six self-predictor variables (i.e., Zt−1,Zt−3,Zt−6,Zt−12,Zt−14, and Zt−24), which were selected by using the FSM of ARIMA modeling to serve as the input variables and Zt as the output variable. After performing SVR modeling, we obtained the parameter settings of {2−9, 27}SVR and associated MAPE = 4.060% for the first design. In addition, we obtained parameter settings of {2−11, 211}SVR and associated MAPE = 4.215% for the second design.For MARS modeling on corn price data with the first design, Table14 illustrates the variable selection results and BFs. As shown in Table 14, we notice that two variables (i.e., yt−12 and yt−11) play important roles in building the MARS forecasting models. Consequently, the MARS forecasting model for corn prices can be described as follows:
(36)Y=262.148−1.354BF1+0.969BF2+0.342BF3.Table 14
The results of MARS1 for corn price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
262.148
yt−12
100.0
BF1 = Max (0, 268.79 − yt−12)
−1.354
BF2 = Max (0, yt−12 − 268.79)
0.969
yt−11
7.2
BF3 = Max (0, 295.29 − yt−11)
0.342For MARS modeling with the second design, Table15 lists the variable selection results and BFs. The MARS forecasting model for the corn prices is expressed as follows:
(37)Y=3.787−0.277BF1−0.745BF2−0.466BF3+0.175BF4.Table 15
The results of ARIMA-MARS for corn price forecasting.
Variable selection results
BF
Coefficient
Variable
Relative importance (%)
Intercept
3.787
Zt−24
100.0
BF1 = Max (0, 15.08 − Zt−24)
−0.277
Zt−12
74.0
BF2 = Max (0, Zt−12 − 13.9)
−0.745
Zt−1
44.8
BF3 = Max (0, Zt−1 − 11.59)
−0.466
Zt−6
17.7
BF4 = Max (0, Zt−6 + 7.36)
0.175


## 4. Forecasting Comparison

Several forecasting models were proposed to forecast the three major food crop prices in this study. These models include four single models (i.e., ARIMA, ANN1, SVR1, and MARS1) and three integrated models with FSMs (i.e., ARIMA-ANN, ARIMA-SVR, and ARIMA-MARS). Tables 16–18 present the forecasting results, as well as the MAE, RMSE, and MAPE values of the forecasting models for rice, wheat, and corn price forecasting, respectively. Low MAE, RMSE, or MAPE values are associated with better forecasting accuracy.Table 16
Performance comparison of various models for rice price forecasting.
Inputs
MAE (%)
RMSE (%)
MAPE (%)
(1) Single models
ARIMA
48.697
53.077
12.615
ANN1
yt−1~yt−12
13.915
15.320
3.557
SVR1
yt−1~yt−12
27.994
34.514
6.912
MARS1
yt−10,yt−11,yt−12
20.392
25.090
5.196
(2) Integrated models
ARIMA-ANN
Zt−1,Zt−2,Zt−7,Zt−13
10.365
12.037
2.650
ARIMA-SVR
Zt−1,Zt−2,Zt−7,Zt−13
7.958
9.786
2.057
ARIMA-MARS
Zt−1,Zt−7,Zt−13
11.244
15.566
2.916Table 17
Performance comparison of various models for wheat price forecasting.
Inputs
MAE (%)
RMSE (%)
MAPE (%)
(1) Single models
ARIMA
33.520
42.250
16.651
ANN1
yt−1~yt−12
11.752
13.230
5.285
SVR1
yt−1~yt−12
10.203
13.586
4.661
MARS1
yt−6,yt−11,yt−12
12.804
13.794
5.667
(2) Integrated models
ARIMA-ANN
Zt−1,Zt−7,Zt−12,Zt−13
8.984
10.883
4.078
ARIMA-SVR
Zt−1,Zt−7,Zt−12,Zt−13
9.779
11.371
4.347
ARIMA-MARS
Zt−1,Zt−7,Zt−13
10.643
13.937
4.809Table 18
Performance comparison of various models for corn price forecasting.
Inputs
MAE (%)
RMSE (%)
MAPE (%)
(1) Single models
ARIMA
14.783
20.047
8.784
ANN1
yt−1~yt−12
7.433
9.589
4.326
SVR1
yt−1~yt−12
7.003
9.058
4.060
MARS1
yt−11,yt−12
7.015
9.337
4.110
(2) Integrated models
ARIMA-ANN
Zt−1,Zt−3,Zt−6,Zt−12,Zt−14,Zt−24
7.286
10.771
4.256
ARIMA-SVR
Zt−1,Zt−3,Zt−6,Zt−12,Zt−14,Zt−24
7.190
9.588
4.215
ARIMA-MARS
Zt−1,Zt−6,Zt−12,Zt−24
6.495
9.102
3.785In comparison to the forecasting performance of the first design, or single models, in Tables16–18, we observe that the three CI models demonstrated better performance than the ARIMA models. The possible reason may be that ARIMA modeling is difficult for capturing nonlinear features in the food crop price data series. By reviewing these tables, we found that there is no best model for food crop price forecasting. For example, while the ANN1 model seems to possess better forecasting accuracy for rice price forecasting, the SVR1 model seems to have better forecasting accuracy for wheat and corn price forecasting.It can be clearly seen in Tables16–18 that the proposed integrated models possess better forecasting accuracy than the single models for most cases. For example, by reviewing Table 16, the proposed integrated ARIMA-ANN model is associated with MAE, RMSE, and MAPE values of 10.365%, 12.037%, and 2.650%, respectively, for rice price forecasting. These three performance values are smaller than the corresponding performance values for any one of the four single models. Take the proposed integrated ARIMA-SVR model as another example. By reviewing Table 17, the ARIMA-SVR model is associated with MAE, RMSE, and MAPE values of 9.779%, 11.371%, and 4.347%, respectively, for wheat price forecasting. These three performance values are also smaller than the corresponding performance values for any one of the four single models in Table 17. Thus, our proposed models, which were integrated with FSM, provide more accurate forecasting results than the single models.In addition, Table19 presents a comparison with respect to the overall percentage improvements (PIs) of forecasting accuracy for the proposed integrated models over the single models. The PIs of the MAE, RMSE, and MAPE are defined as follows:
(38)MAEPI=MAEsinglemodel−MAEintegratedmodelMAEintegratedmodel∗100%,RMSEPI=RMSEsinglemodel−RMSEintegratedmodelRMSEintegratedmodel∗100%,MAPEPI=MAPEsinglemodel−MAPEintegratedmodelMAPEintegratedmodel∗100%.Table 19
Improvement of the proposed integrated models over single models.
(1) Rice price forecasting
MAEPI (%)
RMSEPI (%)
MAPEPI (%)
(ii) Integrated ARIMA-ANN
ARIMA
369.822
340.945
376.038
ANN1
34.250
27.273
34.226
(ii) Integrated ARIMA-SVR
ARIMA
511.925
442.372
513.272
SVR1
251.772
252.688
236.023
(iii) Integrated ARIMA-MARS
ARIMA
333.093
240.977
332.613
MARS1
81.359
61.186
78.189
(2) Wheat price forecasting
(i) Integrated ARIMA-ANN
ARIMA
273.108
288.223
308.313
ANN1
30.810
21.562
29.598
(ii) Integrated ARIMA-SVR
ARIMA
242.775
271.562
283.046
SVR1
4.336
19.483
7.223
(iii) Integrated ARIMA-MARS
ARIMA
214.949
203.152
246.247
MARS1
20.304
−1.029
17.842
(3) Corn price forecasting
(i) Integrated ARIMA-ANN
ARIMA
102.896
86.124
106.391
ANN1
2.018
−10.973
1.645
(ii) Integrated ARIMA-SVR
ARIMA
105.605
109.089
108.399
SVR1
−2.601
−5.525
−3.677
(iii) Integrated ARIMA-MARS
ARIMA
127.606
120.253
132.074
MARS1
8.006
2.579
8.587From Table19, it is obvious that positive PIs can be achieved by using the proposed integrated models. For example, in the rice price forecasting, the MAPEPI of the proposed ARIMA-ANN model over the two single models, ARIMA and ANNS, were 376.038% and 34.226%, respectively. The corresponding MAEPI and RMSEPI were 369.822% and 34.250% and 340.945% and 27.273%, respectively. Apart from the rice price forecasting, most of the MAEPI, RMSEPI, and MAPEPI are positive large numbers for wheat and corn price forecasting. Accordingly, considerable forecasting accuracy improvements are achieved through the use of the proposed integrated models. Furthermore, we also found five negative PIs in Table 19.For wheat price forecasting, theRMSEPI of the ARIMA-MARS over MARS1 was −1.029%. For corn price forecasting, the RMSEPI of the ARIMA-ANN over ANN1 was −10.973%. The three remaining negative values occurred in corn price forecasting when using the ARIMA-SVR model, and the MAEPI, RMSEPI, and MAPEPI were −2.601%, −5.525%, and −3.677%, respectively. A negative PI value implies that the forecasting performance declined. However, the magnitudes of those five negative values are small, and they have minor effects on forecasting performance. Additionally, the associated MAPEPI for those five negative PIs were only one negative. That is, if a forecasting model has the minimum MAPE, this does not mean it has the minimum MAE or RMSE.


## 5. Discussion

While the aforementioned description focuses on accuracy comparisons for each individual food crop price forecast, the following discussion provides the comparison regarding overall performance for food crop price forecasting as a whole.Table20 lists the average PIs of the proposed integrated models over their common element, the ARIMA models. Figure 4 presents the average PIs of MAE, RMSE, and MAPE by employing the proposed integrated models over single ARIMA models. As shown in Figure 4, we notice that considerable accuracy improvements can be reached by using our proposed approaches. Additionally, Table 21 reports the average PIs of MAE, RMSE, and MAPE by using the proposed integrated models over another element, other than ARIMA element. Figure 5 displays the satisfied average PIs by using the proposed integrated models.Table 20
PIs of the whole food crop prices by using the proposed integrated models over the single ARIMA models.
MAE (%)
RMSE (%)
MAPE (%)
(1) ARIMA-ANN over ARIMA
Rice
369.822
340.945
376.038
Wheat
273.108
288.223
308.313
Corn
102.896
86.124
106.391
Average PIs
248.609
238.431
263.581
(2) ARIMA-SVR over ARIMA
Rice
511.925
442.372
513.272
Wheat
242.775
271.562
283.046
Corn
105.605
109.089
108.399
Average PIs
286.768
274.341
301.572
(3) ARIMA-MARS over ARIMA
Rice
333.093
240.977
332.613
Wheat
214.949
203.152
246.247
Corn
127.606
120.253
132.074
Average PIs
225.216
188.127
236.978Figure 4
Average PIs by using the proposed integrated models over single ARIMA models.Table 21
PIs of the whole food crop prices by using the proposed integrated models over another element, other than ARIMA element.
MAE (%)
RMSE (%)
MAPE (%)
(1) ARIMA-ANN over ANN
Rice
34.250
27.273
34.226
Wheat
30.810
21.562
29.598
Corn
2.018
−10.973
1.645
Average PIs
22.359
12.621
21.823
(2) ARIMA-SVR over SVR
Rice
251.772
252.688
236.023
Wheat
4.336
19.483
7.223
Corn
−2.601
−5.525
−3.677
Average PIs
84.502
88.882
79.856
(3) ARIMA-MARS over MARS
Rice
81.359
61.186
78.189
Wheat
20.304
−1.029
17.842
Corn
8.006
2.579
8.587
Average PIs
36.556
20.912
34.873Figure 5
Average PIs by using the proposed integrated models over single ANN, SVR, and MARS models.In addition, Figure6 displays the forecasts of rice prices, obtained by using ARIMA, ANN1, and the proposed ARIMA-ANN models, for the twelve testing records. One can see that the forecasts of the ARIMA-ANN model are closest to the actual observations. The forecasts of the basic ARIMA model are relatively far away from the actual observations. Additionally, the prediction accuracy of ANN1 is better than that of the ARIMA model. Regarding the ANN1 design, we observe that the input vectors I1,I2,…,I12 are equivalent to yt−1,yt−2,yt−3,…,yt−12 and that the output O can be obtained by performing a nonlinear functional mapping, as shown in (10). Regarding the ARIMA-ANN mechanism, we observe that the input vectors I1,I2,I3,andI4 are characterized by Zt−1,Zt−2,Zt−7,andZt−13 and that the output O can be obtained by performing a nonlinear functional mapping, as shown in (10).Figure 6
Forecasts of rice prices versus actual observations for the 12 testing records (with the use of ARIMA, ANN1, and the proposed ARIMA-ANN models).Figure7 presents the forecasts of rice prices, obtained by using ARIMA, SVR1, and the proposed ARIMA-SVR models, for the twelve testing records. Regarding the SVR1 design, we observe that the input vectors I1,I2,…,I12 are characterized by yt−1,yt−2,yt−3,…,yt−12 and that the output O can be obtained by performing a nonlinear functional mapping, as shown in (17). Regarding the ARIMA-SVR mechanism, we observe that the input vectors I1,I2,I3,andI4 are characterized by Zt−1,Zt−2,Zt−7,andZt−13 and that the output O can be obtained by performing a nonlinear functional mapping, as shown in (17).Figure 7
Forecasts of rice prices versus actual observations for the 12 testing records (with the use of ARIMA, SVR1, and the proposed ARIMA-SVR models).Figure8 shows the forecasts of rice prices, obtained by using ARIMA, MARS1, and the proposed ARIMA-MARS models, for the twelve testing records. Regarding MARS1 modeling, we observe that the input vectors I1,I2,…,I12 are characterized by yt−1,yt−2,yt−3,…,yt−12 and that the output O can be obtained by performing a nonlinear functional mapping, as shown in (22). Regarding the ARIMA-MARS design, we observe that the input vectors I1,I2,I3,andI4 are characterized by Zt−1,Zt−2,Zt−7,andZt−13 and that the output O can be obtained by performing a nonlinear functional mapping, as shown in (22). By observing Figures 7 and 8, we can see that the forecasts of the ARIMA-SVR and ARIMA-MARS models are closest to the actual observations. Both figures illustrate that the forecasts of the ARIMA models are relatively far away from the actual observations.Figure 8
Forecasts of rice prices versus actual observations for the 12 testing records (with the use of ARIMA, MARS1, and the proposed ARIMA-MARS models).Similar implications were obtained for the cases of wheat and corn price predictions. That is, the forecasts of the proposed integrated models are closer to the actual observations. The forecasts of the ARIMA model are far away from the actual observations. These findings can be observed in Figures9–14.Figure 9
Forecasts of wheat prices versus actual observations for the 12 testing records (with the use of ARIMA, ANN1, and the proposed ARIMA-ANN models).Figure 10
Forecasts of wheat prices versus actual observations for the 12 testing records (with the use of ARIMA, SVR1, and the proposed ARIMA-SVR models).Figure 11
Forecasts of wheat prices versus actual observations for the 12 testing records (with the use of ARIMA, MARS1, and the proposed ARIMA-MARS models).Figure 12
Forecasts of corn prices versus actual observations for the 12 testing records (with the use of ARIMA, ANN1, and the proposed ARIMA-ANN models).Figure 13
Forecasts of corn prices versus actual observations for the 12 testing records (with the use of ARIMA, SVR1, and the proposed ARIMA-SVR models).Figure 14
Forecasts of corn prices versus actual observations for the 12 testing records (with the use of ARIMA, MARS1, and the proposed ARIMA-MARS models).In this study, we collected 309 records to build various models for the prediction of the prices of rice, wheat, and corn. One assumption for building the ARIMA models is that their structures will remain unchanged over time. Therefore, even if more sample data become available, there is no need to rebuild the ARIMA model for our proposed integrated models. This is a significant advantage over existing models. In other words, the ARIMA feature selection procedure only needs to be performed once. Even if more sample data becomes available, we can still use the same input variables for our proposed ARIMA-ANN, ARIMA-SVR, and ARIMA-MARS models. Some feature selection modeling processes may need to be performed repeatedly when more sample data become available, which may be a time-consuming task. It is true that, after performing feature selection with the ARIMA, we can increase the accuracy of our forecasts by performing ANN, SVR, and MARS modeling with the additional sample data.


## 6. Conclusion

Humans survive based on food crops. We eat significant quantities of rice, wheat, corn, and other simple crops to maintain energy and good health. Accordingly, food crop price forecasting is very important and has drawn considerable attention in recent decades.Typical CI forecasting techniques require proper explanatory variables to make predictions. However, proper explanatory variables are difficult to capture, and it is infeasible to obtain the future values of these variables. Therefore, we proposed integrated ARIMA-ANN, ARIMA-SVR, and ARIMA-MARS models in order to perform forecasting for three important food crop prices. The role of the ARIMA element is as an FSM that can capture important self-predictor variables. Rather than using unavailable explanatory variables, the self-predictor variables serve as the inputs for ANN, SVR, and MARS modeling. The experimental results reveal that the proposed integrated models are desirable alternatives for food crop price forecasting, because they all have excellent forecasting performance. Most importantly, the main contribution of the proposed models is their ability to provide predictions of food crop prices without requiring extensive effort to obtain the future values of explanatory variables.One limitation of our models is that the techniques for identifying the correct ARIMA model from the variety of possible models may be unintuitive and computationally expensive. However, the modeling in this study can be used as a guideline for developing forecasting models for other food crop price data series. Additionally, because MARS is effective for selecting important variables for predicting response variables, attempting to extend the FSM of MARS modeling may be a valuable future research direction. Finally, the integrated models may also be combined with other CI techniques, such as extreme learning machines, time delay neural networks, or artificial immune systems, which may be worthy of investigation in the future.



---

*Source: 1910520-2018-07-10.xml*