Modelling and Forecasting Volatility of Value Added Tax Revenue in Kenya

Taxation is one of the means by which governments finance their expenditure by imposing charges on citizens and corporate entities. Kenya Revenue Authority (KRA) is the agency responsible for the assessment, collection and accounting for of all revenues that are due to government. Volatile government revenue is a challenge for fiscal policy makers since it creates risks to government service provision and can make planning difficult, as revenue falls short of expenditure needs both frequently and unexpectedly. The main objective of this study was to model and forecast the volatility of VAT revenue collected in Kenya as well as computing its value at risk and the expected shortfall. Secondary data on daily VAT revenue collections for a period of 3 years was analyzed. The first step was to model the mean equation of the return series using the ARIMA model and ARIMA(3,0,3) was identified to be the most suitable since it had the least values of AIC and BIC. The Lagrange Multiplier test confirmed the presence of ARCH effects using the residuals of the mean equation. A number of heteroscedastic models were fitted and the TGARCH family (ARIMA(3,0,3)/TGARCH(1,2)) was preferred to fit the volatility of the returns. One step ahead forecasting of volatility of the returns was done using the model which gave a value of 7.212. Estimation of value at risk and expected shortfall involved use of POT method by fitting a GPD function to the return data series. The first step was determination of threshold by use of MRL plot and later fitting a GPD function to the return data series using the threshold. The shape, location and scale parameters were estimated using MLE and they were later used to compute the VaR loss and ES at 95% and 99% confidence intervals. The VaR at 95% and 99% was 1.45% and 1.49% respectively while the ES at both the intervals was 0.04% and 0.1% respectively. This study concluded that volatility is persistent in the daily VAT revenue collections and it can easily be modelled using conditional heteroscedastic models.


Introduction
Revenue volatility can be defined as the extent of revenue fluctuation over the course of the business cycle. It is a serious concern particularly for state policymakers and fiscal administrators operating within the context of balanced budget requirements. It makes it hard to make accurate forecasts for future revenues and the establishment of longterm fiscal plans for the stable operation of public programs and services. With the global economy being liberalized and more tightly interwoven, the tax environment of governments has been increasingly volatile and uncertain over the past years, and as a result, fiscal planning and management have become more challenging particularly for state governments operating under the institutional constraints of balanced budget requirements.
The important implication of revenue volatility for state finance is that in the absence of adequate fiscal reserves, it is hard for states with volatile revenue bases to avoid massive spending cuts and tax hikes in times of economic crisis when governments counter cyclical fiscal actions are needed more than ever. Another implication is that such pro-cyclical austerity measures affect real economies, reducing households and businesses' propensity to consume and consequently creating the vicious circle of economic recession.
Capturing revenue volatility is very important. Looking at the perspective of a prospective planning, if we desire to make decisions about the structure our revenue sources on which we will rely on to provide sufficient revenue for all essential government services, we should consider both the mean rate of growth of the revenue sources and the volatility as well. Also, revenue volatility can be used as an independent variable in developing an understanding about past decisions such as revenue structure, the use of rainy day funds, or state borrowing trends among others.
A study was carried out on the importance of considering tax portfolios and economic conditions in determining the growth rate and volatility of state tax revenues. Using 1989 -2009 state data and simple graphical constructs, various comparative analysis of the long-run growth rate and shortrun volatility in percent changes of state economies and tax revenues and state tax portfolios were constructed. In the analysis, it was found that wide variations in these respects exist among states. Other findings suggested that in the short run, states cannot alter the underlying structure of the economy but can mitigate the impacts of the business cycle on their fiscal conditions by making changes to their tax portfolios [1].
Tax analysis and revenue forecasting are very essential in providing directions as to enhancing tax revenue, improving equity and efficiency of taxes, promoting investments and consequently economic growth with a view of increasing national income. In addition, it is important in monitoring processes and budget planning. It defines the resource envelope forming the basis for effective medium-term and long-term planning [2].
Revenue forecasting is an essential part of budgeting in the public sector. Therefore, it is necessary for a government to forecast the revenue it collects for planning purposes. Budgetary uncertainties have led to increased reliance on economic and revenue forecasting by state governments in recent years. Because of the magnitude of the fiscal problems facing many states, forecasting has assumed a more central role in the policy making process. As a result, revenue forecasts are closely examined and accuracy is essential for planning purposes. To improve accuracy, there is need for analysts to assemble as much information about their respective state economies as possible, including formal and informal consideration of alternative forecasts [2].
In Kenya, the tax revenue analysis and forecast is produced by a single agency, which is the Ministry of Finance in conjunction with the Kenya Institute for Public Policy Research and Analysis (KIPPRA). The ministry formulates the government budget by undertaking revenue projections and at the same time prepares the monthly/annual revenue targets for relevant revenue collection agencies in the country. VAT in Kenya is usually a tax on consumption, therefore its volatile depending on the economic seasons. The purpose of this project was to make use of heteroscedastic models to model and forecast the volatility of the VAT data series. The study also looked at the risks associated with the VAT by estimating its Value at Risk (VaR) and the Expected Shortfall (ES).

