Modelling of Jowar production in Karnataka state - A hybrid approach

Introduction

Sorghum, (Sorghum bicolor), cereal grain plant of the grass family (Poaceae) and its edible starchy seeds. The plant likely originated in Africa, where it is a major food crop and has numerous varieties, including grain Sorghums used for food and grass Sorghums grown for hay, fodder and broomcorn, used in making brooms and brushes. In India Sorghum is known as Jowar, Cholam or Jola, in West Africa as Guinea corn and in China as Kaoliang. Sorghum is especially valued in hot and arid regions for its resistance to drought and heat. Globally, Sorghum produced about 626.9 lakh tonnes during 2024-25. United States stands top with 87.30 lakh tons (14%) followed by Nigeria with 64.21 lakh tons (10%), Brazil with 61 lakh tons (10%) and India with 50 lakh tons (8%) (fas.usda.gov, 2024-25) India ranks fourth in total Sorghum production with 50 lakh tonnes grown in an area of 39.95 lakh hectares in 2024-25, where majority of Sorghum was produced during rabi season. The kharif Sorghum with 17.37 lakh tonnes was grown predominantly in Rajasthan, Uttar Pradesh, Haryana and Madhya Pradesh, while rabi Sorghum with 31.87 lakh tonnes grown predominantly in Maharashtra, Karnataka, Tamil Nadu and Andhra Pradesh (indiastat, 2024-25). Karnataka is the eighth largest state in India, with a total size of 192 lakh square kilometres, or 6.3 per cent of the nation's total land area. It is estimated that 83.63 lakh ha of area under total food crops will be cultivated with the production of 146.31 lakh tonnes. Jowar is mainly grown in Kalaburagi, Belagavi, Vijayapura, Bagalkot, Raichur, Ballari, Yadgir, Bidar and Gadag districts and it is grown in an area of 6.48 lakh ha, with an annual production of 8.15 lakh tonnes and productivity of 1258 kg per hectare in Karnataka state (indiastat,2024-25). Jowar production in Karnataka has changed significantly and declining during the last couple of years, indicating shrinking cultivation area, diversion to more remunerative crops, marketdriven changes and regulatory initiatives. Nonetheless, nonlinear trends impacted by exogenous shocks and environmental conditions are also evident in the year-to-year fluctuation. The necessity for precise forecasting models that go beyond conventional linear techniques is highlighted by this complexity. In order to anticipate its production, Jowar has been chosen as one of the main crops for this study, which investigates sophisticated hybrid forecasting models such as ARIMA, TDNN, NLSVR, ARIMA-TDNN, and ARIMA-NLSVR. For those concerned in agricultural planning, policy creation, and market regulation, these models seek to offer trustworthy insights.

Materials and methodologies 2.1 Data Source

The study has been illustrated with the time series data. The secondary data on the Jowar production in Karnataka state for 63 years (from 1962-63 to 2024-25) was obtained from the Directorate of Economics and Statistics, Dept. of Agriculture, Cooperation and Farmers Welfare, Ministry of Agriculture, Govt. of India.

Statistical Methodologies

In this study, Auto Regressive Integrated Moving Average (ARIMA), Time Delay Neural Network (TDNN), Nonlinear Support Vector Regression (NLSVR) and a proposed hybrid methodology combining ARIMA-TDNN or NLSVR are employed to forecast the Production of Jowar in Karnataka state, India. The analysis is carried out using R software. To examine the stationarity of the time series data, the Augmented Dickey-Fuller (ADF) test is performed. The Brock-Dechert-Scheinkman (BDS) test is used to detect nonlinear dependencies in the residuals, ensuring an appropriate model selection process. Additionally, the Box Pierce/Ljung-Box test is applied to check for autocorrelation in the residuals, validating the adequacy of the fitted models. The forecasting performance of the models is assessed using Root Mean Square Error (RMSE), Mean Absolute Error (MAE) and Mean Absolute Percentage Error which provide measures of predictive accuracy and reliability.

Box-Jenkins ARIMA model

Before ARIMA, the concept of validation of data and Stationery was discussed briefly.

Validation of data

For model development and evaluation, the dataset comprising 63 annual observations was divided into training and testing sets. Out of these, 57 data points (1962-63 to 2018-19) were utilized for training the models, while the remaining 6 data points (2019-20 to 2024-25) were reserved as a testing set to assess the predictive accuracy and generalization performance of the models. This split ensures that the models are trained on a sufficiently large sample while also being validated on unseen data to avoid overfitting and to obtain realistic forecasting performance.

Stationary

