Introduction

The cattle population has increased compared to the nineties, currently there are more than five million head, being the Creole type with the largest presence, followed by Brown Swiss, which represent 17.6 % of the total. This breed is bred mainly in the Andean zone for milk production¹. Dairy production is of great importance due to its weekly, biweekly and/or monthly payment, developing its activities mainly in areas of the Sierra, in these production systems the economic income reaches between 2000 to 15000 S per year, however, these amounts are below the Peruvian minimum wage².

Currently the development of dairy farming is in charge of farmers, societies, professionals and so far, it is still a challenge to increase milk production, studies of production systems in order to characterize and identify factors such as management, facilities, feeding and production records³. The analysis of records allows characterizing and analyzing milk production, being applied statistics an important tool for decision making. Currently, multivariate linear regression models are used to predict milk production, reporting little reliability; however, the application of time series models allows describing and predicting the behavior of a phenomenon that varies in time, and its dependence between successive observations⁴.

ARIMA (Autoregressive Integrated Moving Average) modeling, developed by Box-Jenkins, is a revolution in time series analysis^5,6. It is also a tool to recognize the main structural components in the temporal evolution of milk production and the lactation curve^7,8. This model becomes a simple and and flexible tool for predicting missing data to optimize the data to optimize the impact of genetic improvement programs on the genetic improvement programs in dairy cattle production. In this sense, the objective of the study was to apply the ARIMA model to forecast the milk production of Brown Swiss cows in the Peruvian Altiplano cows under Peruvian Altiplano conditions.

Materials and methods

Place of study and data systematization. The study was characterized as retrospective descriptive. Milk production records from the 2008 to 2016 campaigns of the Chuquibambilla Research and Production Center (CIPC) of the Faculty of Veterinary Medicine and Zootechnic of the National University of the Altiplano Puno (Peru) were used. The daily records of each individual were categorized by third of lactation. They were systematized in a spreadsheet (Microsoft Excel 365^®), the sum of the monthly productions of each year was calculated, including the year 2008 through 2016 (Table 1).

Application of ARIMA model. The model is characterized by three stages: identification, estimation and

diagnostic review. and diagnostic review:

Horizontal plotting of total milk production ("y" axis) and years ("x" axis) was used using the package or "library" of the RStudio program called ggplot and specifically with the command called autoplot (Figure 1).

A seasonal plot of milk production was applied for each production year with the difference of having on the y-axis the months of production using the ggseasonplot command (Figure 2).

The milk production forecasts for 2017 were made using the RStudio commands meanf (average), naive (naive), snaive (seasonal naive) and rfw (model with drift), with total production on the y-axis and production years on the x-axis (Figure 3).

Finally, the autoregressive model that combines the autoregressive (AR(p)) and moving average (MA(q)) processes was applied, making a graph of the forecasts (autoplot) with confidence intervals at 80 and 95 % (Figure 4). All the analyses were performed using the statistical program R v. 4.0.3 with its RStudio extension⁹.

Results

Figure 1 shows the behavior of milk production (L) of Brown Swiss cows from CIPC, Puno, during the years of the study.

Figure 1 Horizontal plot of total milk production in relation to the years 2008-2016 

Figure 2 Total seasonal milk production in relation to months (2008-2016) 

Figure 3 Milk production supplemented with simple forecasting methods 

Figure 4 ARIMA model forecasts and their prediction intervals 

Figure 2 shows milk production in relation to months for the nine years of the study.

The highest yields in three years coincide with the month of May while the lowest peaks are reported for the months of August and September.

Figure 3 expresses a horizontal plane of milk production similar to Figure 1 complemented with simple forecasting methods such as: naïve, seasonal naïve, drift model and the average model.

Figure 4 shows the forecast using the ARIMA models for the following 12 months (year 2017) observing a similar pattern with high peaks in months such as May and drops for September and October, respectively. Figure 4 shows milk production horizontally complementing the respective forecast (blue line). The 80 % confidence interval is shown as blue shading, while the 95 % confidence interval is shown above it as a light blue to gray shaded area.

Discussion

Milk production shows seasonality over time, without expressing any positive or negative trend with high and low peaks in milk production at different times of the year. However, an important growth is observed in the positive peak of production in 2015, which is due to the increase of cows in production. Likewise, significant drops in production are observed in mid-2008 and in the last months of 2013 and 2014. The condition hypothesis of being non-stationary was tested and accepted using the Dickey Fuller test (p=0.02811). Milk production and its behavior were similar to those reported by Sánchez et al.⁴, Sánchez López et al.¹⁰, however, Mishra et al.¹¹ indicate that milk production is stationary after the evaluation of the autocorrelation function (ACF) and partial autocorrelation function (PACF)¹², as shown in Figure 1. The seasonal (Figure 2) clearly shows dependence on climatic characteristics such as: rainy season (December - March), the highest milk production is observed, dry season (April - September), in most years there is a noticeable decrease in milk production. Studies in the Peruvian altiplano characterize the dry season as an extended period with strong winds and extreme temperatures that restrict the development of pastures and forage crops, while the dry season is a period with strong winds and extreme temperatures that restrict the development of pastures and forage crops, the transition period is a short period with relative rainfall that results in a slight improvement in climatic conditions, and finally the rainy season is characterized by higher precipitation and favorable environmental temperatures for the growth and development of pastures and forage crops¹³.

Simple forecasting methods are characterized by being surprisingly effective, such as the seasonal naïve method, which is very useful for seasonal data, where each forecast is set to be equal to the last observed value for the same season and month of the year. Figure 3 shows in light green the "seasonal naïve" forecast, very similar to productions of past years, which indicates that the known data was processed and used accurately, which generates an increase in the probability of obtaining an appropriate and efficient forecast¹⁴.

The global autocorrelation observed in Figure 4 between the residuals was verified by means of a correlogram analysis and the value of the "Q" statistic (Q=14.791 and p=0.6763) indicating that this phenomenon is similar to that reported by other authors⁵ when also analyzing lechera production characteristics. ARIMA models are characterized by short term predictions with the disadvantage of not being able to consider the long term.

However, this model is a good option to represent milk production and to establish a forecast¹². In conclusion, the ARIMA model proposed for the characterization of milk production was convenient, since it allowed forecasting the production for the year 2017, as long as the complete data (production per day, month and year) of milk production is available.