Forecasting tax data for next 12 months
1.Analysis objective
This analysis has an objective to forecast revenue collection for the next 12 month from the last available data in 2021. Based on monthly Pakistan Tax Collection data, these time series cumulate from 2001 to the latest 2020 available data. __ This is still under development and in detail video on how to make this analysis and creating a file like this one will be uploade on Data analysis and econometric
Main purpose of this analysis is to develop a tutorial for students who use time series modeling for forecasting. However, FBR officials can also use this analysis for forecasting and develop scenario for themselves. Our time series is a conventional univariate time series with date and monthly tax variables which we will list in the following chunk. The step applied in this analysis can be used in any forecast analysis.
2. Models used: explanation
This analysis is statistical forecasting (not economic forecasting) based on three classic forecasting models(Seasonal naïve, Exponential smoothing, ARIMA) and experienced one neural network model Seasonal naïve model.
I will start from very basics and will increase complexity gradually. Let me explain each model briefly.
As it will be observed in the following graphs that data have strong seasonality. Therefore, the first more basics models applied is the Seasonal Naive model. The seasonal naive model makes the forecast using the last value from the same season, for example, the year before to forecast the value for the next year.
The second model applied is the exponentials state smoothing method by using ETS model who refer to error, trend and seasonality. This model can perform better in a short-term and on a univariate time series forecast. The model uses the exponentially weighted moving average (EWMA) to “smooth” a time series and trying to eliminate the random effect. The model uses a smoothing constant (a) which is transformed into a damping factor (1-a), the constant (a) define the weight which is applied to each period. ETS model can be applied as additive or multiplicative, but R selects the most optimal.
The third model applied is the autoregressive integrated moving average (ARIMA). ARIMA models are more complex model than the two previous models mainly because of the algorithm which backed this forecast model in R. ARIMA is the combination of two models. First, autoregressive model AR(p), which forecast the variable of interest using a linear combination of past values of the variable, where (p) is the lag number. Second, moving average models MA(q), which is applied as a linear regression of the current value of the series against current and previous white noise error terms or random shocks.
The fourth model applied, are neural network models (NN), which is the most complex model used in this analysis. This NN model performs in nonlinear time series and with big data sets. Because there seems nonlinearity in our data, we decided to test the predictive capacity of this model. The NN model is organised in multiples layers, the simplest networks contain no hidden layers and are equivalent to linear regressions. The coefficients attached to these predictors are called “weights”. The forecasts are obtained by a linear combination of the inputs. The weights are selected in the neural network framework by using a “learning algorithm”.
2.1 Analysis outline
First we upload data from the relevant working directory. How R and R Studio can be installed and data can be uploaded, watch R &R Studio and Data upload
Install R packages, load the data and declare this data series as a time series. Preliminary data observations
Data decomposition: Stationarity and identify lag specification, Seasonal component, Cycle component, Trend component.
Finding the most accurate model 4.1. Seasonal Naive method 4.2. ETS method, exponential smouthing models 4.3. ARIMA model 4.4. Neural network model Make the forecast for the next 12 months. Conclusion Create a report with Markdown
2.2 Package uploaded for this analysis
Some of these packages are not required for this analysis.
2.3 Upload data
This data is about monthly tax collection in billion of rupees (data source FBR).I have uploaded data from the respective directory.
Code
library(readxl)
tax_data <- read_excel("data/tax_billion_rs.xlsx") #Source FBR
tax_dataCode
Code
2.4 Time Plot
We have used “The Economist Themse”, one may use one’s own choice like Financil Times, Fivethirtyeight.com etc.
2.5 Multiple line plot
There are total taxes, direct and indirect taxes, and subcomponents of indirect taxes- sales tax, custom and excise duties. Now we plot direct and indirect taxes. How t
Code
ggplot(tax_data, aes(x=Date)) +
geom_line(aes(y = direct), color = "darkred") +
geom_line(aes(y = indirect), color="steelblue", linetype="twodash")+
scale_color_manual(values = c("darkblue", "darkred", "steelblue")) +ylab("Rs. in Billion") +
ggtitle("Monthly tax collection year wise in billion of Rs. by FBR")+ labs(caption=" Source :FBR")+theme_economist()#Method 1
2.5 Another approach for multiple line chart
Code
Code
# Visualization
ggplot(df, aes(x = Date, y = value)) +
geom_line(aes(color = variable, linetype = variable)) +
scale_color_manual(values = c("darkblue", "darkred", "steelblue")) +ylab("Rs. in Billion") +
ggtitle("Monthly tax collection year wise in billion of Rs. by FBR")+ labs(caption=" Source :FBR")+theme_economist()
We observe a strong positive trend and seasonality. Trend over time seems bit nonlinear. Even though, in this first observation, we can intuitively conclude that these time series have no stationary. In the next section, we will investigate and identify those potential problems in the way to make some correction to this time series and for selecting the most appropriate forecasting model.
3 Preliminary data decomposition
3.1 Why we investigate Stationarity? (unit root test)
A time series is stationary if it’s characteristics like mean, variance, covariance, are time variant: that is, they do not change over the time. Non-stationarity may cause autocorrelation which we explained in the next step.
We will make the Dickey-Fuller Test to check the stationarity in the data.
Augmented Dickey-Fuller Test
data: tax_data$direct
Dickey-Fuller = 4.5963, Lag order = 12, p-value = 0.99
alternative hypothesis: stationaryNull hypothesis of non-stationarity is not rejected. This test confirmed that this series is not stationary.
To correct the non-stationarity problem, we apply the first difference and make the Dickey-Fuller again.
Augmented Dickey-Fuller Test
data: DS
Dickey-Fuller = -3.6786, Lag order = 12, p-value = 0.02607
alternative hypothesis: stationaryThe Dickey-Fuller test allows you to reject the null hypothesis with this small p-value and we can conclude this series is stationary. By taking the first difference, we are making the making the corrections on this initial non-stationary time series.
We can visualize the data, for seeing the impact of the first difference on this time series.
Code
We will take the first difference from the data to remove the trend. With this first difference, we can work with this time series without having the trend influence forecasting .
3.2 What problems are caused by non-stationary (trend)
Also, It’s important because it helps to identify the driving factors. When we detect a change in a time series, we may be able to infer a correlation. But we need both time series to be stationary (no trend and no seasonality); otherwise, the correlation we find will be misleading.
There are three main problems with stochastic trends: 1. AR coefficients can be badly biased towards zero. This means that if you estimate an AR and make forecasts, if there is a unit root then your forecasts can be poor (AR coefficients biased towards zero)
2. Some t-statistics don’t have a standard normal distribution, even in large samples (more on this later)
3. If Y and X both have random walk trends then they can look related even if they are not – you can get “spurious regressions.”3.4 Why we investigate autocorrelation?
Autocorrelation means that there are correlations in the error or lag correlations of a given series with itself, lagged by a number of time units. Which signify:
\[Y_c =x + \beta X_i+u_i\]
\[Cov(u_i, u_s) \neq 0 \forall i\neq s\] Autocorrelation measures the linear relationship between lagged values of a time series. We will see dependence in the data across a range of lag value.
Code
acf(tax_data$direct)
acf plot shows that the data are strongly non-random and autorgressive model maybe appropriate. We can check the autocorrelation by plotting residual and standardized residuals of regression against time and compare if they show a similar pattern which signs for autocorrelation.
3.6 Seasonality
Seasonality is a pattern which occurs when a time series is observed at less than annual frequency and is affected by seasonal factors such as the time of the year or the day of the week. We know that tax collection authorities usually have to meet quarterly target and some accounting lags also make collection a seasonal phenomena. We need to detect seasonality in a time series in the way to make the necessary adjustment or for choosing the appropriate models. Seasonality adjustment has three main reasons: 1. to aid in short-term forecasting
2. to help in relating time series to other series or extreme events 3. to allow series to be compared from year to year or month to month or the day today.
Following two visualizations will help to identify the seasonality.
Code
library(ggthemes)
ggseasonplot(ts_tax,year.labels=TRUE, year.labels.left=TRUE) +
ylab("Rs. in Billion") +
ggtitle("Monthly tax collection year wise in billion of Rs. by FBR")+ labs(caption=" Source :FBR")+theme_economist()
We can observe that those 19 multiples colour lines, already have the same pattern over the years. Those colour lines informed you of the possible presence of seasonal cycles in this time series.
Now, look at this other seasonal plot which isolates the variation for one month at the time
Code
ggsubseriesplot(ts_tax)
The horizontal lines indicate the means for each month. This form of plot enables the underlying seasonal pattern to be seen clearly, and also shows the changes in seasonality over time. It is especially useful in identifying changes within particular seasons.
Both line and seasonal plots indicate that there seems strong trend and seasonality pattern in tax data.
4. Fit a Model
In this section we are going to fit a model on tax data and will find out the most appropriate model. ## 4.1 seasonal naïve method as our benchmark
This model is used as a benchmark and performance of any technical model is compared with this model. This model using the most recent observation as a forecast. We are using seasonal naive model because data seems seasonal.
\[\hat{y}_{T+h|T}=y_{T+h-m(k+1)}\]
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2022 -54.564050 -105.741349 -3.386751 -132.832975 23.70487
Sep 2022 -6.865504 -58.042803 44.311795 -85.134428 71.40342
Oct 2022 105.135952 53.958653 156.313251 26.867027 183.40488
Nov 2022 -40.176730 -91.354029 11.000568 -118.445655 38.09219
Dec 2022 6.499261 -44.678037 57.676560 -71.769663 84.76819
Jan 2023 95.839112 44.661813 147.016411 17.570188 174.10804
Feb 2023 -48.087819 -99.265118 3.089480 -126.356743 30.18111
Mar 2023 -3.683492 -54.860791 47.493806 -81.952417 74.58543
Apr 2023 65.958451 14.781152 117.135750 -12.310473 144.22738
May 2023 -41.394026 -92.571325 9.783273 -119.662950 36.87490
Jun 2023 12.427184 -38.750114 63.604483 -65.841740 90.69611
Jul 2023 125.949914 74.772615 177.127213 47.680989 204.21884
Aug 2023 -54.564050 -126.939680 17.811580 -165.253025 56.12492
Sep 2023 -6.865504 -79.241134 65.510127 -117.554478 103.82347
Oct 2023 105.135952 32.760321 177.511582 -5.553023 215.82493
Nov 2023 -40.176730 -112.552361 32.198900 -150.865705 70.51224
Dec 2023 6.499261 -65.876369 78.874892 -104.189713 117.18824
Jan 2024 95.839112 23.463482 168.214743 -14.849862 206.52809
Feb 2024 -48.087819 -120.463449 24.287812 -158.776793 62.60116
Mar 2024 -3.683492 -76.059123 68.692138 -114.372467 107.00548
Apr 2024 65.958451 -6.417179 138.334082 -44.730523 176.64743
May 2024 -41.394026 -113.769656 30.981605 -152.083000 69.29495
Jun 2024 12.427184 -59.948446 84.802815 -98.261790 123.11616
Jul 2024 125.949914 53.574283 198.325544 15.260939 236.63889Code
checkresiduals(fit_SN)
Ljung-Box test
data: Residuals from Seasonal naive method
Q* = 91.108, df = 24, p-value = 9.441e-10
Model df: 0. Total lags used: 24We have Residual sd: 42.8065 from summary(fit_SN) command. The acf graph of residuals indicate autocorrelation among the residuals. Ljung Box statistics is also very high and it reject the null hypothesis that residuals are uncorrelated.
The naive and Snaive model are fundamental models. Some business uses those model basic forecasting models, maybe because of the lack of internal resources. Producing or maintaining extra stock it is a cost for the company and creates inefficiency. We continue testing the forecast performance of others models.
4.2 Fit ETS method, exponential smoothing method
Second, we will apply ETS model: Error, Trend, Seasonal. The flexibility of the ETS model lies in its ability to trend and seasonal components of different traits. This function ets() automatically optimizes the choice of model and necessary parameters. We present the structure of the additive and the multiplicative form.
Assuming: \[\mu_t = \hat{y}_t=l_{t-1}+ b_{t-1}\] and \[\varepsilon_t = y_{t} - \mu_{t}\]
ETS additive
\[y_t = l_{t-1} + \phi b_{t-1} + \varepsilon_t\]
\[l_t = l_{t-1} + \phi b_{t-1} + \alpha\varepsilon_t\]
\[b_t = \phi b_{t-1} + \beta^*(l_{t}-l_{t-1}- \phi b_{t-1} = \phi b_{t-1}+\alpha\beta^*\varepsilon_t\]
\[\varepsilon_t = (y_t-\mu_t ) / \mu_t\]
\[y_t = (l_{t-1} + \phi b_{t-1}) (1+\varepsilon_t)\]
\[l_t = (l_{t-1} + \phi b_{t-1}) (1+\alpha\varepsilon_t)\]
\[b_t = \phi b_{t-1}+\beta(l_{t-1} + \phi b_{t-1})\varepsilon_t\]
Code
fit_ets <- ets(ts_tax) #residual = 221.1631
checkresiduals(fit_ets)
Ljung-Box test
data: Residuals from ETS(A,N,A)
Q* = 115.63, df = 24, p-value = 5.784e-14
Model df: 0. Total lags used: 24We have a residual sd =40.49 from summary(fit_ets) command, which are more accurate than the seasonal naïve models and what mean for the exact month the years before, which missing on average 40.49 billion.
So this model increases the precision and offer a better fit but: if we look to AFC graph, we observe that there remains autocorrelation because of the bar going out of the 95% confidence line.
The Ljung-Box test Test results also indicate rejection of null hypothesis that residuals are non-correlated.
We realize that with using just a bit more complex forecasting model, we increase accuracy and can make a significant improvement in tax collection. We continue our analysis with a one of the most performant forecasting model the ARIMA.
4.3 Fit on ARIMA model
ARIMA model is a Generalized random walk model which is fine-tuned to eliminate all residual autocorrelation. It is a Generalized exponential smoothing model that can incorporate long-term trends and seasonality.
AR(p) model \[(1-\sum^p_{k=1}\alpha_kL^k)X_t = \varepsilon_t\] MA(q) model \(X_t = (1+\sum^q_{k=1}\beta_kL^k)\varepsilon_t\) First difference is \(\Delta X_t=X_t -X_{t-1} = (1-L)X_t\) where \(\ Y_t = (1-L)X_t\)
ARIMA(p, d, q) full model \[(1-\sum^p_{k=1}\alpha_kL^k)(1-L)^dX_t = (1+\sum^q_{k=1}\beta_kL^k)\varepsilon_t\]
ARIMAX model
\[\Delta y_t=\alpha_0+\sum_{j}\alpha_j \Delta y_{t-j}+\sum_h\gamma_h\epsilon_{t-h}+X\beta+\epsilon_t\]
Code
#fit_ARIMA <- auto.arima(ts_tax, d=1,D=1, stepwise = FALSE, approximation = FALSE)
#checkresiduals(fit_ARIMA)We have a residual standard error (residual) = sqrt(1251), what mean for the exact month the years before, which missing on average 35 billion.
At this stage, we can conclude that the ARIMA model offers the best fit base on residual and standard deviations.
The Ljung-Box test Test results also indicate that model is appropriate fit as null hypothesis of non-autocorrelation among residuals is rejected. So far ARIMA model not only passes the diagnostic testing but also performs better than the seasonal naive and ETS models.
Though for economic theory or behavioral models, simpler the better. But for forecasting having capacity for applying and Understanding complex models is a plus value for any business or forecasting team. Many decisions are based on the capacity to forecast revenues. Accuracy in short-term revenue forecast has the critical impact on all economic planning. Accuracy on long-term forecasting, contribute for supporting the significant fiscal decision. In the end, short-term and long-term forecast accuracy will have a significant impact to optimize the growth of the economy. Lousy forecasting can reduce FBR capacity for effecient decision making.
4.4 Neural network models
This model allows complex nonlinear relationships between the response variable and its predictors. Our time series are mostly linear, and in this case, the neural network may not be the most appropriate model, but we will compare performance against seasonal naive, ETS and ARIMA.
In neural network the inputs of each node are combined using a weighted linear combination. For exemple, the inputs (for 4 inputs) into hidden neuron j are combined linearly to give \[z_j = b_j+\sum_{i=1}^4 w_i,_jx_i\] and in the hidden layer,this is then modified using a nonlinear function such as a sigmoid, \[s(z) = \frac{1}{1=e^-z}\]
Code
NNL <- nnetar(ts_tax)
NN <- nnetar(ts_tax, lambda="auto")
checkresiduals(NNL)
accuracy(NN)There seems some spikes going out of usual range in acf plot of residuals which may indicate little problem of residuals autocorrelation. This is very complex model but often very useful for deep learning.
4.5 Comparative forecast performance
We will compare those forecast models performances with using two strategies.
First strategy: We will compare the indicators produced by the output of the accuracy() function when all the time series is the training set.
Second strategy: Like the first strategy, we will compare the indexes produce by the output of the accuracy() function but we will split the dataset by creating, one training set and one test set.
It is important to use the right tool for an accurate analysis and for understanding the impact on the results when a different strategy is used.
We will base ours performance comparison model analysis on those indexes which will we define it before we use it.
Mean absolute error: MAE
\[MAE=mean(e^2_{t})\] Root mean squared error: RMSE
\[e_{t}={y}_{t}-\hat{y}_{t|N}\] \[RMSE=\sqrt{mean(e^2_{t})}\] Mean absolute percentage error: MAPE \[p_{t}=100e_{t}/y_{t}\] Mean absolute scaled error : MASE \[q_{t}= e_{t}/Q\] Where Q is a scaling statistic computed on the training data. \[p_{t}=100e_{t}/y_{t}\] Autocorrelations of error at lag1 (ACF) \[ACF=\frac{Covariance(x_{t},x_{t-h})}{Variance(x_{t})}\] If we print summary of each model and compare accuracy, we observe that ARIMA model seems the best while NN model seems the second best. Therefore, now we forecast with NN and ARIMA models for 24 time periods ahead.
4.4 Forecast Using ANN model
Code
knitr::opts_chunk$set(echo = TRUE)
NN <- nnetar(ts_tax, lambda="auto")
fcast2 <- forecast(NN, h=24, PI=TRUE, npaths=100)
autoplot(fcast2, include = 60)
Forecast method: NNAR(7,1,4)[12]
Model Information:
Average of 20 networks, each of which is
a 8-4-1 network with 41 weights
options were - linear output units
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 3.265063 29.61893 16.87076 -Inf Inf 0.7405954 -0.1497553
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2022 -5.619309e+01 -115.9590461 -1.677937e+01 -154.03614 -6.2763381
Sep 2022 -3.926319e+00 -47.5600567 -4.083768e-04 -92.85391 0.1923545
Oct 2022 1.300172e+02 57.7724784 2.673309e+02 38.04682 313.6271246
Nov 2022 -5.212019e+01 -141.6192785 -1.212580e+01 -267.93625 -5.2007405
Dec 2022 1.711498e+00 -1.7817644 5.410234e+01 -12.20932 137.0238907
Jan 2023 1.226421e+02 41.1059691 2.314269e+02 29.51041 354.8527644
Feb 2023 -5.386332e+01 -121.4497681 -5.730489e+00 -213.56927 -0.2054690
Mar 2023 -4.573142e-01 -25.1830719 8.063896e+00 -47.92359 28.0007977
Apr 2023 1.011079e+02 36.2166958 2.426411e+02 19.16445 343.6042723
May 2023 -5.145643e+01 -120.0443599 -6.458789e+00 -190.56144 -0.4129747
Jun 2023 3.594046e+00 -0.7570878 7.073977e+01 -12.56524 124.3842624
Jul 2023 1.352299e+02 50.7347396 2.322142e+02 12.33302 324.7208425
Aug 2023 -5.788544e+01 -151.1809694 -1.064100e+01 -221.39087 -4.6071002
Sep 2023 -1.264988e+00 -16.6411900 5.051960e+00 -31.95307 43.0203901
Oct 2023 1.475782e+02 45.8018802 2.627653e+02 26.40376 357.3545812
Nov 2023 -5.431170e+01 -125.2161253 -5.986957e+00 -202.52984 -1.5782961
Dec 2023 3.303619e-01 -7.4224151 3.295390e+01 -41.77100 95.7063215
Jan 2024 1.373830e+02 34.1849492 2.302237e+02 12.78367 293.1243548
Feb 2024 -5.767457e+01 -123.8418655 -5.318603e+00 -211.33518 -0.5965792
Mar 2024 2.205332e-04 -12.3850050 8.360754e+00 -36.26833 63.1119778
Apr 2024 1.396638e+02 41.4789942 2.352737e+02 25.16609 345.7797895
May 2024 -5.312843e+01 -140.7811118 -1.080373e+01 -223.27575 -1.7500398
Jun 2024 4.156411e-01 -8.3394684 5.937000e+01 -21.36912 94.8124719
Jul 2024 1.482162e+02 34.7721363 3.030343e+02 18.32510 435.9820062Forecast Using ARIMA model
Code
#fit_ARIMA %>% forecast(h=24) %>% autoplot()
#auto.arima(tax_data$direct) %>% forecast(h=24) %>% autoplotCode
ggplot(tax_year1)+aes(x=year)+geom_line(aes(y = dir_yr), color = "darkred") +
geom_line(aes(y = indir_yr), color="steelblue", linetype="twodash")+
scale_color_manual(values = c("darkblue", "darkred", "steelblue")) +ylab("Rs. in Billion") +
ggtitle("Indirect (dashed) and direct taxes (solid line) in billion of Rs. by FBR")+ labs(caption=" Source :FBR")+theme_economist()
Code
