Temporal Trends and Future Projections: A Deep Dive into India’s Buffalo Milk Production Through Time Series Modelling

India is the world’s top milk producer, accounting for 24 percent of the world’s milk production in 2021-2022. Over the course of the last eight years, or between 2013-14 and 2021-22, India’s milk output has climbed by 61%, reaching 221.1 million tonnes in 2021-22. Compared to the 2020-21 financial year, there has been a 5.29% increase in milk production. Rajasthan (15.05%), Uttar Pradesh (14.93%), Madhya Pradesh (8.6%), Gujarat (7.56%) and Andhra Pradesh (6.97%) are the top 5 states that produce milk.  Collectively, they account for 53.11% of the nation’s milk production. In 2022-2023 India exported 67,572.99 MT of dairy products to the world, valued at $284.65 million.
       
India is recognized as a rich source of buffalo germplasm resources, hosting all acknowledged high-producing breeds of this species. As the Indian dairy industry undergoes transformative changes, the buffalo sector plays a crucial role in overall milk production, contributing 66.3% of the total global volume of buffalo milk. With Indian buffaloes accounting for nearly 56% of the world’s buffalo population, about 25% of global buffalo meat is produced in India. Buffalo farming has evolved into a livelihood and resource-generating enterprise for diverse segments of our farmers. It plays a major role in poverty alleviation and commercial buffalo enterprises now provide employment opportunities to rural communities.
       
Buffalo farming holds significant importance in our rural economy, contributing not only to milk production but also being valued for meat and draught purposes. Buffalo milk, known for its higher fat content of 7-7.5%, nearly double that of cows, commands a higher market price. Milk obtained from the mammary glands of female buffaloes, known as buffalo milk, is a highly nutritious dairy item that holds significant cultural and economic importance. Differing in composition, taste, and nutritional value from cow’s milk, buffalo milk boasts higher levels of fat, protein, and minerals. Recognized for its creamy texture and robust flavour, it typically contains 7% to 8% fat, making it ideal for producing butter, ghee and other dairy products-significantly higher than cow’s milk. Rich in essential nutrients like protein, calcium, phosphorus, vitamin A and B-complex vitamins, buffalo milk’s elevated fat content also contributes to its energy density (Thakore and Jain, 2018).
       
Buffalo milk extensively utilized in creating various dairy products such as paneer, yogurt and traditional sweets, buffalo milk imparts a luxurious and creamy quality to these items. In regions where buffalo farming is prevalent, buffalo milk plays a pivotal role in the dairy industry, significantly contributing to overall milk production and serving as a livelihood source for many farmers. Buffalo milk indeed often contains a higher fat content compared to cow milk. This higher fat content makes buffalo milk desirable for various dairy products. The pricing of milk is commonly influenced by its fat content, among other factors. In many cases, higher fat content can lead to higher profit margins for dairy producers, as products derived from such milk, like butter and cream, are often priced at a premium. In some regions, buffalo milk plays a significant role in the dairy industry. The fact that buffalo milk accounts for the largest share (55%) of the total milk produced in the country indicates its importance in the overall dairy market. This could be due to factors such as the adaptability of buffaloes to different climates, their ability to thrive on diverse feed sources, or specific preferences for buffalo milk and its products in the local culture (Singh, 2015). Notably, countries like India are major global contributors to buffalo milk production, with the Indian subcontinent home to high-yielding buffalo breeds that substantially augment the world’s buffalo milk supply. To summarize, buffalo milk stands as a valuable dairy product with distinct attributes, making it versatile for a range of culinary applications. Its higher fat content, nutritional richness and economic significance underpin its widespread consumption and utilization across different cultures worldwide.
       
This study explores the temporal patterns of buffalo milk production in India using various forecasting models, including ARIMA, SES, Seasonal Naive and ETS. The aim is to predict future production trends and understand the underlying mechanisms shaping buffalo milk production dynamics. The study delves into the economic and cultural significance of buffalo farming, offering nuanced insights into the multifaceted landscape of India’s dairy industry. By comparing and contrasting the performance of different forecasting models, the study aims to identify the most effective approaches for modelling buffalo milk production, providing actionable insights for policymakers and stakeholders to optimize production strategies and contribute to the sustainable growth of India’s dairy sector.

Data collection and preparation
 
The study utilized an extensive time series dataset comprising annual measurements of raw buffalo milk production from 1961 to 2021. The dataset was meticulously curated from authoritative agricultural records and sources, ensuring its integrity and reliability. Preprocessing steps were undertaken to address missing values, ensure chronological order and aggregate the annual data into consistent time intervals.
 
Forecasting models and equations
 
Various time series models are available in literature which is to be used based on the characteristics of the data. The ability of different models to forecast the time series values is assessed by using forecast evaluation measures. Four distinct forecasting models were considered and applied to capture the temporal patterns in the milk production time series data.
 
