Research Article

Using of Autoregressive Integrated Moving Average (ARIMA) Model for Forecasting Milk Production of Dairy Cattle Farms in Dakahlia Governarate of Egypt

Asmaa A. Badr1*, Eman A. Abo Elfadl1, Mohammed M. Fouda1, Sayed M. Elsayed2

1Department of Husbandry and Development of Animal Wealth, Faculty of Veterinary Medicine, Mansoura University, Gomhoria St., Mansoura, P.O. box 35516, Egypt; 2Department of Applied statistics and Econometrics, Faculty of Graduate Studies for Statistical Research, Cairo University, Egypt.

Received | March 30, 2022; Accepted | May 11, 2022; Published | June 18, 2022

*Correspondence | Asmaa A. Badr, Department of Husbandry and Development of Animal Wealth, Faculty of Veterinary Medicine, Mansoura University, Gomhoria St., Mansoura, P.O. box 35516, Egypt; Email: asmaabadr@mans.edu.eg

Citation | Badr AA, Elfadl EAA, Fouda MM, Elsayed SM (2022). Using of autoregressive integrated moving average (ARIMA) model for forecasting milk production of dairy cattle farms in Dakahlia Governarate of Egypt. Adv. Anim. Vet. Sci. 10(7):1473-1480.

DOI | https://dx.doi.org/10.17582/journal.aavs/2022/10.7.1473.1480

ISSN (Online) | 2307-8316

Copyright: 2022 by the authors. Licensee ResearchersLinks Ltd, England, UK.

This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

INTRODUCTION

Milk production is the most attracting trait of dairy cows and continues to receive significant attention of researchers around the globe (Garamu, 2019). Milk yield is the most important determinant of profit for dairy cattle projects. Maximum milk yield from cows is considered the most important goal that could be achieved through 305 days lactation length and 60 days dry period (Snowder and Glimp, 1991). Milk production of dairy cows is not only affected by their genetic makeup but also is determined by many environmental factors such as breed, season of the year, lactation length, calving interval (CI), age of calving (AC), parity, stage of lactation, nutrition and days open (DO) (Susanto et al., 2019). Shortage of labor, emphasis on increasing farm efficiency and quality of life of the farmers were the driving factors for increasing milk production (Medeiros et al., 2022).

One of the most important benefits of milk is that it can tremendously improve the nutritional levels of children in Africa (Siddiky, 2015). Dairy products enterprises are increasing as the best trial to make profitable margins, so milk has grabbed the attention of governments trying to implement policies that could forecast its production and subproducts. Global milk production reached nearly 906 million tons in 2020, 2% increase from 2019 driven by output increases in all geographical regions except in Africa, where production remained stable. Milk volume increases were the highest in Asia followed by Europe while Africa and Oceania have the lowest milk production in the world (FAO, 2019). In Asia, milk output rose to 379 million tons in 2020, 2.6 % increase from 2019 this was principally driven by increases mainly in India, China, Pakistan and Turkey. As a commodity, India is the largest milk producing country, where milk output reached 195 million tons in 2020, 2 % increase from 2019 underpinned by the continued rise in dairy cattle numbers and improved feed and fodder availability on favourable monsoon rains (June to September). Egypt produced about 6.6 tons in year of 2020 (FAO, 2020). Holstein cows are the most producing breed of milk when compared to other breeds, as it can withstand well under adverse conditions and maintain high production levels (Nawaz et al., 2013).

