Objective: To determine the temporal patterns of cumulative incidence of brucellosis using autoregressive integrated moving average models. Methods: This cross-sectional study employed yearly and monthly data of 1 117 laboratory-confirmed human brucellosis cases from January 2013 to December 2018 using the Yazd brucellosis national surveillance system. The monthly incidences constructed a timeseries model. The trend of cumulative incidence was perceived by tracing a line plot, which displayed a seasonal trend with periodicity. Thus, the ARIMA models were selected. Thereafter, Akaike information criteria (AIC) and Bayesian information criterion (BIC) values among different models indicated a preferable model from models which were expanded by diverse lags [(3, 0, 3), (2, 0, 3), (3, 0, 2), (4, 0, 3) and (3, 0, 4)]. Then, the achieved ARIMA model was applied to the forecasting cumulative incidence of monthly brucellosis incidences. All analyses were performed using Stata, version 11.2. Results: For the ARIMA (3, 0, 4) model, MAPE value was 56.20% with standard error 0.009-0.016, and white noise diagnostic check (Q=19.79, P=0.975) for the residuals of the selected model showed that the data were completely modelled. The monthly incidences that were fitted by the ARIMA (3, 0, 4) model, with AIC (25.7) and BIC (43.35) with a similar pattern of actual cases from 2013 to 2018 and forecasting incidences from January 2019 to December 2019 were, respectively, 0.50, 0.44, 0.45, 0.49, 0.55, 0.58, 0.56, 0.51, 0.46, 0.44, 0.45 and 0.49 per 100 000 people. Conclusions: In summary, the study showed that the ARIMA (3, 0, 4) model can be applied to forecast human brucellosis patterns in Yazd province, supplementing present surveillance systems, and may be better for health policy-makers and planners.