The robustness and accuracy of Box-Jenkins ARIMA in modeling and forecasting household debt in South Africa

This paper adopted the Box-Jenkins methodology to estimate a univariate time series model. Quarterly data collected from the South African Reserve Bank covering the period 1994 to 2014 was used. The initial plot of the series revealed that household debt is explained by an irregular and non-seasonal component. Owing to the non stationarity of the series, first differencing was applied to induce stationarity. The ACFs and PACFs identified six models. Of the six identified models,ARIMA 3, 1, 0 was selected according to the standard error estimates and the information criteria. The proposed model passed all the diagnostic tests and was further used for producing ten period forecasts of household debt. The forecasted household debt rates obtained were above 75% and within confidence bounds of 95%. Insample and out-of-sampling forecasts moved together confirming the reliability of the model in forecasting household debt and vigour in predictive ability. The proposed model exhibited the best performance in terms of Max APE and Max AE and ascertained the robustness and accuracy of the BoxJenkins ARIMA in forecasting. Both a trend of the data captured and non-seasonal peaks were predicted by the model. These forecasts were proven to be realistic and a true reflection of economic reality in the country. The paper recommended a non-seasonalARIMA 3, 1, 0 be used by researchers, policy makers and decision makers of different countries to make forecasts of household debt. The South African authorities were also encouraged to use this model to produce further forecasts of the series when making long term planning.


