Modeling and Forecasting of Milk Production in the Western Zone of Tamil Nadu

Milk is closely linked to our civilization, as the domestication of dairy animals started thousands of years ago. It is now well connected to several sectors like food, health, pharmaceuticals, cosmetics, etc. In most countries, milk ranks among the top five agricultural commodities in terms of monetary returns. Besides their dietary purpose, it plays a vital role in the growth of the Indian economy employment generation, nutritional security and protection of native breeds. Among agricultural products, milk and other dairy products account for 14% of global trade (Mottet et al., 2018). The estimated global milk production was about 906 million tonnes (MT) in 2020 (an increase of 2.0% from 2019). At the beginning of India’s independence in 1950, milk production was around 17 MT annually. However, the launch of Operation Flood in 1970 resulted in a rise in annual production of 23.2 MT in 1973-74 to 209.96 MT in 2020-21, an eight-fold growth in less than five decades (Khera et al., 2022). Now, India is the world’s leading milk producer, accounting for 23% globally and 52% of Asia’s total output. As the single-largest agricultural commodity in India, the dairy industry contributes around 26 % to the entire agriculture GDP and more than 5% to the national economy (Chatellier, 2021). The per capita availability of milk in India has also increased to about 427 grams per day (2021) against the ICMR recommendation of 300 grams per day. This shows the sustained growth in milk availability for the burgeoning subcontinent’s population.
       
The United Nations reported that the global population is expected to reach 9 billion and India will have a population of 1.66 billion by 2050. India should produce around 401.4 million metric tonnes of milk to meet the demand in 2050 (Aatralarasi et al., 2021). The emerging population growth is expected to elevate the need for all agricultural products, especially animal protein. The increase in the production of resources should majorly come from intensification rather than expansion. The increase in milk production has to come from productivity improvements. Still, the number of interactive factors related to breeding, feeding, health and management continue to yield less milk. At least in the short run, there is a scope to raise milk yield with the existing milch stock and its quality (Patel et al., 2022). However, all these studies were confined to a particular zone or region, state or production system, which widely lacks policy implications. To control market demands, stakeholders and authorities had to formulate some effective short- and long-term plans. These plans aid in managing milk supply using precise estimations for milk production derived from well-accepted prediction methods. Hence there is a need to forecast the milk production to foresight if any demand would create in the future.
               
Forecasting is an important aspect in a developing economy like India. It aids in decision-making and planning for sustainable growth, poverty alleviation and overall development (Zhang, 2003). Different forecasting techniques are performed using statistical analysis. These statistical predicting models are used to develop suitable time series models by identifying the trends and patterns of past data. This research objects to model and forecast the milk production of the western zone of Tamil Nadu using ARIMA and ANN techniques.

Data collection
 
The study aims to forecast the milk production of the western zone of Tamil Nadu, India. Yearly data on milk production from 1976 to 2020 (National Dairy Development Board, Anand, Gujarat) was used for the forecasting purpose. To fit the model, the data was split into training and testing data i.e., data from 1976 -2015 was taken for model building and data from 2016-2020 was accepted for model validation purposes. R software was used to fit and forecast the study data.
 
Compound growth rate (CGR)
 
Exponential function was to analyse the compound growth rate, which was initially transformed into a linear form by taking the log function (Kalidas et al., 2020). The final equation for estimating the growth rate was given by:
   
Where,
r = Compound growth rate.
b = Regression coefficient.
 
Time series models
 
Forecasting is a technique where future observations are predicted using past observations. It is made possible using time series models, extending its application to all fields (Hamid et al., 2016). Generally, no time series models are assumed to be the best model for data, as the data don’t follow similar patterns in all cases. So, the suitable method for data is determined based on its nature. Among several time series models, some were tried fitting for the data. The best-fitted model will be used for forecasting future milk production (Demir and Kirisci, 2022). Different time series models employed in the study are detailed below.
 
Auto-regressive integrated moving average model
 
Auto-regressive integrated moving average (ARIMA) is a statistical time series model used to estimate and forecast the time series data (Punyapornwithaya et al., 2022). It is a linear model used for handling univariate data. The model comprises three major parts. The first part of the model is auto-regressive (p), the second part is integration (d) and the third part is moving average (q). The ARIMA model can be expressed as follows:
   
Autoregressive operator: φ (B) = 1 - φ₁ (B) - φ₂ (B)² -...- φ_p (B)^p
Moving average operator: q(B) = 1 - θ₁(B) - θ₂(B)² -... - θ_p(B)^q                  
Where,
ɛ_t = White noise or error term.
d = Differencing term.
B = Backshift operator, i.e.        
Being a linear model, it can be applied only to the data unless it is stationary. If the data is non-stationary, it should be converted into stationarity by differencing it to the respective order (Deshmukh and Paramasivam, 2016). The major steps involved in building the ARIMA model.
 
Identification
 
The values of q and p can be found using the plots of ACF (Auto Correlation Function) plot and PACF (Partial Auto Correlation Function) plot. The p and q denote the model information from past values and errors. The d in Integration indicates the number of times the data should integrate to convert into stationary.
 
Parameter estimation
 
The parameters of ARIMA models are estimated using Maximum Likelihood Estimation (MLE) method by which the values of AR and MA values are found.
 
Diagnostics
 
The model which gives the low AIC and BIC values is chosen to be the best fitted model for the data (Naveena and Subedar, 2017). Finally, the data is forecasted if the residual of fitted model is uncorrelated, i.e., white noise.
 
Artificial neural network model
 