A stochastic process is said to be stationary if its mean and variance are constant over time and the value of the covariance between the two time periods depends only on the disturbance or lag between the two time periods and not on the actual time at which the covariance is computed. It can be identified through Autocorrelation Function (ACF) of actual data, if the ACF does not die out rapidly, it indicates that the data is non-stationary. Under this situation, the auto correlation corresponding to most of lags are statistically significant. For reducing the data to stationary, selected data has to transform by taking first order differences (d=1). If the autocorrelation functions of first differenced data indicate a rapid decrease, then it can be concluded that the transformed data is stationary. If not, again the data has to be transformed by taking second order differences (d=2). Continuing in a similar way as that of d=1, thus the order of differencing i.e.,'d' can be determined.

Test of Stationary:

The order of differencing can be assessed by performing unit root test. The most widely used test for detecting the unit root (non-stationary) of time series is Augmented Dickey Fuller (ADF) test. In conducting Dickey Fuller (DF) test, it has assumed that the error term is uncorrelated. But in most of cases correlated, hence Dickey-Fuller (1979) has developed another test, known popularly known as the ADF test. This test is conducted by augmenting the regression (𝛥𝑌 𝑡 = 𝛽 1 + 𝛿𝑌 𝑡-1 + 𝑒 𝑡 ), by adding the lagged values of dependent variable(𝛥𝑌 𝑡 ). The ADF test consist of estimating following regression,

Where;  m is number of lagged difference terms required so that the error term  et is serially independent.

The test can be carried out by performing a t (= τ) -statistic of ′𝛿 ′ with dickey fuller table values at 5% and 10% LOS. If absolute t (= τ) statistic of ′𝛿 ′ is greater than its table value, it indicates that the series is Stationary.

ARIMA (Auto Regressive Integrated Moving Average)

It's denoted by ARIMA (p,d,q) which is a combination of Auto Regressive (AR) and Moving Average (MA) with an order of integration or differencing (d), where p and q are the order of autocorrelation and moving average respectively (Box and Jenkins, 2015) [2] . The Auto-regressive model of order p denoted by AR (p) is as follows:

The Moving Average (MA) model of order q or MA (q) can be written as:

Where;  c is constant term.  Øp is the p th autoregressive parameter.  θq is the q th moving average parameter and  et is the error term at time 't'.

ARIMA in general form is as follows:

Where Δ denotes difference operator like

Here, Yt-1, …, Yt-p are values of past series with lag 1, …, p, respectively.

Time Delay Neural Network (TDNN)

The ANN for time series analysis is termed as Time Delay Neural Network (TDNN). The time series phenomenon can be mathematically modelled using neural network with implicit functional representation of time, whereas static neural network like multilayer perceptron is presented with dynamic properties (Haykin, 1999). One simple way of building artificial neural network for time series is the use of time delay also called as time lags. These time lags can be considered in the input layer of the ANN. The TDNN is the class of such architecture. In this study, the R package "nnfor" was customized and implemented for developing Time Delay Neural Network (TDNN) models. This package provides a flexible and powerful framework for fitting artificial neural networks specifically designed for time series forecasting problems. In view of the residual nonlinearity detected through the BDS test, Time Delay Neural Network (TDNN) models were employed to model the residual series obtained from the ARIMA models for Jowar production in Karnataka. The TDNN models are well suited for capturing complex, nonlinear and dynamic dependencies in the time series data through their feedback and delay structure.

For each data series, multiple TDNN networks were trained by varying the number of neurons in the hidden layer from 1 to 10, with the sigmoid function employed as the activation function in the hidden layer, while keeping the output layers fixed and input layer varying from 2 to 5. To identify the optimal network configuration based on performance metrics such as RMSE, MAE, MAPE and model adequacy through the Ljung Box test p values. The structure of each network is denoted in the format Input-Hidden-Output. Following is the general expression for the final output Yt of a multi-layer feed forward time delay neural network.

The TDNN can be mathematically expressed as: [11] and it has been extended to the nonlinear regression estimation problems. Modeling of such problems is called as Nonlinear Support Vector Regression (NLSVR) model. The basic principle involved in Non-Linear Support Vector Regression (NLSVR) is to transform the original input time series into a higherdimensional feature space using a nonlinear kernel function, where a linear regression function is then fitted, enabling the model to capture complex, nonlinear relationships present in the data. The general form of NLSVR is:

Where;  𝛼𝑖 are the coefficients,  𝐾 (𝑥𝑖, 𝑥) is the kernel function and  𝑏 is the bias term.