Auto regressive integrated moving average (ARIMA)
 
The Box-Jenkins ARIMA (p, d, q) model, introduced by Box and Jenkins Box et al., (2015), stands as a widely employed technique for constructing univariate time series forecasting models. Developed by George Box and Gwilym Jenkins in the 1970s, the ARIMA model provides a mathematical framework for predicting processes. The Box-Jenkins modelling approach encompasses the steps of identifying a suitable ARIMA process, fitting it to the available data, and subsequently employing the established model for predictive purposes. The ARIMA model represents the time series as a function of its own past values, integrating autoregressive (AR), differencing (I), and moving average (MA) components. The general ARIMA (p, d, q) equation is expressed as:
   
Where
Y_t = Milk production,
ε_t = Independently and normally distributed with zero mean and constant variance for t=1, 2, ….. , n; and Ø_p and θ_q are also estimated (Venu et al., 2023: Deshmukh Surendra Sagar and Paramasivam, 2016: Mahajan Sunali et al., 2020).
 
Exponential smoothing state space model (ETS)
 
The Exponential Smoothing State Space Model (ETS) is a forecasting method that combines exponential smoothing with a state space approach (Mphale and Narasimhan, 2022). ETS models are particularly useful for time series data that exhibit trends and seasonality. Developed to provide a flexible and adaptive framework for time series forecasting, ETS models have gained popularity in both academia and industry. ETS models provide a flexible framework that can handle a wide range of time series patterns, making them suitable for forecasting in diverse applications. The implementation of ETS models is available in various statistical software packages, and they can be adapted to incorporate exogenous variables for improved forecasting accuracy in ETSX models. The flexibility and adaptability of ETS make it a valuable tool for analysts and practitioners dealing with time series forecasting. The basic ETS (A, A, A) model (where “A” stands for additive) can be represented by the following equations:
 
Where
y_t= Observed value at time, t.
ŷ_t+h = Forecast for time t+h.
Where,
h = Forecasting horizon,
l_t = Level at time t.
b_t = Trend at time t.
s_t = Seasonal component at time t, α, β.
γ = Smoothing parameters for level, trend and seasonality, respectively.
s = Length of the seasonal cycle.
 
Simple exponential smoothing (SES)
 
Simple Exponential Smoothing (SES) is a basic and widely used method for time series forecasting (Ostertagova and Ostertag, 2012). It belongs to the family of exponential smoothing methods, which are designed to capture patterns and trends in time series data. SES is particularly useful when the data does not exhibit complex patterns or seasonality. Simple Exponential Smoothing is a good starting point for time series forecasting, especially when the data has a relatively stable trend and does not exhibit seasonality. It is easy to implement and computationally efficient, making it a practical choice for simple forecasting tasks. However, the appropriateness of SES depends on the specific characteristics of the time series being analysed. The SES equation is as follows:
   
Where
ŷ_t+1  = Forecast for the next time period (t+1).
y_t  = Actual observation at time t.
ŷ_t  = Forecast for time t.
α = Smoothing parameter.
       
Also known as the smoothing coefficient, and it takes values between 0 and 1. α determines the weight given to the most recent observation in the smoothing process.
 
Seasonal naive model
 
The seasonal naive model is a simple and intuitive time series forecasting method that relies on the observation that many time series exhibit seasonality, where patterns repeat over fixed intervals (Pala and Atici, 2019). The Seasonal naive model predicts future values based on the most recent observation from the same season in previous periods. While it may not capture all nuances in the data, it serves as a baseline for comparison and is particularly useful when the time series exhibits clear and consistent seasonality. If we denote the time series data as y_t.
Where,
t = Time index and the seasonal period as s, then the equation for the seasonal naive method is given by:
   
Where,
ŷ_t = Forecast for time t.
y_t-s = Observation from the same season in the previous year, as s is the seasonal period.
       
The Seasonal Naive model assumes the next value to be identical to the last observed value at the same seasonal point.
 
Model evaluation and fitness metrics
 
Each model was trained on the designated training set and utilized to generate forecasts for the test set. The performance of these models was assessed using common forecast accuracy metrics:
 
Mean error (ME)
 
Mean error is the average of all the errors in a set of observations. That is the arithmetic mean of the differences between predicted values and actual values. ME is given by:
   
n = Number of data points.
P_i  = Predicted value for the ith data point.
Q_i = Actual value of the ith data point.
       
If the mean error is close to zero, it suggests that the observations are accurate. That is, a non-zero ME indicates a bias in the observations. However, sometimes ME might not give the correct picture of accuracy, especially if there are both positive and negative errors that cancel out.
 
Mean absolute error (MAE)
 