Introduction
A sequential collection of a set of data overtime is known as time series, i.e., hourly, daily, monthly, quarterly, yearly, etc. The kind has the property that neighbouring values are correlated known as autocorrelation, (Etuk et al., 2013). Time series analysis is mainly used as pattern of change detection in statistical information over a consistent time interval. The pattern is projected to do estimation for the future. As expected, time series data is non-stationary. A stationary series has a constant mean and variance and satisfies all the properties of a white noise process such as zero mean and unit variance (Wei, 2006). An autoregressive integrated moving average ( )method is a combination of time series models and is effective in forecasting a value in a response time series as a linear combination of its own past values, past errors (also called shocks or innovations), and current and past values of other time series. This approach was first introduced by Box and Jenkins (1976) and has since become the most popular models for forecasting univariate time series data. Madsen (2008) and Boubaker (2011) are a few of the authors who popularised the model and they have also written extensively on non-seasonal models. The methodaffords a broad set of tools for building a univariate time series model. The procedure for this method includes model identification, parameter estimation and selection, and forecasting. This paper explores the application of time series analysis to forecast household debt in South Africa (SA). To be precise, the paper is aimed at selecting the model that may be suitable in generating the forecasts of the series. The choice of this model is due to its innovativeness and wide use in recent studies. This model is suited to perform duties as highlighted in the previous paragraph. Additionally, the frameworkhas has been found to be reliable in terms of its performance. This framework has advantage over others as it provides an analyst an opportunity to make a choice between candidate models and can effectively handle larger number of time series, the duty which multivariate analyses fails to perform. It also provides a room to critically scrutinise candidate models and decide on the one that is best suited for the data. The framework is also uncomplicated and gives the results that allows the researcher to make a reflection about the data based on its past, present and future state.
The emphasis of is basically on forecast performance and advocates more focus on minimising out-of-sample forecast errors than on maximising in-sample 'goodness of fit' (Meyler, Kenny and Quinn, 1998). Consequently, isdefiantly one of 'model mining' with the aim of optimising forecast performance. There is no prior knowledge or impression about the model to be fitted or chosen. The onus to decide on whether the model is adequate or not is dependent on the measures used and its performance to provide forecasts. The framework is appealing to this study since only one variable (univariate analysis), household debt is analysed. Stockton and Glassman (1987) and Litterman (1986) advocate for this framework as it frequently outperform more sophisticated structural models in terms of short-run forecasting ability. The following authors among others as also cited by Suhartono (2011) used these models to effectively perform forecasting of time series variables (economic and financial data); Haswell et al. (2009) applied this model for forecasting soil dryness index, Modarres (2007) and Abebe and Foerch (2008) used it for drought forecasting, Briet et al. (2008) for short term malaria prediction, Momani (2009) for forecasting rainfall, Wagner (2010) for forecasting daily demand in cash supply chains. Other applications of these models can be found in Vu and Turner (2005;2006), Chu (2008;), Chang et al. (2009), de Oliveira Santos (2009), Coshall (2009), Song et al. (2009.
Forecasting household debt of SA is of special importance for policy makers. Without proper and reliable forecasts, it becomes difficult for policy makers to come up with feasible plans. With accurate forecasts generated, it may be easy for the responsible authorities to formulate and device appropriate strategies to deal with household debt in the country. Also with proper modelling technique applied, improved forecasting accuracy of household debt in the country may be attained. Furthermore, this paper adds to the literature by evidently studying trends in household debt using the Box-Jenkins approach and for assessing the effectiveness of this method in forecasting. Additionally, the paper contributes to the research information on household debt in SA. Both the application and theory that is accessible to a various students, practitioners and researchers is presented in this paper. The remainder of this paper is organized as follows; Section 2 defines the data used and Section 3 describes the procedure used for data analysis. In Section 4, empirical results are discussed and Section 5 provides conclusions and some recommendations for future studies policy modification.

Data:
A twenty year historical household debt to disposable income data which is seasonally adjusted at current pricesin SA are collected from the first quarter of1994to the first quarter of 2014.A total of 80 observations are used. The data is sourced from the South African Reserve Bank website (www.sarb.co.za). This data was subjected to log transformation prior to the analysis. The Statistical Software Analysis (SAS) version 9.3 was used to execute the analysis of data.
Household debt to disposable income: Clark and Daniel (2006) define this variable as a ratio that measures the percentage of the average households' disposable income used to repay debt. Increased households are associated with decreased interest rates which lure people to entering into mortgage debt, increased inflation rates, unemployment rates and other. Additionally, the ratio in household debt increases as the amount of disposable income decreases. However, this study aims at building a univariate model that will be used to forecast this variable on the basis of its past and present values.

Theoretical framework
This section describes the methodological procedure used for data analysis. The methodology involves the four stages according to Box-Jenkins. Identifying ARMA Models: This is the most crucial step in time series analysis. There are different ways to identify the fitted seasonal and/or non-seasonal ARIMA model. A set of autocorrelation function (ACF) and partial ACF may be visually inspected to decide on the most appropriate model for the data (Faraway and Chatfield, 1998;Kihoro, Otieno, Wafula (2004) and Hipel and McLeod, 1994). These are statistical measures which reflect how the observations are related to each other in a time series. If the objective of the study is modelling and forecasting, it is often useful to plot the ACF and PACF against consecutive time lags so as to determine the order of non-seasonal and/or seasonal autoregressive and moving average terms. An appropriate ARIMA and/or SARIMA model is considered in identifying a tentative model. This procedure refers to the methodology in identifying necessary transformations to the data, the decision to include the deterministic parameter and proper orders of p and q and/or P and Q for the said model (Hu et al., 2010). The following are important steps in the identification the process as outlined in Hu (ibid): Step 1: Provide plot of the series and choose proper transformations. The plot will reveal features of a trend, seasonality, outliers, non-constant phenomena, etc. This information provides a basis for postulating possible transformations to the data at hand.
Step 2: Compute and examine the ACF and PACF of the original series. This helps in further confirming a necessary degree of differencing to induce stationarity to the data. Differencing is needed when the ACF decays slowly and the PACF cuts after lag 1.
Unit root tests provided by Dickey and Fuller (1979;1981) may be used in this scenario. The following ARIMA (p, d, q) models are recommended for economic data. It is a modified DF test also known as ADF and is estimated as follows: Where∇ is the first difference operator; t is the time drift; k represents the number of lags used and  is the error term; ′ and ′s are the model bounds. The correct value for k is determined using the Akaike and the Schwarz-Bayesian information criteria. The ADF test includes a constant and time trend. For the decision rule, assuming that the series, − 1 follows the AR (p) process, Hamilton (1994) shows that the rejection or acceptance of the null hypothesis of a unit root is based on running the regression; Where −1 = −1 − − −1 for j = 0, 1, 2, . . . , − 1and is a white noise process. The ADF test statistic is given as; Where − is the test statistic of − , ( ) is the standard error of − . The null hypothesis of a unit root 0 : − is rejected if [3] is less than the appropriate critical value at some level of significance. Alternatively the test statistic rejects the null hypothesis if the corresponding probability value exceeds the level of significance. Table 1 shows characteristics of theoretical ACF and PACF for stationary process. To obtain optimal results, the patterns in the ACF and PACF should be matched with these theoretical patterns. Tails off as exponential decay or damped sine wave Cuts off after lag p MA (q) Cuts off after lag q Tails off as exponential decay or damped sine wave ARMA (p,q) Tails off after lag (q-p) Tails off after lag (p-q) If the visual examination of the ACF and PACF does not provide appropriate conclusions, other widely used measures such as the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) proposed by Faraway and Chatfield (1998) may be useful. Mathematical formulae of these measures are given as: Where the likelihood function is evaluated at a minimum likelihood estimates and T is defined as a sample size (Tsay, 2010). For a Gaussian model, AIC reduces to: Where ℓ 2 denotes the maximum likelihood estimate of 2 also known as the variance of . Tsay (2005) explains the first term of this criterion as a measure of the goodness of fit of the model to the data. The second term he defines it as the penalty function of the criterion. The author explained that this name is given to this criterion as it penalises the candidate model by the number of parameters used. The BIC is presented as: [6] The penalty for each parameter used is 2 for AIC and for BIC. On this basis, BIC has a tendency to select the model with the least lag more specifically when the sample size is moderate or large. In this study, the both the AIC and BIC are computed for ℓ = 0, … , , where P is a pre-specified positive integer. The optimal model order is chosen by the number of model parameters, which minimizes either AIC or BIC.
Model estimation and selection: Suppose a non-seasonal ARIMA (p, d, q) model for time series = 1 , 2 , … , is of the form: Where {at} is a white noise series and p and q are non-negative integers. The AR and MA models are special cases of the ARMA (p, q) model. Using the backshift operator, the model can be written as: Tsay (2010) suggested no common factors between the AR and MA polynomials; otherwise the order (p, q) of the model must be in a reduced form. He also suggested that the AR polynomial should introduce the characteristic equation of an ARMA model, a duty performed by the AR model. Weak stationarity of an ARMA model is satisfied only if all the absolute values of the solutions of the characteristic equation are less than one. If this is the case, the unconditional mean of the model according to Tsay (2010) is: From [8], the AR polynomial is the same as ( ) and the MA polynomial becomes ( ). Invertibility condition is satisfied if the roots of = 0 lie outside the unit circle. To satisfy stationarity, it is required that = 0 also lies outside the unit circle.Non-stationarity problem can be dealt with by reducing the series to stationary time series. This can be done by taking a proper degree of differencing. Introducing a differencing term in the model will result in ARIMA model, which is used to describe various homogeneous non stationary series. Differencing 1 − is applied to [8] and the resulting stationary ( , ) model becomes: Where stationary AR operator = 1 − 1 − ⋯ − and the invertible MA operator = 1 − 1 − ⋯ − share no common factors. According to Wei (2006), the parameter 0 plays very different roles for = 0 and > 0.When = 0, the original process is stationary. When ≥ 0, however, θ0is called a deterministic trend. This term may be omitted from the model unless it is really needed (Hu et al., 2010). The model in [10] is called ( , , ).
Conditional least squares method is used to estimate parameters of the model. The sample size plays an important role in the estimation. At least 50 observation or more is needed according to Box and Jenkins (1976).The error term of the model is assumed to follow a white noise process. Selection of the model that best fits the data is done with the AIC and the BIC as shown in equations [4] and [6]. The model with least values of the statistics is considered the best and will be used to generate the forecasts.