Methodology
Time series analysis is the procedure of fitting a time series to a proper model. It comprises methods that attempt to understand the nature of the data series and its characteristics which include trend, seasonality, outliers and cycle. When modeling the volatility of a return series, a mean equation is needed for testing the arch effects. ARMA or ARIMA models can be used to model the mean equation

Autoregressive Integrated Moving Average (ARIMA)
From application view point ARMA models are not the best to properly describe non-stationary time series, which are frequently encountered in practice. For this reason, the ARIMA model was proposed, which is an extension of an ARMA model to include the case of non-stationary data as well. In ARIMA models a non-stationary data is made stationary by applying finite differencing of the data points [3].

Heteroscedastic Models
Volatility modeling provides an approach of calculating Value at Risk (VaR) of a financial position in risk management. Building a volatility model for a return series involves the following steps: i Specifying of a mean equation by testing for serial dependence in the data series. If need be, build an econometric model for the return series so as to be used to remove any linear dependence. ii Use of the squared residuals of the mean equation to test for ARCH effects. iii Specifying of a volatility model if ARCH effects are statistically significant and performing a joint estimation of the mean and volatility equations. iv Checking the fitted model carefully and refining it if necessary.

Testing for ARCH Effects
Let a = y − u be the residuals of the mean equation. The squared series a is used to check for conditional heteroscedasticity (ARCH effect). This can be done using the Lagrange multiplier test which is equivalent to the usual F statistic for testing α (i = 1,2, … m) in the linear regression. a = α + α a + ⋯ + α ! a ! + e , t = m + 1, ⋯ , T (2) Where; i e denotes the error term. ii m is a pre-specified positive integer. iii T is the sample size. This gives us Which is asymptotically distributed as chi-squared distribution with m degrees of freedom under the null hypothesis. This test is known as Lagrange multiplier test.

The ARCH Model
The first model that provides a systematic framework for volatility modeling is the ARCH model of [4]. The basic idea of ARCH models is that: i The mean corrected asset return a is serially uncorrelated, but dependent. ii The dependence of a can be described by a simple quadratic function of its lagged values. Specifically, an ARCH (p) model assumes that Where i ϵ is a sequence of independent and identically distributed (iid) random variables with mean zero and variance 1. ii α > 0 and α ≥ 0 for i>0.

The GARCH Model
Although the ARCH model is simple, it requires many parameters to adequately describe the volatility process of an asset return. A useful extension known as the generalized ARCH (GARCH) model was therefore proposed [5]. For a log return series y , it's assumed that the mean equation of the process can be adequately described by an ARMA model. Let a = y − u be the mean corrected log return. Then a follows a GARCH (m,s) model if Where i ϵ is a sequence of independent and identically distributed (iid) random variables with mean zero and variance 1. ii α > 0, α ≥ 0 β 3 ≥ 0.

The TGARCH Model
This is one of the models used to handle leverage effect [6]. A TGARC H (m, s) model assumes the form σ = α + ∑ (α + γ N ) 4 2 a + ∑ β 3 σ 3 ! 32 iii α , β 3 and γ are non-negative parameters satisfying conditions similar to those of GARCH models. From the model, it is seen that a positive a B ; contributes α a B ; to σ whereas negative a B ; has a larger impact (α + γ )a B ; with γ > 0 . The model uses zero as its threshold to separate the impacts of past shocks. This model is preferred over the others because it allows the asymmetric effects between positive and negative asset returns [7].
Model Selection Criteria Both AIC and BIC were used in the model selection. These can be expressed as; Where; i L is the maximized value of the likelihood function ii. N is the number of recorded measurements. ii k is the number of estimated parameters. Existing studies shows that AIC is not consistent and has the tendency to choose models which are over-parameterized.

Forecast Evaluation
The mean absolute error (MAE) and the mean square error (MSE) were used as measures for evaluating forecasting performance. The MAE and MSE for n step ahead forecast are denoted as follows Where; i r BJS is the return over horizon n steps ahead at current time t. ii r̅ is the mean of return. iii h O (n) is the forecasted conditional variance over horizon n steps ahead at current time t.