In this study R package, "e1071" was customized and implemented. This package provides several options for SVR model fitting using different types of kernels like linear, radial basis, sigmoid and polynomial. Although the most often used and recommended function is radial basis. After detecting potential nonlinear patterns in the ARIMA residuals using the BDS test, a nonlinear model such as Nonlinear Support Vector Regression was applied to capture the nonlinear structures not modelled by ARIMA. The NLSVR models were fitted to the residuals of the optimal ARIMA models for Jowar production dataset. The radial basis function (RBF) kernel was chosen due to its flexibility and ability to model complex nonlinear relationships. The optimal hyperparameters for the NLSVR models, include Cost, Gamma, Epsilon were determined through iterative tuning to minimize the error while ensuring good generalization ability. A 10-fold cross validation procedure was adopted to ensure the robustness and generalization capability of the selected hyperparameters. This method involves dividing the datasets into ten subsets, iteratively training the model on nine subsets and validating it on the remaining one thereby minimizing overfitting and improving model reliability. The number of support vectors selected in each model reflects the complexity of the fitted regression function. The best performance value represents the minimal error achieved during the cross-validation process for each model. 1. a high dimensional feature space and then build the regression model in a new feature space. 2. Let us consider a vector of data set where is the input vector, yi is the scalar output and N is the size of data set. The general equation of Nonlinear Support Vector Regression estimation function is given as follows:

where, ϕ(.): R n → R nh is a nonlinear mapping function which maps the nonlinear mapping function which maps the original input space into a higher dimensional feature space vector. W ∈ R nh is weight vector, b is bias term and superscript T denotes the transpose.

Hyperparameter Tuning

The performance of NLSVR depends significantly on the choice of three key hyperparameters:  Cost (C): Controls the trade-off between model complexity and the degree to which deviations larger than ε are tolerated.  Gamma (γ): Defines the influence of each data point on the regression function.  Epsilon (ε): Specifies the margin of tolerance within which no penalty is given for prediction errors.

The Proposed Hybrid Methodology

The hybrid method considers the time series 𝑦 𝑡 as a combination of both linear and non-linear components. This approach follows the Zhang's (2003) [12] hybrid approach. Accordingly, the relationship between linear and nonlinear components can be expressed as follows

Where, 𝐿 ̂𝑡 and 𝑁 ̂𝑡 represent the linear and nonlinear component, respectively. In this work, the linear part is modeled using ARIMA model and non-liner part by TDNN or NLSVR. The methodology consists of three steps. Firstly, an ARIMA model is employed to fit the linear component.

Let the prediction series provided by ARIMA model be denoted as . In the second step, the residuals ( ) obtained from ARIMA model are tested for non-linearity by using BDS test (Brock et al., 1996) Once the residuals confirm the non-linearity, then, they are modelled and predicted using TDNN and NLSVR. Finally, the forecasted linear and nonlinear components are combined to generate aggregate forecast.

Where 𝐿 ̂𝑡 and 𝑁 ̂𝑡 represent the predicted linear and nonlinear component, respectively.

Performance Criteria

The model's performance was evaluated using four key metrics: Coefficient of Determination R 2 , which measures goodness-of-fit; Root Mean Square Error (RMSE), Mean Absolute Error (MAE), which quantify prediction errors in absolute terms and Mean Absolute Percentage Error (MAPE), which expresses errors as a percentage. Furthermore, the BDS test was used to check for non-linear structures in the data, while the Augmented Dickey Fuller (ADF) test confirmed if the time series data was statistically stationary or not.

Results and discussions 3.1 Descriptive Statistics of Jowar Production of Karnataka State

A detailed descriptive statistical analysis was carried out to examine the characteristics and distribution of the production (in thousand tonnes) of Jowar in Karnataka for the period from 1962-23 to 2024-25. These measures provide information on central tendency, dispersion, distribution pattern and behaviour of the data. The summary of descriptive statistics is presented in Table 3.1. The production data set been plotted vs time as seen in Figure 3.1 3.2, since the null hypothesis was rejected at the five percent level of significance because the p value was less than 0.05. Consequently, the ACF plots of the differenced series for production was examined in Figure 3.2. The differenced series exhibited a rapid decay of autocorrelations with most of the autocorrelation coefficients falling within the 95 per cent confidence limits, indicating that the data had attained stationarity.

After setting d=1 to achieve stationarity, several ARIMA models were developed. The models were compared using R 2 , RMSE, MAE and MAPE criteria which was shown in Table 3.3. Among them, ARIMA (3,1,1) was identified as best model with high R 2 , low RMSE, MAE and MAPE values. The ARIMA (3,1,1) fitted equation:

Further, a residual analysis was conducted to assess the adequacy of the selected model. The analysis revealed that none of the lag values in the residual ACF and PACF plots were significant, as shown in the Figure 3.2. Additionally, Ljung-Box test and Kolmogorov-Smirnov test were conducted to check whether the residuals were independently and normally distributed or not. The Ljung-Box test statistic was Q = 0.008 (p = 0.92) K-S test statistic was D = 0.06 (p value = 0.96). The p values were greater than 0.05 indicating that the residuals were independently and normally distributed. This confirms the validity and goodness of fit of the selected ARIMA models for forecasting Jowar production in the Karnataka State.

