India is the world’s largest milk producer, contributing nearly 24% of global milk production (FAO, 2026). India is also one of the world’s largest consumers of fluid milk, with consumption projected to increase by approximately 2.5-3% to 90-92 million metric tons, driven by population growth and rising incomes (USDA, 2026). Milk is the primary source of animal protein for millions of lacto-vegetarians and plays a vital role in India’s food safety and nutritional strategy (Deshmukh and Paramasivam, 2016A).
       
Beyond its nutritional significance, the dairy sector is important contributor to the rural economy, supporting the livelihoods of more than 85 million people, mainly small and marginal farmers, with women accounting for nearly 62% of the active dairy workforce (GOI, 2026). Dairy farming contributes to poverty reduction, income diversification, livelihood resilience and women’s economic empowerment (Birthal and Negi, 2012; Tricarico et al., 2020; Lyngkhoi et al., 2022). In addition, dairy cooperatives have promoted inclusive growth and equitable resource allocation (Gaillard and Derville, 2022). Consequently, the dairy sector makes significant contributions toward achieving the United Nations Sustainable Development Goals (SDGs), particularly SDG 1 (No Poverty), SDG 2 (Zero Hunger) and SDG 5 (Gender Equality).
       
India’s dairy sector has expanded rapidly over recent decades due to population progress, urbanization and increasing disposable incomes, reflecting the country’s ongoing “Livestock Revolution” (Delgado et al., 1999; Birthal and Negi, 2012). The industry continues to exhibit sustained growth (GOI, 2026) and milk remains the only livestock commodity with a surplus in per capita availability (Srivastava et al., 2024). Moreover, the demand for milk and dairy products is highly responsive to income growth, exceeding that of other animal-source foods such as eggs, meat and fish (Vasant and Zhou, 2010). These trends highlight the importance of generating reliable production forecasts to support policy formulation, supply planning and investment decisions.
       
Forecasting milk production is challenging because it is influenced by numerous economic, technological, environmental and regional factors (Rao, 2017). In India, substantial interstate variation in socio-economic conditions, feeding practices, infrastructure and production systems contributes to significant differences in milk productivity and output (Rajeshwaran et al., 2014; Kale et al., 2019). At the same time, advances in statistics, econometrics and artificial intelligence have expanded the range of forecasting techniques available for modelling complex agricultural systems (Di Franco and Santurro, 2021). Among these, machine learning algorithms have demonstrated considerable potential in livestock applications, including disease detection (Bhardwaj et al., 2023), milk yield prediction (Kumar et al., 2019), livestock product forecasting (Kim et al., 2024) and methane emission estimation (Perumal et al., 2024).
       
Several studies have applied conventional forecasting approaches, including ARIMA, Holt’s exponential smoothing and related econometric models, to analyse milk production trends in India (Deshmukh and Paramasivam, 2016B; Jirli et al., 2025). Although machine learning techniques, particularly artificial neural networks, have increasingly been adopted in agricultural forecasting, comparative evidence evaluating their performance against traditional time-series models for long-term milk production forecasting remains limited. Therefore, the present study compares the forecasting performance of ARIMA, holt linear trend, exponential smoothing state space (ETS) and neural network autoregression (NNAR) models using annual milk production data for India to identify the most accurate and reliable model for long-term forecasting.
       
The contribution of this study lies in the systematic comparison of four fundamentally different forecasting approaches ARIMA, Holt’s linear trend, ETS and NNAR using India’s annual milk-production data for 1961-2024. The study aims to evaluate the temporal properties of milk production, compare the forecasting performance of the competing models, identify the most suitable approach and generate long-term projections for 2025-2030 with prediction intervals. The findings provide empirical evidence on the relative suitability of conventional and neural-network forecasting approaches and offer a basis for long-term milk-production planning and dairy-sector policy in India.

Data source and transformation
 
The study used annual milk-production data for India obtained from the FAOSTAT database. The FAOSTAT Item was “Milk, total” and the Element was “Production”, measured in tonnes, covering the period 1961-2024 (64 observations). The data were downloaded from FAOSTAT on 20 August 2025. For analysis, production quantities were expressed in million tonnes. The series was transformed using the natural logarithm to reduce scale effects and stabilize variance. No observations were removed from the original series and the transformed series was used for subsequent time-series analysis. All analyses were conducted in R statistical software using appropriate time-series and forecasting packages.
 
Stationarity analysis
 