Forecasting of milk production is required so that necessary policy formations can be done and strategic decision can be taken to enhance dairy development (Mishra et al., 2020a). There are several forecasting methods which might vary from sector to other sector and local need as Simple Average Growth Rate (SAGR), Compound Average Growth Rate (CAGR), Exponential Growth Rate (EGR), Autoregressive Integrated Moving Average (ARIMA) and Holt’s Linear Models (Gooijer and Hyndam, 2006). Predicting milk production is the best tool to adjust its supply due to the importance of milk as a dairy product. Since South-Asian countries are the leading countries in milk production, they try to forecast milk production using ARIMA/GARCH models and Holt’s Linear Model (Oliveros, 2019). In a study conducted by Akhter and Rahman (2010), they forecasted milk supply up to 3 years for a dairy cooperative in the United Kingdom. While Murphy et al. (2014) and Zhang et al. (2020) conducted a study to identify the different modeling techniques for the prediction of total daily herd milk yield and non-linear model was used especially for short-term milk yield predictions. Mishra et al. (2020b) used time series models as ARIMA and VAR methods in milk production and forecasted milk production in India for year 2024-2025. Moreover, ARIMA approach indicates that India would be the leading country in milk production with 91 million tons in the year 2024–2025 among South Asian countries (Pal et al., 2007). The second ranked country is Pakistan, whose milk production would reach 26 million tons in 2024–2025. China is the third country with 3 million tons, while Bangladesh and Sri Lanka seem to be the countries with the lowest milk production (Deshmukh and Paramasivam, 2016). The increase in milk production of dairy cows is going to be low even though the government policies due to a number of reasons such as the low genetic capacity of the indigenous cattle for milk production, low adaptation ability of exotic and hybrid dairy cows, substandard feeding, poor health care and high cost of inputs (Abunna et al., 2018). Therefore, forecasting of milk production is an important strategy to decrease inputs and increase farmers income, which constitute an important portion of dairy industry.

MATERIALS AND METHODS

Ethical statement

The study protocol was approved by the Research Ethics Committee of the Faculty of Veterinary Medicine, Mansoura University, Egypt. Data were collected from farms during the presence of farms owners and farm administration.

Data collection

Data were collected from accurate records in the farms or by research questionnaire methods that were conducted when there were no records in the farms (Atallah, 1997). A total of 1534 lactation records were collected from accurate records in the farms and 266 lactation records were collected by research questionnaire. Incomplete records or pedigree files with unclear information were excluded from the data sets. The original set of data consists of 1800 lactation records of prevalent Holstein Friesian and Holstein German cows from 2013-2021 during COVID-19 occurrence. A total of 900 lactation records were belonged to Holstein Friesian and the other 900 lactation records were belonged to Holstein German. Data were collected from Albayoumi farms in Dakahlia governorate, which is located in Egypt (N 29° and E 25.48°) according to GPS reading. Dakahlia governorate is present in the east of the Delta of the Nile and covers about 3.459 km2. It locates in a very strategic location overlooking Damietta branch of the River Nile and the Mediterranean Sea coast and boarded with El-Sharkia governorate from the east, El-Gharbia Governorate from the west and Damietta governorate to the northwest.

Studied variables and statistical analysis

Total milk production of the farm was calculated from total daily milk production through total 305 days lactation curve during period of 2013–2021. The statistical analysis by ARIMA was performed using statistical software (STAT GRAPHICS centurion, version 17).

ARIMA model

ARIMA model is known as Box-Jenkins method, who developed a coherent versatile three-stage iterative cycle for time series identification, estimation and verification (Box and Jenkins, 1976). ARIMA method explained the movement of a variable by its past or lagged values. It produces predictions based on the synthesis of time series data. It helps to analyze both probabilistic and stochastic properties of time series data. A time series containing records of a single variable is termed as univariate. But if records of more than one variable are considered, it is termed as multivariate. A time series can be continuous or discrete. In a continuous time, series, observations are measured at every instance of time as temperature reading, flow of a river and concentration of a chemical process. Meanwhile, discrete time series contains observations measured at discrete points of time as production of farms and exchange rates between two different currencies. In discrete time series the consecutive observations are recorded at equally spaced time intervals such as hourly, daily, weekly, monthly or yearly time separations as mentioned by Kantz and Schreiber (2004). ARIMA can be done on single and multiple variables (Beck and Katz, 2011). ARIMA model is most widely used for forecasting milk production elsewhere. ARIMA (p, q, d) model where ‘p’ is the order of the autoregressive part (AR), ‘d’ donates the degree of differencing involved and ‘q’ is the order of the moving average part (MA).

Autoregressive model (AR)

AR is a linear regression model that uses its own lags as predictors. AR equation is:

Yt = 𝜇 + Ø1Yt-1 + Ø2Yt-2 + …… + ØpYt-p + 𝜀t