Diagnostic checking:
The overall fit of the model is evaluated by assessing the residuals. According to Box and Jenkins (1976), the residuals should follow a white noise process with zero mean and unit standard deviation. A graphical representation of the standardized residuals should look random and homoscedastic. Again an ACF is used to evaluate the assumption of randomness, or alternatively the Ljung-Box test (1978) may be used. This test is defined as: The statistic is asymptotically distributed as 2 with m degrees of freedom (Vogelvang, 2005). Where ρ is called the sample ACF of (Tsay, 2005), T is the number of observations, k is the highest order for autocorrelation for which to test and 2 is the ℎ autocorrelation. The square of this test used to check the presence of autocorrelation in the model is as follows: * > ; − − − − . [12] The null hypothesis is rejected if the observed ≤ 5%significance level or alternatively if the observed statistic is in excess of the critical value. This implies that autocorrelation exists up to order k.The assumption of normality of the residuals is evaluated with a histogram and QQ-plots.
Forecasting using an ARMA Model: Forecasts of an ( , ) model have characteristics similar to those of an AR (p) model after adjusting for the impacts of the MA component on the lower horizon forecasts. Tsay (2010) denote the forecast origin by h and the available information by ℎ . The author suggested the 1-step ahead forecast of ℎ +1 from ARMA model as follows: The associated forecast error is calculated as: The variance of 1-step ahead forecast error is: = . [15] For the ℓ − step ahead forecast, we have: Note that ℓ − = +ℓ− if ℓ − ≤ and ℓ − = ifℓ − > 0 and ℓ − = ℎ +ℓ− if ℓ − ≤ . As a result, the multistep ahead forecasts of an ARMA model can be computed recursively (Wei, 2006). The associated standard error is computed as: Performance Comparison: Once the model has been fitted, subjected to diagnostic checks and used to produce forecasts, it is evaluated with forecast fit measures. This study uses the maximum absolute percentage error (MaxAPE) and the maximum absolute error (MaxAE) to evaluate the amount of this forecast error. These measures are computed using the proposed model both for the original and forecasted series. The model with least forecast errors is considered accurate.

