Unfolding the potential of the ARIMA model in forecasting maize production in Tanzania

Abstract

Purpose

This paper was set to develop a model for forecasting maize production in Tanzania using the autoregressive integrated moving average (ARIMA) approach. The aim is to forecast future production of maize for the next 10 years to help identify the population at risk of food insecurity and quantify the anticipated maize shortage.

Design/methodology/approach

Annual historical data on maize production (hg/ha) from 1961 to 2021 obtained from the FAOSTAT database were used. The ARIMA method is a robust framework for forecasting time-series data with non-seasonal components. The model was selected based on the Akaike Information Criteria corrected (AICc) minimum values and maximum log-likelihood. Model adequacy was checked using plots of residuals and the Ljung-Box test.

Findings

The results suggest that ARIMA (1,1,1) is the most suitable model to forecast maize production in Tanzania. The selected model proved efficient in forecasting maize production in the coming years and is recommended for application.

Originality/value

The study used partially processed secondary data to fit for Time series analysis using ARIMA (1,1,1) and hence reliable and conclusive results.

Citation

Lwaho, J. and Ilembo, B. (2023), "Unfolding the potential of the ARIMA model in forecasting maize production in Tanzania", Business Analyst Journal, Vol. 44 No. 2, pp. 128-139. https://doi.org/10.1108/BAJ-07-2023-0055

1. Introduction

Zea mays is regarded as the major food crop produced and consumed worldwide (Nyaligwa, Hussein, Laing, Ghebrehiwot, & Amelework, 2017). It can be consumed by humans and livestock and used as raw materials for biofuel production. In Sub-Saharan Africa (SSA), maize is the most crucial primary cereal crop grown in over half of the countries and one of the top two cereals in over three-quarters of these countries (Suleiman & Kurt, 2015; Faostat, 2021). The crop flourishes on soils with pH between 5.0 and 7.0; nevertheless, a moderately acidic environment of pH 6.0 and 7.0 is mostly favorable (Baijukya et al., 2020).

Tanzania has been among the world's top 25 maize-producing countries (Twilumba, Ahmad, & Shausi, 2020). According to Suleiman and Kurt (2015), Tanzania is among the significant maize producers in SSA. It is a famous and prominent staple food (Laudien, Schauberger, Makowski, & Gornott, 2020), and both rural and urban population consumes it (Baijukya et al., 2020).

Food crop production, particularly maize, has a direct effect on the status of food security (Mkonda & He, 2017), which is one of the major areas of concern all over the world due to its contribution to all forms of human survival (Quaye, Yawson, Ayeh, & Yawson, 2012; Ngongi & Urassa, 2014). The significance of food security cannot be overemphasized, and in recent decades, the developing world has experienced food shortage, which led to food insecurity (Rwanyiziri et al., 2019). It is estimated that about 355 million people in SSA will be under food shortage by 2050 (Rwanyiziri et al., 2019).

In recent years, forecasting food crop production has become more challenging (Liu & Basso, 2020) due to several major drivers, particularly climate extremes like heavy rains, storms and floods. It is clear that major players' efforts to ensure the world has enough food are shattered by several major drivers, particularly climate change (WHO, 2021). Despite the efforts from stakeholders and the government, the country still has not done well regarding crop yield sustainability and food security (Mkonda & He, 2017). Although modelling cereals production has attracted intensive research due to the vitality of major food crops, little is known, especially in maize yield forecasting in response to food insecurity.

Time series methods have been widely used in forecasting future values based on past observations (Enders, 2015). Time series analysis involves studying the variables on which observations are arranged sequentially over time. In most cases, forecasting concentrates on univariate time series models pioneered by Box and Jenkins, including autoregressive (AR), moving average (MA) and autoregressive integrated moving average (ARIMA) models (Box & Jenkins, 1970). This method has been successful in many applications, including economics (Petrevska, 2017; Yildiran & Fettahoğlu, 2017), agriculture (Uwamariya & Ndanguza, 2018; Mgaya, 2019; Bezabih, Wale, Satheesh, Fanta, & Atlabachew, 2023), climate and environment (Fayaz, Meraj, Khader, & Farooq, 2022).

Few studies exist in Tanzania that model and forecast maize production using different methods. For example, Ogutu, Franssen, Supit, Omondi, and Hutjes (2018) forecasted maize production in Tanzania using dynamic ensemble seasonal climate forecasts. Laudien et al. (2020) used LASSO regression to forecast maize yields before harvest. Liu and Basso (2020) forecasted maize yields for smallholder farmers in three selected regions in Tanzania by integrating a field-based survey with a process-based mechanistic crop simulation model. In Mkonda and He (2018), the analysis focused on the trend of yields by plotting graphically, and the results were unreliable and inconclusive. This paper attempts to model annual maize production and forecast future production using the ARIMA model for the next ten years. The ARIMA method is chosen due to its proven forecasting reliability and ability to predict sequential series accurately.