Stationarity was first assessed to determine the order of integration of the annual milk production series. The augmented dickey-fuller (ADF) and kwiatkowski-phillips-schmidt-shin (KPSS) tests were employed to examine the stationarity of the data prior to model estimation (Pandey et al., 2025).  The ADF model is represented as:
 

 
Where,
Y_t= The logarithm of milk production at time t.
Δ= The first difference.
α= The intercept.
βt= The deterministic time trend.
γ= The coefficient associated with the lagged level of the series.
k= The number of lagged differences.
δ_i= Short-run autoregressive coefficients.
ε_t= The error term.
       
The null hypothesis of the ADF test is H0:γ=0, indicating a unit root, against H1:γ<0, indicating stationarity.

The KPSS test is calculated via:

 
Where,
S_t = The cumulative sum of the residuals up to time t.
n= The total number of observations.
The estimated long-run variance of the residuals.
       
The cumulative residual is defined as:

 
Where,
e_i= The residual obtained after removing the deterministic component from the series.
Where the hypotheses are formulated as:
H0 : Series is stationary.
H1 : Series is non-stationary.
       
The ADF test was estimated with an intercept and deterministic linear trend and the lag length was selected using the akaike information criterion (AIC) to account for serial correlation in the residuals. The KPSS test was specified with an intercept and deterministic trend and the newey-west bandwidth was selected using the automatic bandwidth selection procedure implemented in the R package. These specifications were maintained consistently across the level and differenced series.
 
Model identification and estimation
 
The autocorrelation function (ACF) and partial autocorrelation function (PACF) plots were examined to assess the serial dependence of the time series and assist in model identification. Significant autocorrelations in the ACF and PACF plots were used to determine the initial values of the non-seasonal autoregressive (p) and moving average (q) orders for the ARIMA model. The general ARIMA(p,d,q) model can be stated using the backshift operator notation:
 

ϕ (B) (1 - B)^d Y_t = θ (B) ε_t

 
Where,
B= Backshift operator.
d= Differencing order.
ϕ (B)= Autoregressive polynomial.
θ (B)= Moving average polynomial.
Y_t= The logarithm of milk production.
ε_t= The random error term.
       
The selected optimal framework layout was ARIMA (0,2,1), expressed as:
 

(1 - B)² Y_t = (1 - θ₁B)ε_t

 
Where,
θ₁= The MA(1) coefficient.
 
Training and validation framework
 
The 64 annual observations were divided into a training period from 1961-2011 (51 observations) and an out-of-sample validation period from 2012-2024 (13 observations). The competing forecasting models were estimated using the training period and forecasts for 2012-2024 were generated without using the actual validation observations during model estimation. Forecast accuracy was subsequently evaluated by comparing the predicted and observed values over the validation period.
 
Alternative forecasting models
 
The selected ARIMA model was compared with holt’s linear trend, exponential smoothing state space (ETS) and neural network autoregression (NNAR) models. The ETS model was specified as ETS(A,A,N), representing additive errors, additive trend and no seasonality. Holt’s linear trend model was used to capture level and trend components without a seasonal component because the data are annual.
 

The NNAR model was specified as: NNAR(p,k)

 
Where,
p= The number of lagged observations used as input variables.
k= The number of neurons in the hidden layer.
       
The neural network autoregression (NNAR) model was estimated using the nnetar() function in the R forecast package. The lag order (p) and hidden-layer size (k) were selected automatically by the function based on the training data.

Diagnostic checking
 
Model adequacy was assessed through diagnostic tests applied to the residuals. The ljung-box test determined whether residual autocorrelation remained:

 
Where,
n= The number of residual observations.
m= The number of lags tested.
r²_k= The sample autocorrelation of the residuals at lag k.
The null hypothesis is that the residuals are independently distributed up to lag m.
The jarque-bera test examined residual normality based on skewness (S) and kurtosis (K):

Where,
n= The number of observations.
S= The sample skewness.
K= The sample kurtosis.
The null hypothesis is that the residuals are normally distributed.
Conditional heteroscedasticity was monitored via the ARCH-LM test layout:

 
Where,
ε²_t = The squared residual at time t.
α₀= The intercept.
α_i= The ARCH coefficients.
q= The number of lags.
μ_t= The error term.
 
Forecast validation and alternatives
 
Forecast accuracy was assessed using root mean square error (RMSE), mean absolute error (MAE) and mean absolute percentage error (MAPE) (Singh et al., 2025).:

 
Where,
Y_t= The annual milk production in different years.
Ŷ_t= The forecasted annual milk production in the corresponding years.
n= The number of years used as forecasting period.
 