Mean Absolute Error (MAE) (Willmott, 2005: Chai et al., 2014) is a metric used to evaluate the performance of a regression model. It measures the average absolute difference between the actual value and predicted values. The mean absolute error is calculated by:
   
Where,
n = Number of data points.
P_i= Predicted value for the i^th data point.
Q_i= Actual value of the i^th data point,
   
Root mean squared error (RMSE)
 
Root Mean Squared Error (Willmott and Matsuura, 2005; Chai and Draxler, 2014) is commonly used to evaluate the performance of a predictive model, particularly in regression analysis. It is the quadratic mean of the differences between observed values and actual values, giving more weight to larger errors. RMSE is given by:
   
n = Number of data points.
P_i = Predicted value for the i^th data point.
Q_i = Actual value of the i^th data point.
 
Mean absolute percentage error (MAPE)
 
MAPE (Moreno, 2013; Goodwin et al., 1999) is used to assess the accuracy of a forecasting or prediction model, especially in the context of time series analysis. MAPE is expressed as a percentage and measures the average absolute percentage difference between the predicted and actual values. If n is the number of samples.
P_i = Predicted value for the i^th data point.
Q_i= Actual value of the i^th data point, then MAPE is given by;
   
MAPE provides a clear and intuitive way to assess the accuracy of forecasting models, especially in situations where understanding errors in percentage terms is important. Lower MAPE values indicate better accuracy, with a value of 0 indicating a perfect match between predicted and observed values. However, MAPE is sensitive to zero values in the actual data, which can result in undefined or infinite percentages. Despite such limitations, MAPE is probably the most widely used forecasting accuracy measurement. These metrics were instrumental in evaluating the accuracy and reliability of the models in predicting future milk production values.
 
Residual analysis
 
A thorough residual analysis was performed on the ARIMA model to confirm its dependability and assumptions (Revathi et al., 2023: Jaiswal Priyanka and Bhattacharjee Mahua, 2022). The study included visual inspections, Q-Q plots, histograms and the Autocorrelation Function (ACF) to confirm the distribution's normality and look for any temporal trends or systematic errors. The residuals' absence of evident trends or correlations supported the model's future predictions of oilseed yield accuracy. This research confirmed that the independence, normality, and constant variance assumptions were met, which improved the model's prediction reliability.

Ljung-box test
 
The Ljung-Box test is a statistical test used in time series analysis to assess whether the residuals from a time series model exhibit significant autocorrelation at different lags Dare (2022). The test is particularly useful for evaluating the adequacy of the chosen model by examining whether there are any remaining patterns in the residuals that the model has not captured. The Ljung-Box test statistic is based on the sum of the squares of the sample autocorrelations of the residuals at different lags. It is calculated as follows:
   
Where,
Q = Test statistic, n is the sample size.
p̂k = Sample autocorrelation at lag k.
h = Maximum lag considered.

The study utilized an extensive time series dataset comprising annual measurements of raw buffalo milk production from 1961 to 2021 taken from FAOSTAT (Fig 1). The stationarity of buffalo milk production time series data was assessed using the Augmented Dickey-Fuller (ADF) test. The test fitted a regression model with lags and a time trend to the differenced series. The test-statistic value of 0.7954, which is very much greater than critical, indicates that the data is non-stationary. Some lagged differences had statistically significant regression coefficients, indicating the persistence of temporal patterns. Non-stationarity implications will be considered for modelling.
 

 
ARIMA
 
Using R software, the best ARIMA model for time series data was found using the forecast package. As we already found that the data is not stationery from ADF test differencing technique is applied. For d=2 the data become stationery. ARIMA (1, 2, 1) was found to be the best model with the lowest AIC value among the several ARIMA models considered. The model equation is as follows:
   
Where,
Y_t = Differenced and stationary time series,
ε_t = White noise error term at time t, and the autoregressive and moving average coefficients are given by 0.5423 and -0.9117, respectively.
       
The autoregressive coefficient of 0.5423 indicates a positive correlation between current buffalo milk production and previous values, implying that production levels will persist over time. In contrast, the moving average coefficient of -0.9117 indicates that the model rapidly adjusts to deviations from the predicted values, reflecting a strong tendency for fluctuations in buffalo milk production to revert to the long-term average. The Ljung-Box test was used to examine the residuals of the ARIMA (1, 2, 1) model for time series data on buffalo milk production. The test produced a high p-value of 0.6995, indicating that there was no significant evidence of autocorrelation at lags of up to ten. As a result, the ARIMA (1, 2, 1) model appears adequate for capturing and explaining temporal patterns in buffalo milk production time series data. The residual plots of ARIMA (1, 2, 1) is given in Fig 2, which also suggest the same conclusion as above.
 

 
Simple exponential smoothing
 