2. Materials and methods

2.1 Data and data source

The paper used time series data of maize yield (hg/ha) from 1961 – 2021 obtained from FAO STAT calculated annually. Maize was selected due to its potential as the primary food crop for most households in Tanzania (Laudien et al., 2020). The dataset contained 60 observations, meeting the requirement of the generalization ability of time series analysis. Hyndman and Athanasopoulos (2018) argued that there is no justification for the minimum number of observations for ARIMA modelling, the only theoretical limit being that there should be more observations than the number of parameters in the forecasting model. Time series analysis was done using forecast, ggplot2 and ur packages ensembled in R software. The function auto. arima ( ) was used to obtain the appropriate ARIMA model with estimated parameters automatically instead of specifying using autocorrelation function (ACF) and partial autocorrelation function (PACF), as suggested by Hyndman and Athanasopoulos (2018).

2.2 The model formulation

The ARIMA method was considered appropriate to achieve the maize production forecasting objective. ARIMA is a family of time series models introduced by Box and Jenkins (1970) and is a widely used method in modelling and forecasting stationary and no-stationary time series with non-seasonal components. The reason for choosing the ARIMA technique is its ability to generate high predictive accuracy compared to the AR and MA as standalone models for short-run forecasting (Box & Jenkins, 1970). This methodology involves three successive phases: identification, which determines the order of the model required (p, d and q) to capture the data's salient dynamic features. This mainly leads to the use of graphical procedures (plotting the series, the ACF and PACF, etc.), estimation involves estimating the parameters of the different models. It proceeds to a first selection of models (using information criteria), and the diagnostic checking involves determining whether the model(s) specified and estimated is adequate. Notably, one uses residual diagnostics.

2.2.1 Autoregressive (AR)

An AR model is the one in which Yt depends only on its own past values Yt−1,Yt−2,Yt−3 e.t.c

A typical representation of an AR model where it depends on p of its past values called AR (p) is as shown below:

(1)Yt=α+ψ1Yt−1+ψ2Yt−2+ψ3Yt−3+....+ψpYt−p+εt

(2)Yt=ψ1Yt−1+ψ2Yt−2+ψ3Yt−3+….+ψpYt−p+εt

the Yt−1,Yt−2,…,Yt−p is predictor variable at lags t−1,t−2,…,t−p, α=μ(1−ψ1−ψ2−ψ3−…−ψp), ψ1,ψ2,ψ3,…,ψp are the parameters to be measured, and εt is the error term which follow white noise process with zero mean and a constant variance, i.e. εt ∼N(0,σ2). This AR (p) can also be presented as follows using backward shift operator;

(3)(1−ψ1B−ψ2B2−ψ3B3−ψpBp)Yt=εt

where B is the backshift operator defined by

BmYt=Yt – m (m =0,1,2 . . . p)

2.2.2 Moving average (MA)

An MA model is the one when Yt depends only on the error terms which follow a white noise process. The general form is given by:

(4)Yt=δ+εt−θ1εt−1θ2εt−2−θ3εt−3−...−θpεt−p

(5)Yt=εt−θ1εt−1θ2εt−2−θ3εt−3−...−θpεt−p

For δ=0

εt ∼WN(0,σ2)

Using the backshift operator, the above equation (5) can be rewritten as

(6)(1−θ1B−θ2B2−θ3B3−θpBp)εt =Yt

2.2.3 ARMA model

ARMA (p, q) is defined as the combination of AR(p) and MA(q) for a stationary time series and is given by the following equation:

(7)Yt =α+ψ1Yt−1+ψ2Yt−2+....+ψpYt−p+εt−θ1εt−1−θ2εt−2−θ3εt−3−...−θpεt−p

For α=0

(8)Yt =ψ1Yt−1+ψ2Yt−2+....+ψpYt−p+εt−θ1εt−1−θ2εt−2−...−θpεt−p

Equation (8) can be defined using backshift operator as follows,

(9)ψ(B)Yt=θ(B)εt

Although ARMA (p, q) is useful in modelling the time series process, it is special only if the series is stationary. Most time series data are not stationary, making ARMA unsuitable for those circumstances. Box-Jenkins methodology outlined a solution for the case of non-stationary series where data transformation (differencing) is needed to attain normality. Under this circumstance, the ARIMA (p, d, q) model is appropriate.

2.2.4 ARIMA model