Where; Yt is milk production, 𝜇 is constant, Ø1, Ø2.... Øp are the parameters of the model and 𝜀t is independently and normally distributed with zero mean and constant variance for t = 1.

Moving average model (MA)

The notation MA (q) refers to the number of lagged forecast errors that should go into the ARIMA Model. MA equation is:

Yt = 𝜇 + θ1 𝜀t-1 + θ2 𝜀t-2+………...+ θp 𝜀t-p + 𝜀t (Fan and Yao, 2008)

Where;­ Yt is milk production, θ1, θ2, ..., θq are the parameters of the model, μ is the expectation of Yt (often assumed to equal 0).

The general form of ARIMA model of order (p, q, d)

Yt = 𝜇 + Ø1 Yt-1 + Ø2 Yt-2 + ……+ Øp Yt-p + 𝜀t + 𝜇 + θ1 𝜀t-1 + θ2 𝜀t-2+………...+ θp 𝜀t-p + 𝜀t (Tsay and Tiao, 1984)

RESULTS AND DISCUSSION

ARIMA model includes the following steps for Holstein Friesian

Model identification

At first, the data is checked for stationarity with the help of the autocorrelation (ACF) and partial autocorrelation function (PACF). As shown in Figures 1 and 2, the data is located between 0.5 and – 0.5, indicating the existence of stationarity assumption. A stationarity process can be defined in precise mathematical process is a flat looking series without trend, constant variance over time, constant autocorrelation over time and no periodic fluctuations. Consequently, parameters such as mean and variance also do not change over time.

 

 

Model estimation

It means estimation of model fitness through estimation of best coefficients as ACI, MSE, MAPE and SBIC. The results presented in Tables 1 and 2 indicated that ARIMA

 

Table 1: Model fit measures.

SBIC	AIC	MPE	ME	MAPE	MAE	RMSE	Model
21.24	21.43	0.2427	1579.74	2.66175	14986.6	28935.2	ARIMA (2, 1, 2)

 

Table 2: Diagnostic checking.

SBIC	AIC	MPE	ME	MAPE	MAE	RMSE	Model
21.52	21.43	0.2427	1579.74	2.66175	14986.6	28935.2	ARIMA (2, 1, 2)
21.79	21.72	1.03746	5779.02	3.71488	20844.4	37468.9	ARIMA (1, 1, 2)
22.16	22.07	0.07011	1941.76	3.70796	20690	39813.3	ARIMA (2, 1, 2)
22.38	22.36	2.20922	15591.6	6.98537	43238	64267.7	ARIMA (2, 1, 2)

 

(2, 1, 2) model is the best fit model because it is the model with lowest values of fit measures as RMSE, MAE, MAPE, ME, MPE, AIC and SBIC.

AIC = -2 log L + 2m

Where, L is the likehood function, m = p + q

A previous study conducted by Taye et al. (2020) differs from our study as the most suited model suggested by their study was ARIMA (1, 2, 1). On the other hand, the results suggested by Sankar and Prabakaran (2012) and Chaudhari and Tingre (2013) showed that ARIMA (1, 1, 0) is the most fit model. ARIMA (1, 2, 1) for the series of Culture purebred milk production, was the best fit model and the best fit model for the series of indigenous milk production was ARIMA as results suggested by (Yonar et al., 2022).

 

Diagnostic checking

For adequacy of the model, the residuals are examined from the fitted model as indicated in Figures 3 and 4. Different ARIMA models are considered, if necessary. If the first models are tried until a satisfactory model fits to the data.

Table 3 indicated ARIMA model contents of AR and MA and also indicated AR (1) (- 0.365357) with lower P value (0.000562) than AR (2) and MA (1) (0.91723) also has the lower P value (0.002591) than MA (2). So, ARIMA model (2, 1, 2) equation is described as follows:

Yt = - 0.365357 Yt-1 +0.91723 𝜀t-1+ 𝜀t

 

Table 3: ARIMA model contents

Parameter	Estimate	Standard error	P-value
AR (1)	-0.365357	0.0365291	0.000562*
AR (2)	0.328785	0.0939647	0.024918*
MA (1)	0.91723	0.137025	0.002591*
MA (2)	-2.65767	0.442846	0.003879*

 