Empirical Analysis
The analysis of data was informed by the methodology reviewed and the results are provided in this section. The initial step in time series analysis is to plot the series against time so as to help identify the properties of a model to be fitted. The results are presented in Figure 1. The effects of this crisis spilled over to many economies and resulted in increased household debt levels in many countries of differing financial sophistication. SA was one of the countries which suffered these effects. An insignificant decline is seen from the last quarter of 2009. This could be as a result of national credit regulations which were legislated to help reduce debts in the country. SA also experienced a recession during this period. Due to the increasing irregular component displayed by the figure, it is concluded that the series is nonstationary process. Furthermore, no seasonal variation is exhibited by the figure. A substantial number of the data points are far away from the mean (i.e. the line runs horizontally in a maroon line) confirming that the series is not stationary in the mean. Presented next are the ACFs and PACFs of the log transformed series. The purpose of these results is to examine if the series follows a white noise process or not. This will also help in determining the level of differencing needed. The results are presented in Figure 2.

Figure 2: ACF and PACF of log transformed household debt
The figure shows the ACF (located on the upper right corner of the figure) dying down slowly and all the lags up to about seven are statistically different from zero. All the lags lie outside a set confidence bounds of 95%. In addition, the PACF (lower left corner) cuts off dramatically after lag one and all other lags appear to be insignificant. The results are in accordance with theoretical guidelines as summarised in Table 1. This provides a good basis to conclude that the series is not stationary according to Box and Jenkins (1976). To be able to apply Box-Jenkins methodology, differencing was applied to further transform the series to induce stationarity and the results are displayed in Figure 3.

Figure 3: ACF and PACF of first order differencing of household debt
Inspection of Figure 3 shows that the series became fairly stationary at first difference. The observations in the plot of first differencing depict an irregular movement but revert to their mean value with an approximately constant variability. Only lags 1, 4 and 5 of the ACF and lags 1 and 3 of the PACF are outside the constraints after first order differencing. These lags are confirmed to be the only lags which are statistically different from zero. This fairly implies that the transformed series became stationary at first difference and allows for the analysis to be continued. The mean and variances of the series are constant and therefore follow a white noise process (Hamilton, 1994). No sign of seasonal difference is revealed in the series even after data transformation. The ADF test in Table 2 was computed to confirm the results displayed by Figure 3. Only the observed probabilities of ADF associated with zero mean are significant at 10% level of significance but insignificant at all levels in the presence of a single mean and a trend. Therefore the ARIMA models will be fitted without a mean and trend components restrictions. The results also just as revealed on Figure 3 confirm stationarity of the series after first difference. It also follows from Table 3 that the ACF follows a white noise process with a mean equal to zero and a unit variance. Therefore the properties of stationarity are satisfied and allows for the forecasting model to be estimated with first order differencing and no mean and trend components restrictions imposed. The results of this analysis are in accordance with Box and Jenkins (1976), Dickey and Fuller (1979;1981) and Hu et al. (2010). Since the correlograms given in Figure 3 are stationary, they are a good basis to be further used to identify the competing models. It is important to note that the patterns of ACFs and PACFs of AR(p) and MA (q) processes are different and may as a result identify different models. The following models are identified by trial and error from the figure; 1, 1, 0 ; 3, 1, 0 ; 0, 1, 1 ; 0, 1, 3 ; 1, 1, 1 ; 3, 1, 3 .
The parameters of the identified models are estimated using the conditional least squares method. The models are fitted and the one that fits the data well according to the tests are chosen to produce the forecasts of household debt in the country. The results of these models are summarised in Table 4.  Table 4 above, parameter estimates of 3, 1, 0 appears to be more significant according to the corresponding approximated probabilities than those of other models. Further note that the observed t-value of this model parameter is larger that of other models. 1, 1, 0 and 0, 1, 1 also have significant parameters and larger t-values. Other models have mixture of parameters (significant and insignificant). The standard error of 3, 1, 0 is less than those of competing models followed by that of 3, 1, 3 .The AIC and SBC also choose 3, 1, 0 as a model that best fits the data (Faraway and Chatfield, 1998). Therefore this model meets all the requirements and will further be used for forecasting of household debt. Mathematically, this model according to [7] with no constant, mean and trend according to the results in Table2 is given as:
Before forecasts can be produced, the model was subjected to the diagnostic checking in order to ensure its suitability for forecasting and the results are given as Figure 4.