Diebold-mariano test
 
The diebold-mariano (DM) test was used to determine whether differences in forecast accuracy between competing models were statistically significant. The loss differential was defined as:
 

d_t = L (e_1t) - L (e_2t)

 
Where,
e_1t and e_2t= The forecast errors from the two competing models.
L(.)= The selected loss function.
The DM statistic is given by:

 
Where,
= The mean loss differential.
Its estimated variance.
       
The null hypothesis is H0:E(dt )=0, indicating equal predictive accuracy between the two competing models.
 
Volatility analysis
 
To investigate variability and persistence in milk production growth a generalized autoregressive conditional heteroskedasticity [GARCH (1,1)] model was estimated using log-derived growth rates. The conditional variance equation is:

 
Where,
h_t= Conditional variance,
ω= Constant variance,
α= ARCH parameter.
β= GARCH parameter.

Descriptive statistics
 
Descriptive statistics of India’s milk production (1961-2024) show an average annual production of 80,513.95 thousand tonnes, ranging from 19,247 to 247,843.90 thousand tonnes (Table 1). This substantial increase reflects the long-term expansion of the dairy sector, supported by initiatives such as the Operation Flood Programme, improved animal husbandry practices and the adoption of crossbred dairy animals (Dalal and Pathak, 2010).


       
The standard deviation (64,222.70×10³ tonnes) indicates considerable variation over the study period, largely driven by a sustained upward trend rather than cyclical fluctuations. A positive skewness of 1.110 suggests that higher production levels are concentrated in recent years, while the kurtosis value of 3.214 indicates a slightly leptokurtic distribution. The estimated compound annual growth rate (CAGR) of 4.05% demonstrates consistent long-term growth, exceeding historical benchmarks reported for other major milk-producing regions (Hyndman and Athanasopoulos, 2018).
 
Stationarity and model selection
 
The stationarity of the log-transformed buffalo milk production series was examined using the ADF and KPSS tests. The original series was non-stationary, while the first-differenced series showed mixed results (Table 2). After second differencing, both tests confirmed stationarity (ADF: p = 0.0100; KPSS: p = 0.1000). Therefore, the series was considered stationary at the second difference (d = 2) and used for ARIMA model estimation.


       
ARIMA (0,2,1) achieved the lowest AIC (-220.12) and BIC (-216.33) among the candidate ARIMA specifications (Table 3) and was therefore selected as the final ARIMA model based on in-sample information criteria. Although ARIMA (2,2,1) recorded a slightly lower in-sample RMSE (0.02380) than ARIMA (0,2,1) (0.02398), this difference did not alter the model selection based on AIC and BIC.


 
Parameter estimation and residual diagnostics
 
The MA (1) parameter was estimated at -0.7417 and was statistically significant at the 1% level (Table 4). The negative MA (1) coefficient indicates an inverse relationship between the current error and the preceding-period error, suggesting that short-term forecast errors tend to be partially offset in the subsequent period (Perez-Guerra et al., 2023). However, this parameter should be interpreted as a statistical error-correction property of the time-series model rather than as evidence of biological or economic self-correction in milk production.


       
Residual diagnostics confirmed model adequacy (Table 5). The ljung-box test confirmed the absence of significant serial correlation (P=0.8812). Residuals were approximately normally distributed supported by jarque-bera test (P=0.0877), Normal Q-Q Plot (Fig 1) and histogram of residuals (Fig 2). The ARCH LM test verified the absence of significant heteroscedasticity (P=0.2147), validating stable residual variance across time.






 
Comparative forecasting performance
 
Out-of-sample forecasting performance was evaluated over the validation period 2012-2024. ARIMA (0,2,1) achieved the lowest RMSE (12,693.94), MAE (10,407.38) and MAPE (4.97%) among the four competing models (Table 6 and Fig 3), indicating the highest predictive accuracy during the validation period. Simpler ARIMA specifications often match or outperform over-parameterized machine learning structures when capturing long-term agricultural production series (Balan et al., 2025; Pachava et al., 2025, Mishra et al., 2020).




       
Holt’s linear trend and ETS (A,A,N) produced identical accuracy measures (RMSE = 29,205.51; MAE = 23,207.55; MAPE = 10.86%), which were verified from the corresponding forecasts. The NNAR model recorded the highest forecast errors (MAPE = 19.29%). Thus, while AIC and BIC were used for in-sample ARIMA model selection, RMSE, MAE and MAPE were used to evaluate out-of-sample forecasting performance. Holt’s linear trend and ETS (A,A,N) are equivalent under the additive-error, additive-trend, no-seasonality specification and therefore produced identical forecasts and accuracy measures.
       
