Prediction of Fruit Production in India: An Econometric Approach

Descriptive statistics are generally associated with the central tendency and dispersion measure used to summarize the data. Mean, standard deviation (SD), standard error (SE), coefficient of variation (CV%), skewness, kurtosis, and cumulative growth rate (CGR%) were considered in our study.

A time series analysis is a process of analyzing a set of data that was compiled over time. The research contributes to constructing a statistical model that predicts the behavior of data series. Our work was driven to build a time series model, namely, an ARIMA model, and Holt's nonlinear ETS model, to forecast fruit production (apple, banana, grape, pineapple, and mango) in India.

The Box–Jenkins method consists of four steps: identification, estimate, diagnostic testing, and forecasting. Before implementing the Box–Jenkins model, it is crucial to assess the stationarity of the time series and the significant seasonality that must be modeled (Ray & Bhattacharyya 2020a; Ray et al. 2021; Raghav et al. 2022). To check the stationarity of time series variables, the autocorrelation function (ACF) and partial autocorrelation function (PACF) are used (Ray & Bhattacharyya 2020b), as well as some popular statistical tests such as augmented Dickey–Fuller (1979), Phillips–Perron (1988), and Kwiatkowski–Phillips–Schmidt–Shin (1992). The ACF and PACF functions estimate the model order in the identification step and build ARIMA (p, d, q); where order p is denoted as the autoregressive model, d is integrated, and q is denoted as the moving average model. The ARIMA (p, d, q) model is denoted as ∅(β)(1 − β)^dZ_t = θ(β)ɛ_t; where ∅(β) = 1 − ∅₁β − ∅₂β² − ∅₃β³ − ⋯ … − ∅_pβ^p (AR model with order p); θ(β) = 1 − θ₁β − θ₂β² − θ₃β³ − ⋯ … − θ_qβ^q (MA model with order q); Zt=(1−β)dεt∇d {{Z}_{t}}=\frac{{{\left( 1-\beta \right)}^{d}}{{\varepsilon }_{t}}}{{{\nabla }^{d}}} (ARIMA model; β is the backshift operator); and ɛ_t ~ WN(0, σ²) (error term follows a normal distribution with zero mean and constant variance).

After determining the model order, the parameter estimation process is calculated using maximum likelihood approaches. The optimal model is chosen using model selection criteria such as Akaike (1974) information criteria, root mean square error (RMSE), and mean absolute error (MSE) (Mishra et al. 2022). Finally, the in-sample and out-of-sample forecasts are used to assess the validity and prediction of the chosen model in the last stage (Ray et al. 2016b). Because most real-time series have a growing and decreasing tendency, they are classified as nonstationary. The first-order differencing approach was applied to make the data series stationary. The normality of residuals is checked using the Jarque–Bera test (Ray et al. 2016b), the independence of residuals is checked using the Box–Ljung test, and the constant variance of residuals is checked using the ARCH LM test. Root mean squared error was used to assess the accuracy of the anticipated values (RMSE). The same procedure was carried out to predict the rainfall in Western Himalayan regions using monthly, seasonal, and annual data series.

Since most time series contain volatility and do not follow a linear trend, it may be beneficial to incorporate components that consider these characteristics in the estimated model. The nonlinear models of Holt cover systemic development that combines exponential smoothing (ETS) models to form a nonlinear dynamic model. Analysis of these models using state-space-based probability calculations supports model selection and prevision standard error calculation (Hyndman et al. 2002).

The model uses three primary time series components: trend (T), seasonal (S), and error (E). It reflects the long-term trend of time series movement, which is the unpredictable part of the time series. The annual data series was used in this study. The components we need are combined in our model in various additive and multiplicative combinations to produce y_t. We have an additive model y_t = T+E or a multiplicative model such as y_t = T×E; where the individual components of the model are given as follows: E [A, M]; T [N, A, M, AD, MD]; S [N, A, M]; where N – none, A – additive, M – multiplicative, AD – additive dampened, and MD – multiplicative dampened (damping uses an additional parameter to reduce the impacts of the trend over time). Accordingly, the models that we are interested in estimating can be found (after selecting S [N]) in Table A.

