*1.2. Time Series and Classical Methods*

Time series analysis studies the past behavior of historical series using different methods (Table 1). It verifies trends, seasonality, and randomness in a dataset in two ways: stationary, when observations oscillate around a central horizontal axis; and non-stationary when oscillates around changing values [42,43]. The most appropriate model for a specific dataset is the coefficient of determination (R), the mean absolute error (MAE), and the mean squared error (MSE) [42,43].



Source: adapted from [42,43].

The coefficient of determination (Equation (4)) measures the linear regression adjustment, which aims to explain the relationship of the variables. The closer this number is to one, the more fitted is the model. However, a measure higher than 0.7 is satisfactory [25,42,43]:

$$\mathbf{R} = \frac{\sum \begin{pmatrix} \hat{y} - \ \boldsymbol{y} \end{pmatrix}^2}{\sum \begin{pmatrix} \boldsymbol{y} - \ \overline{\boldsymbol{y}} \end{pmatrix}^2}. \tag{4}$$

The coefficient of determination is calculated based on the ratio between the explained and the total variance where *y* represents the real value of the series, *y*ˆ is the expected value (value of the regression line approaching the actual value), and *y*¯ is the average value of the series.

Note that the variance is the difference between the expected value and the mean, and the total variance is the difference between the original and mean value [25,42,43]. The MAE and MSE are calculated according to Equations (2) and (3), where *n* is the number of elements in the series.

$$\text{MAE} = \frac{\sum |y - \hat{y}|}{n} \,\text{}\tag{5}$$

$$\text{MSE} = \frac{\sum \left( y - \hat{y} \right)^2}{n}. \tag{6}$$

Finally, functions with error values close to 0 are the most effective in predicting future values. These time series applications are described in the Results and Discussion section.

### **2. Materials and Methods**

### *2.1. Dataset*

To perform this study, we collected data from the Food and Agriculture Organization of the United Nations (FAO) [44] regarding harvest area (million hectares), yield (tons per hectare), and production (million tons) between 1961–2016. The dataset was imported from MS© Excel 2016 spreadsheet to Matlab© R2017b arrays. However, the period from 1961–1966 was used only for delay configuration, and it was not plotted on the time series [29,45].

Firstly, we conducted Time Series Analysis. The historical series was extracted and processed in MS© Excel 2016 spreadsheet format generating graphs with trend lines. Tables 2–4 present the formulas.

**Table 2.** Harvest area.


*x* = timestep(year). *y* = million hectares.

**Table 3.** Yield.


*x* = timestep(year). *y* = tons per hectare.


**Table 4.** Production.

*x* = timestep(year). *y* = million tons.

Secondly, we used neural networks toolbox of the Matlab©R2017b software to create, train, and validate the ANN model—we tested with 10 neurons and six delays (Figure 2).

**Figure 2.** Neural network created in Matlab 2017b.

We adopted the Nonlinear Autoregressive Network with External Input-NARX type because it has proven to be the most effective and accurate solution for multivariable data series [27,46–48]. The NARX network applies historical input data with time delay operators [9]. We used 70% of data for training, 15% for validation, and 15% for testing. We defined the percentage based on k-cross validation that utilizes efficiently the learning abilities of the ANN model [49], and data are distributed randomly by NARX [46]. Moreover, we adopted the Levenberg–Marquardt algorithm for backpropagation due to being the fastest supervised algorithm for training and widely used for time series prediction in the ANN model [25,30,46].

For harvested area (target), we used yield and production as input variables; for yield (target), we used harvested area and production as input variables; for production (target), harvested area and yield were used as input variables.

After that, Matlab© R2017b provided algorithms for closed-loop form simulation (named multistep prediction). This type of simulation is important to verify the ability of the networks to make predictions (calculation of errors) [25]. Figure 3 shows the overall flowchart of the ANN model.

### *2.2. Model Classification*

The differences between the original and predicted values were computed using MAE and MSE, even for ANN. We compared classical models and neural networks where the errors of each model were sorted from lowest to highest. However, regression measures were sorted from highest to lowest. Depending on the use of two measures of error the weighted average was used (Equation (7)):

$$\text{Rank} = \frac{(\text{MAE} \times 0.5) + (\text{MSE} \times 0.5) + (\text{R} \times 1)}{2} \tag{7}$$

where MAE is mean absolute error, MSE is mean squared error, and R is the coefficient of determination.

**Figure 3.** Flowchart of the ANN model.

### **3. Results and Discussion**

*3.1. Time Series Analysis Using Classical Predictive Methods (Functions)*

### 3.1.1. Harvested Area

The first application of classical methods for prediction uses the time series for harvest area (target). Figure 4 illustrates the 1967–2016 timespan in million hectares.

The harvested area raised continuously, mainly after the harvest of 1997 (Timestep 31), and reached around 33 million hectares in 2016 (Timestep 50). Looking back over 1967 year (Timestep 1), the planted soybean area was 2% (around 0.6 million hectares) of the area planted in 2016. There has been a more than a 50-fold increase while the US, the main Brazilian competitor, had in 1967 54% of the current planted area [22].