ARIMA model is a combination of AR, i.e. Autoregressive (lagged observations as inputs), whereby I stand for Integrated (differencing to make series stationary) and MA, i.e. moving average (lagged errors as inputs).

According to Box and Jenkins (1970), the ARIMA model is denoted by ARIMA (p, d, q) where p is the order of autoregressive process, d is the order of integration, i.e. the number of differences to make the series stationary and q is the order of MA process. The general form of the ARIMA (p, d, q) is;

(10)Yt′ =α+ψ1Yt−1′+ψ2Yt−2′+....+ψpYt−p′+εt−θ1εt−1−θ2εt−2−θ3εt−3−...−θpεt−p

For α=0

(11)Yt′ =ψ1Yt−1′+ψ2Yt−2′+....+ψpYt−p′+εt−θ1εt−1−θ2εt−2−...−θpεt−p

Equation (8) can be defined using backshift operator:

(12)ψ(B)(1−B)dYt=θ(B)εt

where Yt′ is differenced response variable, d is the number of difference and εt ∼WN(0,σ2).

The differencing of the response variable can be calculated by using the following relation for non-stationary data:

Yt′ =∇dYt=(1−B)dYt=∑j=0d(dj)(−1)jYt−j

Box and Jenkins methodology (step-by-step procedures) were used to achieve the study's objective, which involves model identification, estimation, diagnostic checking and forecasting.

2.3 Model identification

Before ARIMA modelling, the time series data structure was checked for stationarity. Here, the focus was to observe the behavior of the mean and variance of a stochastic process to identify the existence of trends or seasonal patterns. As Montgomery, Jennings, and Kulahci (2015) suggested, checking for stationarity is essential because it brings equilibrium and stability to data. A series is considered stationary if its mean and variance are constant over time, and covariance depends only on lags (Enders, 2015).

2.3.1 Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test

Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test was used to determine if the series is stationary. The null hypothesis was that the times series level is stationary (Kwiatkowski, Phillips, Schmidt, & Shin, 1992). Reject null hypothesis for a large value of KPSS test statistic compared to critical values (Müller, 2005). This suggests that data transformation is required, i.e. log-transformation or differencing. The KPSS test statistic is given by

KPSS=1T2σ^T2∑t=1TSt2

Where σ^T2 is a consistent estimator of the long-run variance of the residuals, St=∑s=1te^s is the sum of the residuals, T is the number of observations (Kwiatkowski et al., 1992).

2.3.2 Correlograms

The autocorrelation function (ACF) and partial autocorrelation function (PACF) provided essential information to identify if a time series is stationary. The ACF defines the order of the AR process, while the PACF defines the order of the MA process. These two correlograms provided useful information concerning the stationarity of the time series process. For the time series process to be stationary, the ACFs for the AR process should be characterized by decaying exponentially tails off towards zero while the MA process cuts off to zero after lag q. In terms of PACFs, the AR process cuts to zero at lag p while the MA process decays exponentially tails off towards zero to show that the process is stationary (Wei, 2006).

2.3.3 Autocorrelation function (ACF)

The ACF measures the correlation between two observations in a series over the corresponding variable lags, i.e. Yt and Yt+h, using the following equation:

ρh=∑t=1T−h(Yt−Y¯)(Yt+h−Y¯)∑t=1T(Yt−Y¯)2

Where Yt is the current observation and Yt+h is the observation h after the current one.

2.3.4 Partial autocorrelation function (PACF)

Used to measure the degree of association between Y_t and Yt−p when the effects at other time lag Yt−1, Xt−2 ,…., Yt−p+1 are removed.