The best model was selected based on the following goodness of fit criteria: MAE=∑t=1n|yt^−yt|n; RMSE=∑t=1n(yt^−yt)2n;MAPE=1n∑t=1n|yt^−ytyt|×100;Theil=∑t=1n(yt^−yt)2n∑t=1ny^t2n+∑t=1nyt2n;BP=(y^¯−y¯)2∑t=1n(yt^−yt)2n;VP=(sy^−sy)2∑t=1n(yt^−yt)2n; CVP=2(1−r)sy^sy∑t=1n(yt^−yt)2n; \begin{array}{*{35}{l}}\text{MAE}=\frac{\sum\nolimits_{t=1}^{n}{\left| \widehat{{{y}_{t}}}-{{y}_{t}} \right|}}{n};\ \text{RMSE}=\sqrt{\frac{\sum\nolimits_{t=1}^{n}{{{\left( \widehat{{{y}_{t}}}-{{y}_{t}} \right)}^{2}}}}{n}}; \\ \text{MAPE}=\frac{1}{n}\sum\nolimits_{t=1}^{n}{\left| \frac{\widehat{{{y}_{t}}}-{{y}_{t}}}{{{y}_{t}}} \right|\times 100;} \\ \text{Theil}=\frac{\sqrt{\frac{\sum\nolimits_{t=1}^{n}{{{\left( \widehat{{{y}_{t}}}-{{y}_{t}} \right)}^{2}}}}{n}}}{\sqrt{\sum\nolimits_{t=1}^{n}{{}^{\hat{y}_{t}^{2}}\!\!\diagup\!\!{}_{n}\;}}+\sqrt{\sum\nolimits_{t=1}^{n}{{}^{y_{t}^{2}}\!\!\diagup\!\!{}_{n}\;}}};\text{BP}=\frac{{{\left( \bar{\hat{y}}-\bar{y} \right)}^{2}}}{\sum\nolimits_{t=1}^{n}{\frac{{{\left( \widehat{{{y}_{t}}}-{{y}_{t}} \right)}^{2}}}{n}}}; \\ \text{VP}=\frac{{{\left( {{s}_{{\hat{y}}}}-{{s}_{y}} \right)}^{2}}}{\sum\nolimits_{t=1}^{n}{\frac{{{\left( \widehat{{{y}_{t}}}-{{y}_{t}} \right)}^{2}}}{n}}};\ \text{CVP}=\frac{2\left( 1-r \right){{s}_{{\hat{y}}}}{{s}_{y}}}{\sum\nolimits_{t=1}^{n}{\frac{{{\left( \widehat{{{y}_{t}}}-{{y}_{t}} \right)}^{2}}}{n}}}; \\ \end{array} where MAE is the mean absolute error index; RMSE is the root mean square error index; MAPE is the mean absolute percentage error; Theil is the Theil index; BP is the bias proportion; VP is the variance proportion, and CVP is the co-variance proportion. Notably, the BP index determines how far the mean of forecast values is from the mean of the actual values, the VP index measures the proportion of variation in forecast values from the variation in actual values, and the CVP measures the remaining unsystematic errors.

Moreover, ŷ_t and y_t are the prediction and actual series, respectively; y^¯ {\bar{\hat{y}}} and y¯ {\bar{y}} are the mean of forecast and actual series respectively; s_ŷ and s_y are the standard deviations of forecast and actual series, respectively.

To support the decision concerning validity of forecasting, Diebold–Mariano, Giacomini–White, and Clark–West tests were used to compare ARIMA and Holt's nonlinear ETS model using Eviews 12 and 13, and Stata 17 software.