An artificial Neural Network (ANN) is a computational model inspired by the structure and functionality of biological neural networks. They are the network of interconnected neurons mimicking the function of the human brain. A feed-forward neural network is one of the basic neural networks and serves as a non-linear time series model for forecasting purposes (Rathod et al., 2017). They are made of input, hidden and output nodes. Every unit in a layer is related to every unit in the previous layer. The data is given through the input node and the result is obtained from the output node. The in-between hidden layer is the place where processing is done. Each layer consists of weights and biases. The numbers of input and hidden nodes are determined by experimentation, as there is no theoretical base for finding these parameters (Mishra et al., 2021). The mathematical representation of the ANN model is given as:
 
Where,
j (j = 0,1,2,...,q), ij (i = 0,1,2,..., p and j = 1,2,...,q) = Weights.  
0_j = Bias terms. 
ɛ_t = White noise.
 
Estimation of model performance
 
Model validation helps estimate the performance of different models by comparing the actual and fitted values. The various error measures used for model selection are Root Mean Square Error (RMSE), Mean Absolute Error (MAE) and Mean Absolute Percentage Error (MAPE). The model with a low error rate is the optimum model for the study data.
  Y_t = Actual values.
Y_t = Predicted values.
n = Number of observations.

Descriptive statistics is calculated for the study, which is given in Table 1. The maximum milk production was recorded in 2020 (833.12 MT), whereas the minimum milk production was recorded in 1982 (165.20 MT). The mean milk production for the last 45 years is 423.49 MT, with a standard deviation of 175.30. The kurtosis is -0.04 with positively skewed. The coefficient of variation is 41.39%.
 

       
The compound growth rate was calculated to inspect the growth rate percentage in every prefixed period (Fig 1). This study estimated the growth rate for every five years of milk production using an exponential function (Anjum, 2018). The highest growth rate was recorded during 1991-95 (6.31%), followed by 2011-15 (6.20%), whereas the lowest growth rate was recorded during 2001-05. During 1976-80 a negative growth rate (-1.83%) was observed and an overall growth rate of 2.71% was registered over the last 45 years.
 

       
This study used two different time series models for model building and model validation. The model that fits better for the data with a low error rate was used to predict future milk production. Both the models were fitted to the data initially split into training and testing datasets. The Autoregressive Integrated Moving Average (ARIMA) is a linear model, fitted only when the datasets are stationary. If the datasets are not stationary, they are converted into stationary using differencing (Kour et al., 2017). Augmented Dickey-Fuller (ADF) test and Philips-perron test were two autocorrelation tests used for testing stationarity in the data. Table 2 revealed that the p-value of the two tests was greater than 0.05, i.e., they are not stationary. The graph (Fig 2) shows that the spikes came out of significant lines in both ACF and PACF plots. So, the data was differenced one time to convert them into stationarity.
 

 

       
Table 3 revealed that the ARIMA model (1, 1, 2) was fitted to the training datasets with the lowest AIC and BIC values, i.e., 416.18 and 422.84, respectively. The test statistics of AR1 was -0.16 with a standard error 1.10, MA1 was 0.12 with a standard error of 0.07 and MA2 was 0.08 with a standard error of 0.21. In addition, the box-pierce test and Shapiro-wilk test were conducted to test the autocorrelation and normality of the fitted model (Pardhi et al., 2018). The p-value of both tests was greater than 0.05, meaning that the residuals are white noise and normal.
       

 
Feed-forward neural network (FFNN) is a type of ANN used to fit the data in this study (Table 4). The initial step is to find the network of the ANN model. Different combination input nodes (1:10) and hidden nodes (1:10) were tried for the training dataset and the appropriate node was found after several iterations (Tealab et al., 2017). The output node is always one. Finally, the network 6-2-1 was found to be the best fit for the data with weights 17. Linear sigmoidal activation function was used in model building. The Box-pierce test and Shapiro-wilk test show that the fitted ANN model’s residuals are not autocorrelation and are normal.
       

 
Since both the fitted models are checked for white noise and normality tests at the model-building phase, they can be taken for model validation. Table 5 indicates the predictive performance of ARIMA and ANN models. Error measures like RMSE, MAPE and MAE were used for comparing the model (Sankar and Prabakaran, 2012). The ARIMA model was found to have a low error in training and testing, followed by ANN with minor differences. Even though the error rate is calculated for both training and testing datasets, the testing datasets, i.e., model evaluation, are deciding factor for the model section (Aslanargun et al., 2007). Based on that criterion, the ARIMA model’s error rate is comparatively lower than ANN in all aspects. Therefore, the ARIMA model is selected as the best fitted model for study and used for forecasting purposes.
 
 

Since ARIMA models fit well for the study data with a low error rate, satisfying all necessary conditions of a time series model (Fig 3), it can be used for forecasting the future milk production (MT) (Yonar et al., 2022). Table 6 indicates that the forecasted values of milk production from 2021 to 2025. The CGR for the forecasted milk production is 0.02%. Fig 4  gives the graphs showing actual and fitted values, along with predicted values obtained using the ARIMA model.
 

 

 

The Indian dairy sector covers a sizeable portion of the world’s dairy resources. The livestock industry supports the national economy and the nation’s socio-economic development. The overall growth rate for milk production obtained by CAGR for 45 years  is 2.10%. Further statistical modelling procedures employed for milk production reveal that the selection of a suitable time series model will always depend on the nature of the data. Likely ARIMA model is selected as the best model in spite of ANN even if it is considered the powerful model. The CAGR for forecasted milk production from 2020-2025 was 0.02%. When considering the approximate projected population growth from various studies, the forecasted growth rate could be insufficient to meet our demand. Thus, we need more and more research in the dairy sector to increase milk production.

None.