Fitting of Time Delay Neural Network (TDNN) for Jowar Production in Karnataka State

From Table 3.4, it is shown that the best performing network was 4-4-1, which achieved low RMSE, MAE and MAPE values. Further, residual analysis was also carried out to check the adequacy of the selected neural network model and it was discovered that none of the lags of residual ACF chart were found to be significant as per Figure 3.3. Furthermore, the p-value (0.92) for the Ljung-Box test was greater than 0.05, indicating that the residuals were independently distributed. The residuals of the fitted ARIMA model for Jowar production were tested for nonlinearity using the Brock-Dechert-Scheinkman (BDS) test. As seen in Table 4.6, a significant p-value was found at Embedding Dimensions 2 and 3, suggesting the existence of a nonlinear structure not captured by the fitted ARIMA model. These findings showed that while certain nonlinear patterns persisted, ARIMA models successfully captured the linear structure in the data series. Therefore, using a nonlinear model to represent the nonlinear structures was justified.

Table 3.6. BDS Test on ARIMA (3,1,1) Residuals

In the present study, to enhance the forecasting accuracy nonlinear models are integrated with the ARIMA structure. Further the ARIMA residuals are fitted by nonlinear models such as Time Delay Neural Network (TDNN) and Nonlinear Support Vector Regression (NLSVR). This hybrid ARIMA-TDNN and ARIMA-NLSVR approach combines the strengths of ARIMA for capturing linear trends and nonlinear models for addressing remaining nonlinear structures in the data.

Fitting of Time Delay Neural Network (TDNN)

Models on ARIMA Residuals Table 3.7 demonstrates that 4-7-1 was the best-performing network, achieving low RMSE, MAE, and MAPE values. In order to assess the suitability of the chosen neural network model, residual analysis was also performed. According to Figure 3.4, none of the delays in the residual ACF chart were found to be significant. Additionally, the Ljung-Box test's p-value (0.93) was higher than 0.05, suggesting that the residuals had an independent distribution. These results imply that the chosen model provides a good match to the ARIMA residuals.

Acknowledgement

International Journal of Agriculture and Food Science https://www.agriculturaljournals.com

References

1. Ali AM, Mushi J. Investigation of a suitable hybrid time series model for predicting clove rice. J ICT Syst. 2023;1(2):1-16.
2. Box GE, Jenkins GM, Reinsel GC, Ljung GM. Time series analysis: Forecasting and control. 5th ed. Hoboken (NJ): John Wiley & Sons; 2015.
3. Devra SJ, Patel DV, Shitap MS, Raj SR. Time series forecasting of price for oilseed crops by combining ARIMA and ANN. Int J Stat Appl Math. 2023;8(4):4054.
4. Kumar S, Arora K, Singh P, Gupta AK, Sharma I, Vatta K. Performance comparison of ARIMA and time delay neural network for forecasting of potato price in India. Agric Econ Res Rev. 2022;35(2):1-12.
5. Kumar TM, Prajneshu. Development of hybrid models for forecasting time-series data using nonlinear SVR enhanced by PSO. J Stat Theory Pract. 2015;9(4):699711.
6. Jha GK, Sinha K. Time-delay neural networks for time series prediction: An application to the monthly wholesale price of oilseeds in India. Neural Comput Appl. 2014;24(3):563-571.
7. Rathod S, Mishra GC. Statistical models for forecasting mango and banana yield of Karnataka, India. J Agric Sci Technol. 2018;20(4):803-816.
8. Rathod S, Singh KN, Patil SG, Naik RH, Ray M, Meena VS. Modelling and forecasting of oilseed production of India through artificial intelligence techniques. Indian J Agric Sci. 2018;88(1):22-27.
9. Ray M, Rai A, Ramasubramanian V, Singh KN. ARIMA-WNN hybrid model for forecasting wheat yield time-series data. J Indian Soc Agric Stat. 2016;70(1):63-70.
10. Raghu KN. Statistical analysis of jowar production in Karnataka [doctoral dissertation]. Jabalpur: Jawaharlal Nehru Krishi Vishwa Vidyalaya; 2020.
11. Vapnik V, Golowich S, Smola A. Support vector method for function approximation, regression estimation and signal processing. In: Mozer M, Jordan M, Petsche T, editors. Advances in Neural Information Processing Systems. Cambridge (MA): MIT Press; 1997. p. 281-287.
12. Zhang GP. Time series forecasting using a hybrid ARIMA and neural network model. Neurocomputing. 2003;50:159-175.