The Diebold–Mariano test proved that the predictive superiority of the ARIMA model was statistically significant (Table 7). Rejection of the null hypothesis occurred at the 1% level against both ETS (P=0.00013) and NNAR (P=0.00002). Given that holt’s linear trend and ETS (A,A,N) produced equivalent forecasts, their DM comparisons with ARIMA were also expected to be closely aligned.


       
Year-wise validation results (Table 8) show that ARIMA (0,2,1) maintained relatively good forecasting performance throughout the validation period. However, a tendency toward underestimation became more evident from 2016 onward, with the forecast error increasing from 2.68% in 2016 to a maximum of 9.09% in 2021 before declining to 6.58% in 2024. Despite this later-period deviation, ARIMA (0,2,1) achieved the lowest overall RMSE, MAE and MAPE among the competing models, indicating its superior overall out-of-sample forecasting performance. The observed underestimation suggests that the model may have been less responsive to the acceleration in milk-production growth during the later validation years and represents an important consideration when interpreting its long-term forecasts.

 
The GARCH (1,1) analysis was conducted as a supplementary assessment of conditional volatility in milk-production growth rates (Table 9). The diagnostic results indicated no significant residual serial correlation (Weighted ljung-box, p = 0.6730), no remaining ARCH effects (Weighted ARCH LM, p = 0.2378) and no evidence of sign bias (p = 0.4227). These findings provide no statistical evidence of persistent conditional heteroskedasticity in the growth-rate series. Therefore, GARCH modelling does not provide additional evidence for improving the primary milk-production forecasting model and is treated as supplementary to the main forecasting analysis.


       
Following the validation exercise, the selected ARIMA (0,2,1) model was re-estimated using the complete dataset and used to generate milk-production forecasts for 2025-2030 (Table 10). The model projects milk production to increase from 259,075.2 thousand tonnes in 2025 to 323,343.5 thousand tonnes by 2030. This trajectory supports long-term industry expansions if infrastructure changes match the increasing volume (Perez-Guerra et al., 2023; Balan et al., 2025). The corresponding 80% and 95% prediction intervals widen with the forecast horizon, reflecting increasing uncertainty in the projections. Therefore, these estimates should be interpreted as statistical baseline projections conditional on the historical data-generating process rather than deterministic predictions.

This study has several limitations. First, the analysis relies on annual milk-production data, resulting in a relatively limited number of observations for estimating complex nonlinear forecasting models. Second, the forecasting models are based primarily on the historical dynamics of milk production and do not explicitly incorporate potentially important explanatory factors such as animal population, productivity, feed availability, climatic conditions, technological change, input prices and policy interventions. Third, although ARIMA (0,2,1) provided the best overall validation performance, the year-wise results indicated some underestimation during the later validation period, highlighting the possibility of changing growth dynamics. Finally, uncertainty increases with the forecasting horizon, limiting the interpretation of long-term projections as precise predictions. Future research could therefore investigate hybrid and multivariate forecasting approaches incorporating climatic, technological, economic and livestock-sector variables, as well as higher-frequency data where available, to improve predictive performance and capture structural changes in milk production.

This study evaluated the forecasting performance of time-series econometric and machine learning models using annual milk production data for India covering the period 1961-2024. Among the competing models, ARIMA (0,2,1) demonstrated the highest forecasting accuracy and significantly outperformed ETS and NNAR models. The forecasts indicate that India’s milk production is projected to increase from 259.08 million tonnes in 2025 to 323.34 million tonnes by 2030. These estimates should be interpreted as statistical baseline projections rather than deterministic predictions, particularly because forecast uncertainty increases with the projection horizon. The findings highlight the need for continued investment in dairy infrastructure, animal health services, genetic improvement programs and feed and fodder development to sustain future growth. Strengthening dairy value chains, processing facilities and market access mechanisms will also be essential to accommodate rising production and enhance farmers’ incomes. Furthermore, incorporating forecasting-based planning into dairy sector policies can support timely decision-making and contribute to the long-term sustainability and competitiveness of India’s dairy industry.

The authors declare that there is no conflict of interest regarding the publication of this manuscript.