P-value* is significant at 0.05 (p ≤ 0.05).

 

Model forecasting

Seven year forecast from 2022 to 2033 was done. Forecasting milk production for Holstein Friesian indicates that milk production will be increased at year 2022 and 2033 as indicated in Table 4. Milk production will be increased till reach 574313 kg in 2025 and 574183 kg in 2033. This indicates presence of variations in milk production forecasting between increase and decrease in future years.

 

Table 4: Milk production forecasting table for Holstein Friesian.

Forecast	Period
570192	2022
573511	2023
573511	2024
574313	2025
573703	2026
574339	2027
573906	2028
574274	2029
573997	2030
574219	2031
574047	2032
574183	2033

 

Table 5 indicates the milk production of Holstein Friesian from 2013-2021 and forecasted milk which calculated depended on its lagged value, also indicates the residual milk production which calculated from forecasted and collected data of milk.

 

Table 5: Holstein Friesian milk production (Kg).

Period	Data	Forecast	Residual
2013	500000
2014	499321	493666	5654.8
2015	587934	542917	45017.1
2016	499123	529073	-29950.0
2017	698450	707817	-9366.91
2018	532456	525419	7036.74
2019	612980	627290	-14310.4
2020	567812	560811	7001.14
2021	567891	566336	1555.41

 

Figures 5 and 6 showed time sequence plot of milk production combining between actual and forecasted milk.

 

 

ARIMA model for Holstein German in the second farm

Model identification

At first, the data were checked for stationarity with the help of the autocorrelation function (ACF) and partial autocorrelation function (PACF). Looking out to Figures 7 and 8, we found that the data located between 0.5 and – 0.5, which indicates that data are stationarity.

Diagnostic checking

For adequacy of the model, the residuals are examined from the fitted model as shown in Figures 9 and 10.

 

Table 8 indicated ARIMA model contents of AR and MA and also indicated AR (1) ( -0.0331079) and MA (1) with a coefficient -0.997. So, ARIMA model (0, 1, 2) equation is described as follows:

Yt = - 0.0331079 Yt-1 – 0.997 𝜀t-1+ 𝜀t

 

Table 6: Model fit measures.

SBIC	AIC	MPE	ME	MAPE	MAE	RMSE	Model
23.02	22.98	0.029	-493.002	7.05789	50031.1	7845.3	ARIMA (0, 1, 2)

Table 7: Diagnostic checking.

SBIC	AIC	MPE	ME	MAPE	MAE	RMSE	Model
23.02	22.98	0.029	-493.002	7.05789	50031.1	78458.3	ARIMA (0, 1, 2)
23	23	-2.667	-12634.3	13.02	90519.8	98934.4	ARIMA (0, 1, 0)
23.22	23.20	-0.925	-988.26	11.42	79997.1	97742.5	ARIMA (1, 0, 0)
23.32	23.30	-2.752	-12152.4	12.08	83237.0	102908	ARIMA (1, 1, 0)

 

Table 8: ARIMA model contents for Holstein German.

Parameter	Estimate	Standard Error	P-value
AR (1)	-0.0331079	0.105807	0.00940*
MA (1)	-0.997	0.267407	0.004701*

 

P-value* is significant at 0.05 (p ≤ 0.05).

 

 

Model estimation

It means estimation of model fitness through estimation of best coefficients as ACI, MSE, MAPE and SBIC. The results showed in Tables 6 and 7 indicated that ARIMA (0, 1, 2) model is the best fit model because it is the model with lowest values of fit measures as RSMSE, MAE, MAPE, ME, MPE, AIC and SBIC. Uddin et al. (2020) have forecasted the volume of milk in Andassa dairy farm in Ethiopia using ARIMA (1, 1, 1). Moreover, a study mentioned by (Yonar et al., 2022) suggested that ARIMA (1, 2, 1) is the best fit model for cross breed milk production.

Model forecasting

Forecasting milk production for Holstein German indicated that milk production will be increased at year 2022 and 2023 and then steadily increase at the following years as shown in Table 9. Table 10 indicates the milk production of Holstein German from 2013-2021 and forecasted milk which calculated depended on its lagged value, also indicates the residual milk production which calculated from forecasted and collected data of milk.

 

 