Regarding the fit of functions, polynomial and power were more effective in predicting harvest area considering R, MAE, and MSE (Table 5).


**Table 5.** Effective functions for forecasting harvested area.

R is the coefficient of determination (value between 0 and 1), MAE is the mean absolute error (value in millions of hectares), and MSE is the mean squared error (value in millions of hectares).

*Agriculture* **2020**, *10*, 475

**Figure 4.** Original time series for harvested area.

## 3.1.2. Yield

The second application of classical methods of predication verifies soybean yield (target). Figure 5 depicts the time series analysis results considering data from 1967–2016 in hectares.

**Figure 5.** Original time series for yield.

The average yield over the 50-year period was approximately two tons per hectare. The lowest value was 0.9 hectares in 1968 (Timestep 2), and the highest value 3.1 tons per hectare in 2011 (Timestep 45). During the 1967 (Timestep 1) crop season, the national yield was approximately 1.2 ton per hectare, which corresponds to less than half of the yield in 2016 (Timestep 50). In addition, compared with the US in the same period, the lowest value was 1.6 ton per hectare in 1974 and the highest value 3.5 in 2016 [22]. This means that the Brazilian soybean yield varies 158% against 119% of the US.

The yield increase in soybean production was affected by changes in production processes and the use of new technologies. Moreover, genetic improvements of soybean created a variety of grain with better adaptation to the climate that affected productivity [12]. According to Pereira [50], the last 30 years have demonstrated innovative solutions in food production, such as new crop varieties and new irrigation techniques. However, Dani [12] argue that the evolution has generally been technological without improving governance processes causing logistics issues throughout the supply chain. Logistics has a huge impact on soybean production and directly affect the trade [51].

Related to the fit of the functions, linear and polynomial predict more precisely soybean yield given that R, MAE, and MSE (Table 6).


**Table 6.** Effective functions for forecasting yield.

R is the coefficient of determination (value between 0 and 1), MAE is the mean absolute error (value in tons per hectare), and MSE is the mean squared error (value in tons per hectare).

### 3.1.3. Production

Finally, we use classic methods to predict production. Figure 6 shows the time series analysis for soybean production from 1967–2016 in millions of tons.

**Figure 6.** Original time series for production.

Brazilian soybean production raised 130-fold, moving from 700 thousand tons in 1967 (Timestep 1) to 96.3 million tons in 2016 (Timestep 50). In the same period, the US soybean production moves from 26.5 million tons to 116.9 million tons [22]. In 2019, the Brazilian soybean production was 25% higher [21] than 2016.

Furthermore, we identify that in 1970 (Timestep 4) the production was 1.5 million tons, but, in 1977 (Timestep 11), it reached 12.5 million tons. This great expansion of soybean cultivation occurs due to the expansion of international demand and the national soybean oil industry [52].

Regarding soybean production, polynomial and power were more effective considering R, MAE, and MSE (Table 7).


**Table 7.** Effective functions for forecasting production.

### *3.2. ANN Model*

### 3.2.1. Training, Validation, and Testing of Neural Network

Harvested area training reached an optimal value for the regression and correlation among variables after nine interactions (Figure 7). The training procedure stops when the performance on the test data does not improve following a fixed number of training iterations [39]. The main purpose of the training phase is to find the optimal set of weights for the ANN model where the error is minimized [35].

**Figure 7.** Regression and correlation of the harvested area using Matlab R2017b.

The training, validation, and testing indicate that the network learned from the data (R > 0.99). Moreover, the fit was well-aligned, which means the model has a good capacity for generalization and prediction.

Yield training reached an optimal value for the regression and correlation among variables after 12 iterations and pose an R higher than 0.9 (Figure 8).

R is the coefficient of determination (value between 0 and 1), MAE is the mean absolute error (value in millions of tons), and MSE is the mean squared error (value in millions of tons).

**Figure 8.** Regression and correlation for yield using Matlab R2017b

Yield presents correlation and regression results similar to the other two networks. However, fit shows a reasonable alignment representing a capacity of the network generalize and predict.

Finally, the production network was trained and after nine interactions reached an optimal value for the regression and correlation among variables (Figure 9).

**Figure 9.** Regression and correlation for production using Matlab R2017b.

The network shows an excellent rate of learning, with reasonable values of alignment. However, the validation and test pose deviations on fit. The overall results present proper alignment confirming its ability to generalize and make predictions.

Given that there are three networks, the harvested area indicated the best results, followed by production and yield.

3.2.2. Time Series Results with an Artificial Neural Network

Figures 10–12 depict the results of the time series generated by the neural networks in closed-loop form (multistep prediction). The blue line (target) represents the original data, and the red line (prediction) represents the obtained values for each period.

Neural network prediction shows better adjustment to the original data than time series analysis. In other words, these predictions show a smoother follow-up. The trend lines of the classical models follow the randomness of the series increasing the error between the original and predicted values. The base graphic Figure 10 presents the error of the prediction in millions of hectares, in Figure 11 in tons per hectare and Figure 12 in millions of tons.