The Simple Exponential Smoothing (SES) model, when applied to the buffalo milk production data, reveals a notably high smoothing parameter (a = 0.9999\). This parameter signifies that the model heavily relies on recent observations to generate forecasts, with minimal weight attributed to older data points. In essence, the model prioritizes the most recent data, implying a strong preference for short-term trends over long-term patterns in buffalo milk production. The initial state (l = 11,809,985.2045\)) represents the starting level of the smoothed series, indicating that the SES model begins its forecast from this baseline value. The standard deviation of the residuals (s=1,988,051) provides a measure of the variability of the observed data around the smoothed values, indicating the level of uncertainty inherent in the model's predictions. Fig 3 represents the residual plots of the SES model. Evaluation metrics such as the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) suggest that the SES model provides a reasonably good fit to the data, with lower values indicating better model performance. However, despite the relatively low AIC and BIC values (2024.055 and 2030.388 respectively), the error measures reveal a Mean Absolute Percentage Error (MAPE) of 4.07%, indicating that, on average, the model's forecasts deviate by approximately 4.07% from the actual values. Moreover, the Ljung-Box test yields a very low p-value (2.397ε-12), indicating that the residuals from the SES model are not completely free of autocorrelation, suggesting potential room for further model refinement or exploration of alternative modelling techniques to enhance predictive accuracy.
 

 
Seasonal naïve model
 
A seasonal model is developed for the milk production data and tested for model fitness. The Ljung-Box test provided a p value (1.939e-10) significantly smaller than the conventional significant levels which implies high autocorrelation among residuals in the model. Fig 4 depicts the residual plots which in turn agree with the above conclusion. As a result, like SES model, Seasonal Naïve Model also fails to describe the temporal patterns in the data.
 

 
ETS model
 
ETS (A, Ad, N) model was developed for the data using R software and model equations are as follows:
 
The model exhibits significant forecasting implications due to its emphasis on recent observations and high a parameter (0.9999), indicating adaptability to short-term fluctuations, whereas the b parameter assumes a moderate value of 0.8058, indicating a balanced consideration of recent trends and historical patterns. The introduction of damping to the trend component, as denoted by a f value of 0.8843, suggests a controlled reduction in the amplitude of the trend over time. This nuanced approach to trend dynamics is particularly relevant in capturing both short-term variability and the broader trajectory of the time series. Fig  5 represents the residual plots of the ETS model. The fact that the ACF lags are clearly within the confidence interval lines indicates that there is no significant autocorrelation between the residuals. The Ljung-Box test produced a p value of 0.3443, which is greater than the significant levels and implies the same conclusion.
 

 
Model fitness metrics
 
Based on the residual analysis, we discovered that the ARIMA and ETS models explain the data's behaviour better than the SES and Seasonal naive models. We will now consider different model fitness matrices for these four models to validate these findings. Table 1 displays the ME, MAPE, MAE and RMSE values for all four models that we developed.
 

       
When the forecast accuracy of ARIMA, SES, Seasonal Naive and ETS models is compared, distinct patterns emerge. The ARIMA model performed well, with low Mean Absolute Percentage Error (MAPE) values on both the training (2.08%) and test (3.18%) sets, indicating its ability to capture underlying data patterns. While competitive, the SES model had higher MAPE values (4.12% on training and 12.26% on test), indicating potential limitations in handling more complex temporal dynamics. Seasonal Naive demonstrated competitive accuracy (4.07% on training and 8.64% on the test set), excelling in capturing seasonal patterns. The ETS model, which is like ARIMA, demonstrated strong performance (2.18% on the training set and 3.26% on the test set). Table 2 gives forecast of buffalo milk in India using SES, ARIMA, Seasonal Naïve and ETS models with 80% and 90% confidence interval and Fig 6 shows a diagrammatic representation.
 

 

This study investigates the dynamics of buffalo milk production in India over six decades, focusing on data from 1961 to 2021. The data showed non-stationarity, prompting the use of different modelling techniques. Four forecasting models were used: ARIMA, SES, Seasonal Naïve and ETS. The ARIMA model was found to be most suitable, capturing both short-term fluctuations and long-term trends. The SES model, based on recent observations, showed competitive performance but higher Mean Absolute Percentage Error (MAPE) values. Seasonal Naïve and ETS models failed to capture temporal patterns. The ARIMA and ETS models demonstrated superior forecast accuracy, with lower MAPE values indicating higher forecast accuracy. The forecasts, with 80% and 90% confidence intervals, provide valuable insights for policymakers and stakeholders in optimizing production strategies and fostering long-term growth in the dairy industry. The study emphasizes the importance of robust forecasting models in strategic decision-making and advancing the sustainability and resilience of the dairy sector. Further research and refinement of modelling techniques are needed to address the complexities and uncertainties in dairy production dynamics.

The authors declare that they have no conflict of interest or competing interests.