Forecasting plot in comparison with actual milk production was showed in Figures 11 and 12.

 

Table 9: Milk production forecasting table for Holstein German.

Forecast	Period
688840	2022
715748	2023
715748	2024
715748	2025
715748	2026
715748	2027
715748	2028
715748	2029
715748	2030
715748	2031
715748	2032
715748	2033

 

Table 10: Milk production of Holstein German(kg).

Period	Data	Forecast	Residual
2013	805678
2014	768450	844969	-76519.0
2015	645123	645382	-259.477
2016	567324	555626	11698.2
2017	601732	567408	34324.2
2018	745671	616549	129122
2019	678901	790088	-111187
2020	812096	826227	-14130.9
2021	704604	681596	23008.3

 

 

 

Conclusion and Recommendations

This study applied ARIMA model by using statistical graphics program to forecast milk production of Holstein Frisian and Holstein German. ARIMA forecasting results indicated that ARIMA (2, 1, 2) is the best fit model for Holstein Friesian in the first farm. Meanwhile, the ARIMA (0, 1, 2) is the best model for Holstein German in the second farm. The results showed that milk production in 2024 will be decreased in two farms in Egypt so that farms should take attention to increase milk production.

Acknowledgements

We would thank farms owner and farm administration for their help in our study. We also thank the staff at the Department of Husbandry and Development of Animal Wealth, Faculty of Veterinary Medicine, Mansoura University, Egypt for providing materials used in this study.

Novelty Statement

Our study is the first study to emphasizes the use of Autoregressive Integrated Moving Average (ARIMA) Model for Forecasting Milk Production of Dairy Cattle Farms in Dakahlia Governarate of Egypt.

Author’s Contribution

Asmaa A. Badr designed the study protocol, supervised data collection and analysis of data. Eman A. Abo El-fadl and Sayed M. Elsayed analyzed the data and shared in study protocol. Mohammed M. Fouda, Eman A. Abo El-fadl and Sayed M. Elsayed shared in writing the manuscript. All authors have finalized the experimental design and revised the manuscript and then contributed to, edited, and approved the final manuscript as submitted.

Conflict of interests

The authors have declared no conflict of interest.

References

Abunna F, Merid B, Goshu G, Waktole H, Mammo G (2018). Assessment of major reproductive health problems, their effect on reproductive performances and association with brucellosis in dairy cows in Bishoftu town, Ethiopia. J. Dairy Vet. Anim. Res., 7: 14-20. https://doi.org/10.15406/jdvar.2018.07.00183

Akhter SH, Rahman S (2010). Milk supply forecasting for a dairy cooperative in the UK using univariate time series modelling techniques, in book. Agric. Econ. New Res., pp. 227-242.

Atallah ST (1997). Economic and productive efficiency of veterinary management in dairy farms. Ph.D. thesis, Fac. Vet. Med., Alex. Univ. Egypt. Poult. Sci. J., 1: 435-447.

Beck N, Katz JN (2011). Modelling dynamics in time series cross section political economy data. Ann. Rev. Polit. Sci., 14(1): 331-352. https://doi.org/10.1146/annurev-polisci-071510-103222

Box GE, Jenkins GM (1976). Time series analysis: Forecasting and control San Francisco. Calif: Holden-Day.

Chaudhari DJ, Tinger AS (2013). Forecasting of milk production in India. An application of ARIMA model. Indian J. Dairy Sci., 66: 72-78.

Deshmukh SS, Paramasivam R (2016). Forecasting of milk production in India with ARIMA and VAR time series models. Asian J. Dairy Food Res., 35(1): 17-22. https://doi.org/10.18805/ajdfr.v35i1.9246

Fan J, Yao Q (2008). Nonlinear time series: Non parametric and parametric methods., pp. 89-94.

FAO Food outlook (2019). Dairy market review report. Rome, Italy.

FAO Food outlook (2020). Biannual report on global food markets. Italy, Rome.

Garamu K (2019). Significance of feed supplementation on milk yield and milk composition of dairy cow. Vet. Sci. J., 13(2): 555–860.