The classical models are based on elements dependent on the analysis of their predecessors. On the flipside, ANN is a generalization of the classical models, where an element to be predicted also depends on the previous elements of other related time series [27].

**Figure 10.** Multistep prediction of harvested area using Matlab R2017b.

**Figure 11.** Multistep prediction of yield using Matlab R2017b.

**Figure 12.** Multistep prediction of production using Matlab R2017b.

3.2.3. Comparison between Artificial Neural Networks and Time Series Classical Models

Considering R, MAE, and MSE, Tables 8–10 present the ranking of ANN model versus classical functions for forecast harvested area, yield, and production, respectively.

Considering the results, Artificial Neural Networks ranked first for predicting harvested area and production and third place to yield. The polynomial model ranked second in all three series showing the reliability of the model to estimate future values. The logarithmic model is the least fitted and should be discarded for these series.


**Table 8.** Effective functions versus ANN for forecasting harvested area.

R is the coefficient of determination (value between 0 and 1), MAE is the mean absolute error (value in millions of hectares), and MSE is the mean squared error (value in millions of hectares).

**Table 9.** Effective functions versus ANN for forecasting yield.


R is the coefficient of determination (value between 0 and 1), MAE is the mean absolute error (value in tons per hectare), and MSE is the mean squared error (value in tons per hectare).


**Table 10.** Effective functions versus ANN for forecasting production.

R is the coefficient of determination (value between 0 and 1), MAE is the mean absolute error (values in millions of tons), and MSE is the mean squared error (values in millions of tons).

Based on these results, it is possible to infer that predictive capabilities of the developed ANN model are efficient to soybean prediction with short data time series. This fact confirms the superior performance of the ANN model against classical methods. Similar results are obtained for Nedic et al. [53] when compared to an ANN model with classical statistical models to predict traffic noise.

ANN has been recognized as a valuable predictive tool due to its ability to learn, adapt, and generalize the results of a sample of noise data and are more effective and flexible than conventional statistics for dealing with nonlinearity [54]. Moreover, there is a tendency towards the adopting of artificial intelligence in decision models.

### **4. Conclusions**

This study compares classical methods of time series prediction with Artificial Neural Networks using Brazilian soybean harvest area, yield and production from 1961–2016. The results indicate that ANN is the best approach to predict soybean harvest area and production while classical linear function remains more effective to predict soybean yield. However, ANN is a reliable model to predict using time series and can help farmers, government, and trading companies anticipate the soybean world offer to organize efficiently logistics resources and public policies.

Our results confirm the important role of neural networks in dealing with agriculture issues as showed in previous studies in the literature [8,10,35,39]. The R value above 0.9 confirms the high performance of the model. Nevertheless, regarding the agriculture concerns about low availability for planting areas, yield, and production [4–6], your results demonstrated that, at least in case of the soybean, this is not a concern.

Furthermore, we can conclude that the ANN model can be effective even using a short time series—that, in our case, was 50 years. This fact reveals a robustness of the model. However, despite the advantages of the ANN model, classical methods also can produce very good models. A comparison in other agriculture commodities can be made to confirm or refuse the behavior presented in a soybean case.

Finally, we also suggest for further studies to combine neural networks in hybrid systems using, for example, ANN and Fuzzy Logic, similar to that proposed by [41]. Literature has shown that hybrid systems are more efficient. The goal is to achieve a synergy between hybrid systems to compensate for the disadvantage of one by the advantage of another [55–57].

**Author Contributions:** Conceptualization, E.R.A., J.G.M.d.R., O.V., P.L.d.O.C.N., and R.C.T.; methodology, E.R.A., J.G.M.d.R., O.V., and P.L.d.O.C.N.; software, E.R.A.; validation, J.G.M.d.R., O.V., P.L.d.O.C.N., R.C.T., and A.E.d.S.; formal analysis, J.G.M.d.R., O.V., P.L.d.O.C.N., M.d.O.M., and R.C.T.; investigation, E.R.A., J.G.M.d.R., and R.C.T.; resources, E.R.A. and J.G.M.d.R. Data curation, E.R.A. and J.G.M.d.R.; writing—original draft preparation, E.R.A. and J.G.M.d.R.; writing—review and editing, J.G.M.d.R., O.V., P.L.d.O.C.N., M.d.O.M., A.E.d.S., and R.C.T.; visualization, E.R.A. J.G.M.d.R., O.V., P.L.d.O.C.N., M.d.O.M., A.E.d.S., and R.C.T.; supervision, J.G.M.d.R., O.V., and P.L.d.O.C.N.; project administration, E.R.A.; funding acquisition, M.d.O.M., A.E.d.S. and R.C.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior Grant No. 0001.

**Acknowledgments:** The authors would like to thank Sunsetti Treinamentos e Serviços and Universidade Paulista UNIP for the financial incentives.

**Conflicts of Interest:** The authors declare no conflict of interest.