γhh={ρ1,for h=1ρh−∑j=1h−1(γh−1,j∙ρh−j)1−∑j=1h−1(γh−1,j∙ρj),for h>1

where γhj=γh−1,j−γhhγh−1,h−j,j=1,2,3,… ,h−1.

In conditions where the time series is non-stationary, data transformation has to be done by applying appropriate differencing, hence obtaining suitable values (Enders, 2015). The suitable values of p and q will be selected by observing the autocorrelation function (ACF) and partial autocorrelation function (PACF) of the time series data. The appropriate ARIMA models are selected by observing the behavior of ACF spikes and PACF based on the order identified (Hyndman & Athanasopoulos, 2018).

2.4 Model estimation

The parameters α,ψ and θ of the selected ARIMA model with specified values of p,q and d d need to be estimated. Maximum Likelihood Estimation (MLE) Method is proposed to estimate the parameters. The final model for forecasting is selected based on Information Criteria.

2.4.1 Maximum likelihood estimation (MLE)

This technique estimates ARIMA model parameters that maximize the likelihood of getting the observed series data. The method is based on providing the log-likelihood, which describes the logarithm of the probability of observed time series data of the estimated model. The function is given by:

log⁡L=−TL(ln(2π)+ln⁡σ2)−12σ2∑t=1Tεt

where T is the time of the time series data, L is likelihood function, σ2 is the constant variance and εt is the error term.

The model with the maximum log-likelihood value is considered for selection for subsequent scrutiny based on Information Criteria.

2.4.2 Information criteria

The ARIMA model subjected to forecasting is selected based on a threshold of Information Criteria. Akaike Information Criterion (AIC), Akaike Information Criterion corrected (AICc) and Bayesian Information Criteria (BIC) are assessed to specify the model with the lowest index (Hyndman & Koehler, 2006). The AIC is suitable for obtaining the appropriate ARIMA model based on the smallest values compared to other competing models. It is given by

AIC=−2⁡ln⁡L+2p

where, L is the likelihood function of time series data, and p is the number of parameters of the fitted model.

Since AIC does not consider the effect of sample size, AICc makes adjustments to allow the criteria to be used in the presence of a small sample. The following equation defines the AICc:

AICc=AIC+2p(p+1)n−p−1

The BIC is another information criterion that is useful in deciding the required model. It extends the AIC by penalizing free parameters stronger than AIC as a standalone criterion. The BIC is obtained as follows:

BIC=AIC+p(ln(n)−2)

Although all of the above criteria are provided after model estimation, this study selected AICc as a reference to choose the suitable ARIMA model. The model with the lowest AICc must be considered the best forecasting model.

2.5 Diagnostic checking

After choosing a relevant ARIMA model and estimating the corresponding parameters, model adequacy was checked by observing if the fitted model residuals were normally distributed. One of the methods of investigating whether the distribution of residuals from the fitted model is random is using the correlograms. The ACF and PACF can be used to verify if the time series displays a white noise innovation by considering that the model’s residuals should not display significant lags.

2.5.1 The Ljung-Box test

The Ljung-Box test was also used to check if the residuals followed a white noise process. The residuals were analyzed using Ljung-Box Statistic to check if the autocorrelation of the time series is significantly different from zero (Ljung & Box, 1978) and is calculated as:

the Q(k)=n(n+2)∑j=1krj2(n−j)

where n is the sample size, rj is the sample autocorrelation at lag j and k is the lag order.

The null hypothesis under this test is that the residuals are white noise, and the hypothesis is rejected if Q(k) is larger than (1−α) quantile.

2.6 Forecasting

At this stage, the idea is to use the selected ARIMA model to forecast future maize production using the past series. Forecasting should be based on the final selected model derived from the diagnostic stage. The plotted graph of the forecast and actual series will inform if the forecast is good or not before a conclusion is made on the ability of the ARIMA model to forecast future values of time series. To measure the accuracy of the prediction model, Mean Absolute Percentage Error (MAPE), Root Mean Squared Error (RMSE), Mean Absolute Scaled Error (MASE) and Mean Absolute Error (MAE) are suitable (Hyndman & Koehler, 2006). The essence here is to determine the magnitude of errors and bias, and the qualified model should register as minimum errors as possible.

MAE=1n∑t=1n|Yt−Yt^|

the RMSE=1n∑t=1nYt−Yt^2

MAPE=1n∑t=1n|Yt−Yt^Yt|*100

MASE=1n∑t=1n|Yt−Yt^|1n−1∑t=2n|Yt−Yt−1|

Where n is the sample size, Yt is the response variable at time t, Yt−1 is the response variable at time t−1, and Yt^ is the forecast value.

3. Results

3.1 Time series plot and trend analysis

The original and annual maize production data are plotted to assess the trend. The plotted graph in Figure 1 shows the volatility of the time series data. Maize production was relatively constant from 1960s to 1970s before it increased rapidly for the next ten years. The production decreased abruptly through the midpoint between 2000 and 2010 before increasing slightly from 2010 onwards. The time series under investigation is characterized by fluctuation behavior, which is a feature of many time series.

3.2 Stationarity test of time series

Table 1 shows the results of the KPSS Unit Root Test.

Table 1 shows that, the test statistic is bigger than the 1% critical value, indicating that the null hypothesis is rejected and that, the data are not stationary. The KPSS result concurs with the ACF and PACF correlograms shown in Figure 1.

3.3 Stationarity transformation

Since the results above conclude that the time series data is not stationary, differencing and natural log transformations were used for transforming the non-stationary data. This is in agreement with (Enders, 2015), that when time series data is not stationery performing differencing will make the data stationary, and therefore further analysis can be carried out.

Results in Table 2 show that the test statistic is insignificant, meaning the differenced data are stationary. This can be supported by the time series plot, as shown in Figure 2.

3.4 Model selection

Since the series tends to be stationary after transformation, the next step is identifying the order of the model's AR (p) and MA (q) to be estimated. To achieve this, ACF and PACF plots are going to be used. Based on the ACF plot in Figure 2, the order of AR is 1 because of the appearance of one significant spike compared to others. On the other hand, the MA order appears to be 1 due to one significant spike, as the PACF shown in Figure 2. Thus, the proposed model is the combination of AR (1) and MA (1), which resulted to the ARIMA (1, 1, 1) model since the order of differencing for our data is 1.

3.5 Model estimation

ARIMA (1,1,1) model was selected because of the lowest Akaike Information Criterion Corrected (AICc) of 1163.14 and the largest log-likelihood of −578.36 among other models, and it was considered the best model for forecasting. The parameters were estimated by using the maximum likelihood estimation method. The estimated ARIMA (1,1,1) model is given as:

MP=0.3413Yt−1+εt−0.8272εt−1

The established ARIMA (1, 1, 1) model is then scrutinized for suitability in forecasting maize production. The ACF plot of the residuals from the ARIMA (1,1,1) model showed that all autocorrelations were within the threshold limits and were white noise. The Liung–Box test statistic was 12.568, and the p-value was 0.1276, which is greater than 0.05 suggesting that the residuals are white noise and, thus, the model is declared fit for forecasting.

3.6 Forecasting maize production

Based on the ARIMA (1, 1, 1) model, the maize production forecast and its 95% confidence interval for the next ten (10) years are provided in Table 3. In addition, Figure 3 shows the trend of forecasted maize production of actual and forecasts, which suggests that the production has stagnated for a long time. That means the production of maize for the ten (10) years ahead is decreasing slightly with a nearly constant movement with no sign of bouncing back soon.

4. Conclusion

This paper used time series analysis to forecast maize production. Specifically, the examination used the ARIMA model to forecast the future values based on the observed series. ARIMA (1, 1, 1) model was used as it is declared fit to be used in forecasting for such data points, which were collected from 1961 to 2021. Results indicate that maize production for the next ten years (2022 - 2031) is decreasing slightly with a nearly constant movement and no sign of returning soon. The study concludes that this model performs well in forecasting maize production in Tanzania in the short run. As limitations, the present study used data on the production of maize from FAOSTAT for the period 1961 to 2021. Hence, findings may not be necessarily the same when other sources of data are used. Also, the study used maize as a single variable in the ARIMA model; the use of any additional variable would have influenced the results. Lastly, the data used are measured in hectograms per hectare and not kilograms per hectare. Further research may use more than one variable in maize forecasting, for example, the amount of rainfall recorded in mmHg for a specified period of time.

5. Policy implication

Findings from this study will enable policymakers in Tanzania and government officials to make a well-informed decision to improve maize production, which showed a slightly declining trend for the next ten years. An informed decision can include improving the delivery of new farming technologies, promoting and increasing the use of fertilizer and other complementary practices to achieve yield potential and closing gaps in technology adoption and productivity among males and females in maize production. These, together with other well-thought-out interventions, will address the declining trend of maize production.

Figures

Figure 1

Time plot and ACF and PACF plots for maize production (hg/ha) data

Figure 2

Time plot and ACF and PACF plots for differenced maize production data

Figure 3

Forecasts for the maize production in hg/ha

Table 1

The results of KPSS unit root test

Test is of type: mu with 3 lags
Value of test-statistic is: 0.9219
Critical value for a significance level of
	10pct	5pct	2.5pct	1pct
Critical values	0.347	0.463	0.574	0.739

Source(s): Created by authors

Table 2

KPSS unit root test after first differencing

Test is of type: mu with 3 lags
Value of test-statistic is: 0.0368
Critical value for a significance level of
	10pct	5pct	2.5pct	1pct
Critical values	0.347	0.463	0.574	0.739

Source(s): Created by authors

Table 3

Ten years forecast of maize production (hg/ha)

Year	Forecasts	95% confidence interval
		Lower	Upper
2022	15910.95	8541.473	23280.43
2023	15881.24	7595.026	24167.45
2024	15871.10	7196.617	24545.58
2025	15867.64	6930.918	24804.36
2026	15866.46	6707.220	25025.69
2027	15866.05	6499.859	25232.25
2028	15865.92	6300.619	25431.21
2029	15865.87	6106.655	25625.08
2030	15865.85	5916.876	25814.83
2031	15865.85	5730.785	26000.91

Source(s): Created by authors