Gooijer JG, Hyndam RJ (2006). 25 years of time series of forecasting. Int. J. Forecast., 22: 443-473. https://doi.org/10.1016/j.ijforecast.2006.01.001

Kantz H, Schreiber T (2004). Nonlinear time series analysis (Vol. 7). Cambridge university press. https://doi.org/10.1017/CBO9780511755798

Medeiros I, Fernandez-Novo A, Astiz S, Simoes, J (2022). Historical evolution of cattle management and herd health of dairy farms in OECD countries. Vet. Sci., 9(3): 125. https://doi.org/10.3390/vetsci9030125

Mishra P, Faith C, Niranjan HK, Tiwari S, Dubey A (2020a). Modelling and forecasting of milk production in Chhattisgarh and India. Indian J. Anim. Res., 54(7): 912-917. https://doi.org/10.18805/ijar.B-3918

Mishra P, Faith C, Vani GK, Tiwari S, Ramesh D, Dubey A (2020b). Time series investigation of milk production in major states of India using ARIMA modeling. J. Anim. Res., 10(1): 77-84. https://doi.org/10.30954/2277-940X.01.2020.10

Murphy DM, Zhang F, Upton J, Shine PA (2014). Review of milk production forecasting models, in book: Dairy farming operations management, animal welfare and milk production, 1st edition, Nova Science Publishers, pp. 14–61.

Nawaz A, Nizamani AH, Marghazani IB, Fatih A (2013). Influence of genetic and environmental factors on lactation performance of Holstein-Friesian cattle in Baluchistan. J. Anim. Plant Sci., 23(1): 17–19.

Oliveros RM (2019). The dairy industry in South east Asia. Persp. Chall. Opport. Earth Environ. Sci., 372: 1-9. https://doi.org/10.1088/1755-1315/372/1/012068

Pal S, Ramasubramannian V, Mehata SC (2017). Statistical Models for forecasting milk production in India. J. Ind. Soc. Agric. Statist., 61: 80-83.

Sankar TJ, Prabakaran R (2012). Forecasting milk production in Tamil Nadu. Int. Multidisciplin. Res. J., 2: 10-15.

Siddiky MNA (2015). SAARC dairy outlook. SAARC Agriculture Centre, Dhaka, Bangladesh, pp. 1-160.

Snowder G, Glimp H (1991). Influence of breed, number of suckling lambs, and stage of lactation on ewe milk production and lamb growth under range conditions. J. Amin. Sci., 69(3): 923-930. https://doi.org/10.2527/1991.693923x

Stat graphics centurion, version 17.1.06 (32 bit). Edition: multilingual.

Susanto A, Hakim L, Suyadi S, Nurgiartiningsih VMA (2019). Environment (year and season of birth) effects on first-lactation milk yield of dairy cows. November 2019. IOP Conf. Ser. Earth Environ. Sci., 372(1). https://doi.org/10.1088/1755-1315/372/1/012010

Taye BA, Alene AA, Nega KA, Yirsaw GB (2020). Time series analysis of cow milk production at Andassa dairy farm West Gojam Zone Amhara Region. Ethiopia, Modeling Earth System and Environment, pp. 364. https://doi.org/10.1007/s40808-020-00946-z

Tsay RS, Tiao GC (1984). Consistent estimates of autoregressive parameters and extended sample autocorrelation function for stationary and non-stationary ARMA models. J. Am. Stat. Assoc., 79(385): 84-96. https://doi.org/10.1080/01621459.1984.10477068

Uddin MM, Akter A, Khaleduzzaman ABM, Sultana MN (2020). Forecasting milk production in Bangladesh toward achieving self-sufficiency. Livest. Res. Rural Dev., 32(5).

Yonar H, Yonar A, Mishra P, Abotaleb M, Alkhatib AMG, Makarovskikh T, Cam M (2022). Modeling and forecasting of milk production in different breeds in Turkey. Indian J. Anim. Sci., 92(1): 105-111.

Zhang F, Shine P, Upton J, Shaloo L, Murphy MD (2020). A review of milk production forecasting models: Past and future methods.