The first step in this analysis is to look at the behavior of the data series through descriptive statistics (Table 1) from 1961 to 2020, which show that bananas and mangoes were the most abundant fruits in India during the study period. From Table 1, one can see that the production of apples has a range from 90 to 2891 thousand tons, with an average 1153.465 thousand tons. Banana production has increased from 2257 to 31504 thousand tons from the study period, registering the average as 12282.65 thousand tons. Grape production has ranged from 70 to 3125 thousand tons with an average of 889.9539 thousand tons. The mango production has registered from 6988 to 25631 thousand tons with an average and standard deviation 11091.62 and 4863.780. The pineapple production has increased from 76 to 1984 thousand tons with an average of 904.2272 thousand tons. Furthermore, according to the table, the grape variation series has the most significant change (104.37%), followed by bananas (82.68%). Besides, the Skewness index value is positive in all cases, showing a steady increase in production over the research period. Furthermore, except for mangoes, the Kurtosis index found a platykurtic distribution with index values near 3, indicating that our data has no outliers. The Jarque–Bera test was used to check the normality behavior of the series: banana, grape, and mango series do not follow a normal distribution, as the p-value is less than 0.05. On the other hand, visualized actual data (Figs. 1, 2) shows all the series had a linear trend with little volatility, which means that all the series are not stationary at their level series; for this reason, the next step in this study is to examine the stationarity behavior of the series using different tests.

Figure 1.

Figure 2.

Before creating the model, the stationarity of the data must be assessed after examining the nature of the data series. To test the stationarity behavior in our series, we use two different tests to avoid any spurious results; the first test is the Ng–Perron test (Ng & Perron 2001), which gives us four different tests, and the second test is the Dickey–Fuller test with bootstrapping critical values depending on the Park (2003) procedure (Table 2). Accordingly, the results show strongly that all our variables are not stationary at their levels, as the calculated statistic was less than the critical value, so we cannot reject the null hypothesis (H₀: Data series is not stationary). Still, they are stationary at their first differences, as the calculated statistic was greater than critical value, so the null hypothesis was rejected. It also confirmed that the ARIMA integrated will model one or more. This result means that we should use ARIMA models in our forecasting. After the stationarity test, based on the objective, we tried to build ARIMA and ETS models for forecasting.

After determining the stationarity of the data series, we tried to build a statistical model. The results inspired by Table 3 revealed clearly that ETS (Holt's nonlinear) models are the best in all the cases according to all the statistics (AIC, RMSE, MAE, MAPE, Theil, BP, VP and CVP) where we have the smallest values; this result means that the predictions using the ETS model will provide the minimum errors and deviations between forecasting and actual values (Mishra et al. 2022). However, we will forecast using the two methods for more accurate results.

After building the model, we tried to estimate the forecasting nature obtained from the best model of the data series. Data in Figure 1 clearly show significant deviations between forecast and actual values using confidence intervals and graphs, except for the PIN series (ARIMA model). These results revealed that the ARIMA forecasting is inappropriate for our series and that it is necessary to depend on ETS forecasting. The prediction obtained from the ARIMA model is that apple, banana, grape, pineapple, and mango production may reach up to 2790.8844, 35096.4473, 3465.5880, 1997.4287, and 44378.1405 thousand tons, respectively, in 2027.

The data in Table 4 below shows the parameters of the ETS technique. Hence, the smoothing factor α indicated a random walk in the PIN and MAN series since the α index is close to 1, in contrast to APP, BAN, and GRA series. Furthermore, it is clear from the β index that the values are equal to zero in most cases, meaning there are no trend changes.

Moreover, the damping index ϕ is held in APP, BAN, and GRA series with values close to one, decreasing the effect of the trend over time. The forecasting obtained from the ETS model indicated that apple, banana, grape, pineapple, and mango production may reach up to 4741.6803, 45968.0357, 5574.1481, 2044.5578, and 33998.0305 thousand tons in 2027 (Fig. 2).

The last step in this study is to provide a forecast evaluation using three different tests: Diebold–Mariano (2002), Giacomini–White (2006), and Clark–West (2007), to examine the forecasting accuracy between the two techniques in our study. Are they the same or not? As mentioned before, the ETS results are the best in our study. This result is confirmed using the three tests in Table 5, where we reject the null hypothesis of similarity between the two forecasts, especially the Clark–West test.