Value at Risk (VaR) and Expected Shortfall (ES)
One way of modeling VaR and ES by is applying extreme distribution based on methods such as Extreme Value Theory (EVT). EVT refers to branch of statistics which normally deals with the extreme deviations from the mean of a probability distribution. Results in EVT can be obtained either using Block Maxima model (BMM) or Peak Over Threshold model (POT). POT method models a distribution of excess over a given threshold. EVT shows that the limiting distribution of excess is a Generalized Pareto Distribution of GPD [8].
The POT method uses available data more efficiently compared to BMM. In POT, all the data which exceeds a particular threshold level is used while in BMM only the maximum from a block length is retained for distribution estimation.

GPD and EVT
The generalized Pareto distribution (GPD) is a family of continuous probability distributions. It is often used to model the tails of another distribution. It is specified by three parameters: location µ, scale β, and shape parameter ξ.

Peaks Over Threshold (POT)
To fit a GPD to our data set, adoption of the POT method that focuses on the distribution of excess above some high threshold is made. For y − u ≥ 0 the excess distribution function can be rewritten as and hence deduce the following reverse expression which allows us to apply the POT method. There are two steps in applying the POT method. First is to choose an appropriate threshold and then fit the GPD function to data [9].

Mean Excess Function and Threshold Selection
For a random variable Y, the mean excess function is defined as i.e. the mean excess over a threshold u. If the underlying distribution Y > u follows a GPD, then the corresponding mean excess is provided ξ < 1. From the equation it's clearly seen that the mean excess function must be linear in u. More precisely, Y > u follows a GPD if and only if the mean excess function is linear in u. This gives us a way of selecting an appropriate threshold.

Model Validation
Quantile plots can be used to assess the quality of a fitted generalized Pareto model. If GPD is a reasonable fit for the excess above u, then the Q-Q plot should depict points that are approximately linear. Furthermore, the goodness-of-fit of GPD can be confirmed by utilizing the excess distribution plot and plot of the tail of underlying distribution. For a good fit, the excess should lie close to the theoretical curves. Lastly, a scatter plot of residuals should not depict any visible pattern to indicate independence of the excess.

VaR
For a random variable Y (usually the return in some risky financial instrument) with distribution function F over a specified time period, the VaR, for a given probability p, can be defined as the p th quantile of F, i.e., where F is the quantile function. VaR is a common measure of extreme risks and GPD is used to approximate this measure. In particular, using Equation of POT the following is obtained.
Where i β and ξ are the maximum likelihood estimates of the GPD parameters. ii n is the total sample size. iii Nu is the amount of observations above the chosen threshold.

ES
ES was proposed as a better measure of risk, which is subadditive and also informs us about the likely magnitude of excess [10]. In contrast to VaR, ES is used measure the riskiness of an instrument by considering both the size and likelihood of losses above a particular threshold. ES gives the expected size of return that exceeds VaR.