Figure 4: Residual normality diagnostics
Firstly, the model residuals were checked and they were found to be approximated by a normal distribution as the histogram and Q-Q plots reveals. This validates the normality of the residuals assumptions and confirms 3, 1, 0 as a sufficient candidate model.

Figure 5: Residual correlation diagnostics
Next, the residuals of the model were assessed for randomness. All the residuals were also not significant as noted in Figure 5. They were all less than 5% significance level. The ACFs and PACFs have all lags within the 95% confidence bounds except for lag one of the ACF. This implies that the residuals of the selected model are random and homoscedastic. The model is therefore adequate and can be used for forecasting household debt. The results are also in accordance with Ljung-Box (1978) as described by [12]. Forecasts of household debt in SA were generated using 3, 1, 0 model. The associated confidence intervals are provided and the results are given as a summary of Figure 6. Note that the forecasts in Figure 6 are based on the 95% confidence limits. Past values of the series (shown as circles) are in par with the forecasted values (in a solid line). This confirms the reliability of the model ( 3, 1, 0 ) in forecasting household debt and its vigour in predictive ability. Both a trend of the data is captured and non-seasonal peaks are predicted by the plot. Forecasts were produced using the selected model for the next two years and the remaining three quarters of 2014. Extrapolated values reveal a wide confidence interval and somehow sluggishly rising trend. A two-period ahead forecast was yielded using historical data of 1994 first quarter to 2014 first quarter. These short-term period forecasts of household debt in SA are between 2014 second quarter to 2016first quarter and shows that about 75% or more of debts are expected in the next two years. These estimates are still high according to South African standards. In the short run, it is important for the authorities in SA to come up with policies that may be useful in addressing these high expected trends. Urgent attention needs to be paid to this problem and possible solutions must be sought sooner to try preventing another recession.
, , models have more advantage when used in forecasting in that they simultaneously provide the results of the forecast based on the past and future values of the series as depicted by Figure 6. The results of the two series do not drift too far from each other explaining the rigidity of the model. The values of the original and forecasted data were compared to further assess the robustness and validity of 3, 1, 0 model. The results are presented in Table 5. It is important to note that the MaxAPE and the MaxAE for the model of forecasted data are less than those of the original data. MaxAPE represents the largest forecasted error expressed as a percentage. This measure is as being useful for imagining a worst-case scenario for the forecasts. Ljung-Box tests for these models are both greater than 0.05 with one for forecasted model even bigger. This proves that the residuals of the second model are random and heteroscedasticity. Given these results, 3, 1, 0 is confirmed to be robust and valid in modelling and forecasting household debt for the selected period in SA.

Conclusion
In the main, the paper used Box-Jenkins methodology to estimate the model used for forecasting household debt in SA. A time series data collected for the first quarter of 1994 to the first quarter of 2014 from the SARB was used. Time series plot revealed that household debt are explained by an irregular and non-seasonal component. Owing to the non stationarity of the series, first differencing was applied to the series to induce stationarity. The ACFs and PACFs identified six models. These models were compared and 3, 1, 0 was chosen according to standard error estimates and the information criteria. The model was subjected to diagnostic checks which confirmed its suitability to produce forecasts. The forecasts produced were above 75% with confidence bounds of 95%. The forecasted values from the proposed model are much more realistic and are a true reflection of economic reality in the country. To evaluate forecasting accuracy of this model, past and future values of the series were compared. The MaxAPE and MaxAE proved the robustness and accuracy of the proposed model. The paper therefore recommends a non-seasonal 3, 1, 0 be used by researchers, policy makers and decision makers of different countries to make forecasts of household debt. The South African authorities are also encouraged to use this model to produce further forecasts when making long term planning. It cannot be emphasised further that preliminary analysis of time series data is very important as it guards against choosing irrelevant model and also endeavors to reveal hidden characteristics of the series.The use of quarterly data is encouraged to help avoid violation of non-normality of residuals and autocorrelations.