ES s = E(Y |Y > VaR
(20) And equivalently where the second term above represent the mean of the excess distribution. Proceeding as before, if the threshold VaR is sufficiently large then the expected shortfall can be computed as follows.

Results and Discussion
Secondary data of the daily VAT revenue collections from January 2016 to March 2019 was used for this research.  Figure 1 shows volatility clustering of the VAT revenue collections. It is noticeable that periods of low volatility are followed by periods of low volatility and periods of high volatility followed by periods of high volatility as well. The periods of high volatility are few days prior to the due date of payment. During these days, payments tend to increase with the highest payment being made on the due date of the month. After the due dates, the amounts of payments reduces drastically and the collections made in the proceeding days are usually low.
These VAT payments were converted in to their log returns and a time plot of the return series is as shown in Figure 2.

Stationary Test
The aim of this test was to check the presence or absence of unit root, where the null hypothesis is stated by default as the presence of unit root and the alternative hypothesis as absence of a unit root. Augmented Dickey Fuller (ADF) was used to test if data was stationary and produced the following results. Since the p-value is less than 0.05 significance level, we reject the null hypothesis and hence conclude that the returns data series is stationary. The series was later modeled using ARIMA models and checked for the arch effects.

Model Selection
Various ARIMA models were fitted to the differenced series and selection was made using the AIC and BIC values. The model with the lowest values of AIC and BIC is preferred. From table 2, the ARIMA (3,0,3) is preferred to be the best model since it has the least values of AIC and BIC.

Testing for ARCH Effects
This was done using the residuals of the mean equation. Using the squared residuals, Lagrange Multiplier (LM) test was used to test for the arch effects and the results obtained are as follows. The null hypothesis is that there is no arch effects. From table 3, the p-value obtained was less than 0.05. Therefore, the null hypothesis is rejected and a conclusion is made that there is presence of arch effects.

Model Selection
After confirming the presence of arch effects in our return data series, various ARIMA (3,0,3)/GARCH (p,q) as well as ARIMA(3,0,3)/EGARCH(p,q) and ARIMA(3,0,3)/TGARCH (p,q) models were considered for modelling the VAT returns and selection was made based on the values of the AIC and BIC.  From table 4 above, it is evident that the best model is ARIMA(3,3)/TGARCH(1,2) since it had the lowest value of AIC and BIC using student-t distribution.

Parameter Estimation
Using maximum likelihood estimation, the parameters of our chosen model were estimated and the results obtained are as shown in table 5. From table 5, the p-values of most of the parameters is less than 0.05 which indicates that they are statistically significant while few are greater than 0.05. The value of γ is positive which implies that negative events will have more impact on the volatility of the returns than positive events.

Forecasting
Using the TGARCH model, one step ahead forecast of the VAT returns was done and the results are as shown in table 6

Estimation of Value at Risk (VaR) and Expected Shortfall (ES)
This was done using the extreme value theory by using GPD function.

Determination of Threshold
Before fitting the GPD function to the data, it is first necessary to choose a threshold. Mean Residuals life (MRL) plot of the return data series was used to estimate the threshold as shown in figure 3 Since the idea is to select a threshold whereby the graph is linear, a threshold of 0.1 appears to be a reasonable value since a reasonably straight line could be placed within the uncertainty bounds from this point up.

Estimation of GPD Parameters
The next step was to estimate the scale parameter σ and the shape parameter ξ by fitting a GPD function to the data series. Maximum likelihood estimate of the parameters was obtained and the results are as shown in table 7. These values obtained were used to estimate the values of VaR and ES losses.

Model Adequacy
To assess the quality of the fitted model, various plots were used. These are excess distribution plot, tail of underlying distribution plot, scatter plot of residuals and Q-Q plot of residuals. From figure 4, the top plots from the left are excess distribution plot and plot of tail underlying distribution respectively. These plots show that the excess lie close to the theoretical curves. Scatter plot of the residuals does not depict any visible pattern in the residuals which indicate independence of the excess. Again the Q-Q plot depicts points which approximately linear. Therefore, all of the above plots confirmed the goodness of fit of our model.

VaR Computation
After confirming the goodness of fit of our model, the VaR was computed making use of the GPD parameters that were obtained through MLE. VaR was computed at 95% and 99% confidence interval and the results obtained were as shown in table 8.  Table 8 shows that we are 95% and 99% confident that the worst loss of daily VAT collection will not exceed 1.45% and 1.49% respectively.

Estimation of Expected Shortfall
The ES was estimated in a similar way as the VaR by fitting GPD function to the data series. The scale parameter, shape parameter and the location parameter obtained by MLE from the GPD function together with the VaR estimated were used to compute the ES at different confidence intervals and the results obtained as in the table below Since ES is normally the expected loss of a return that exceeds VaR, we can conclude from Table 9 that if the VaR loss of our return becomes 1.45% at 95% confidence interval, then the worst loss on daily VAT collection will not exceed beyond 0.038%. Similarly if the VaR loss at 99% confidence interval is 1.49%, then the worst loss on daily VAT collection will not exceed beyond 0.096%.

Conclusions
This research study investigated the best model to capture the volatility of daily VAT revenue collected in Kenya from January 2016 to March 2019. Different statistical tests were carried out on the data so as to determine the best model to fit the returns. The analysis showed presence of components such as trend, seasonal variation and random changes.
Because of these, it was possible to model the data using time series ARIMA(3,0,3) model. The results showed that volatility was persistence in the VAT returns. ARIMA(3,0,3)/TGARCH(1,2) was identified as the best model to forecast the volatility of the revenue returns. Going further, the study estimated the market risks associated with the revenue returns. These were Value at Risk and the Expected Shortfall. Using Peak Over Threshold method of the Extreme Value Theory, various parameters that are required were computed by fitting a Generalised Pareto Distribution function in the data. These were then used to estimate the value of the risks (losses) at various confidence intervals.

Recommendations
The general recommendation goes to the future researchers to examine the effectiveness of other heteroscedastic models like EGARCH, IGARCH and GARCH-M in modeling and forecasting of revenue volatility. Also, other methods of risk measures like use of block maxima model (BMM) and GARCH model can also be looked at and a comparison made between the different methods. Future research can also look at other tax revenues and expenditure as well.