Application of Vector Autoregressive (VAR) Process in Modelling Reshaped Seasonal Univariate Time Series
Chepngetich Mercy^{1}, John Kihoro^{2}
^{1}Jomo Kenyatta University of Agriculture and Technology, School of Mathematical Sciences, Nairobi, Kenya
^{2}Cooperative University College of Kenya, Department of Computing and elearning, Nairobi, Kenya
Email address:
Chepngetich Mercy, John Kihoro. Application of Vector Autoregressive (VAR) Process in Modelling Reshaped Seasonal Univariate Time Series. Science Journal of Applied Mathematics and Statistics. Vol. 3, No. 3, 2015, pp. 124135. doi: 10.11648/j.sjams.20150303.15
Abstract: Seasonal Autoregressive Integrated Moving Averages (SARIMA) model has been applied in most research work to forecast seasonal univariate data. Less has been done on Vector Autoregressive (VAR) process. In this research project, seasonal univariate time series data has been used to estimate a VAR model for a reshaped seasonal univariate time series for forecasting. This was done by modeling a reshaped seasonal univariate time series data using VAR. The quarterly data is reshaped to vector form and analyzed to vector form and analyzed using VAR for the year 1959 and 1997 to fit the model and the prediction for the year 1998 is used to evaluate the prediction performance. The performance measures used include; mean square error (MSE), root mean square error (RMSE), mean percentage error (MPE), mean absolute percentage error (MAPE) and Theil’s U statistic. Forecasting future values from the fitted models in both SARIMA and VAR using Box Jenkins procedures, (Box and Jenkins; 1976) was done. The results showed that both models are appropriate in forecasting but VAR is more appropriate model than SARIMA model since its predictive performance was shown to be the best. Other data sets were also used for analysis and comparison purpose.
1. Introduction
1.1. Background of the Study
The goal of a time series analysis is to obtain the parameters of the underlying physical process that governs the observed time series and used it to forecast future values. The modelling and predictions is meant to determine the best model for forecasting. This is useful to many areas of study such as meteorological forecasting, electricity consumption, demographic forecasting and Economic growth. In these forecasts, both short term and long term seasonal trends can be forecasted. To model these variations for forecasting, usually a seasonal ARIMA model of Box and Jenkins (1976) is used. ARIMA time series is a popular models used for studying weather and market data
Seasonal ARIMA is a model of ARIMA class suitable for forecasting seasonal data variations and trends. SARIMA model exhibit some disadvantages that; constructing it is expensive. Also, adapting it for use requires expertise and consumes a lot of time. Most of the univariate data are seasonal. Hence, fitting the model that gives the best performance prediction helps a lot in future predictions and analysis. This will be of great importance to economic and planning purposes. For seasonal univariate time series data, SARIMA model is usually used to do forecasting, with its advantages and limitations. In this project SARIMA model is reshaped and forecast using Vector Autoregressive (VAR) model.
This Project applied VAR in the reshaped seasonal univariate time series data, test and distinguish between the two models in terms of its predictive performance results. The results show that although both models are appropriate in forecasting, VAR model gives the best performance compared to the SARIMA model.
1.2. Literature Review
Sims (1980) introduced VAR as a technique that could be used by microeconomics to characterize the mutual dynamic behavior of collection of variables without requiring strong restrictions of the kind needed to identify underlying structural parameters. Sims developed a VAR model with p lags, VAR (p) for expressing a set of variables as a weighted linear combination of each variable's past values and the past values of the other variables in the set. However, multivariate data analysis in the context of VAR has evolved as standard instrument in econometrics.
The VAR (p) models are more flexible than AR models. Mei, Liu and Jing (2011) constructed a multifactor dynamic system VAR forecast model of GDP by using six important economic indicators. This included the fiscal revenue, social retail goods, secondary industry output, and investment in fixed assets, employment rate and tertiary industry output which based on data from the shanghai region in china. In their analysis, the model showed high significance and the results show that the relative error forecast is quite small. The authors therefore conclude that the VAR model has a considerable practical value.
Clarida and Friedman (1984) use a VAR model to forecast the United States shortterm interest rates during April 1979 to February 1983. A constantcoefficient, linear VAR model is generated to estimate the preOctober 1979 probability structure of the quarterly data , which takes six important united states macroeconomic factors into consideration. The result shows that shortterm interest rates in the United States have been too high since October 1979. Based on their VAR model, the prediction results of conditional and unconditional forecast are both lower than the actual United States shortterm interest rates during this period.
The VAR model extends univariate autoregressive models to dynamic multivariate and provides better forecasts than univariate time series models (Zivot and Wang, 2006). VAR models are used to describe and forecast multivariate time series for stationary time series. For nonstationary time series a vector error correction term is added to form a vector error correction model (VECM) and it is necessary to test for the existence of a stationary linear combination of the nonstationary terms (co integration). It must be transformed into vector error correction model (VECM) by taking the first difference.
Autoregressive integrated movingaverage (ARIMA) models are simple time series models that can be used to fit and forecast univariate data such as fisheries landings. With ARIMA models, data are assumed to be the output of a stochastic process, generated by unknown causes, from which future values can be predicted as a linear combination of past observations and estimates of current and past random shocks to the system (box et al., 2008).
Xiao Han Cai, (2008), uses both VAR model and SARIMA model to analyze time series data of air pollution co in California south coast area. Their results showed that VAR model is a better model to forecast multiple variables data sets though not easy to find the order to fit an accurate VAR model. On the other hand, SARIMA model presents how the current month air pollutant concentrations depend upon the previous month’s air pollutant concentrations.
1.3. Statement of the Problem
Modeling a seasonal univariate data, usually a SARIMA model is applied. This is because of the seasonal variations exhibited. This seasonal variation could either be reshaped into vector form, making it a multivariate time series then forecasted using VAR
VAR has been seen as a more advanced most successful, flexible and easy to use model for forecasting (Zivot and Wang, 2006). This calls for a research on this to examine if VAR can be a better tool to use in the reshaped time series than the SARIMA model for forecasting a seasonal univariate time series.
2. Methodology
In this study we will first examine the appropriateness of VAR model in forecasting a reshaped seasonal univariate time series by predicting one step ahead vector observation. Then compare the forecast with those obtained with univariate SARIMA model.
2.1. Vector Autoregressive Process
In order to build a VAR model, (box and Jenkins; 1976) steps can be followed. This includes model identification, estimation of constants, diagnostic check and finally forecasting. Conditional heteroscedasticity and outliers of the residual series is also checked. The existence of co integration is used to check the presence of any common trends. Error Correction Model (ECM) can be developed in case of co integration. This improves the forecasting of long term.
2.1.1. Model Specification and Identification
Just like any other univariate time series, the first step in building a VAR (p) model involves model identification. This helps in identifying the order of the appropriate model. The most common information criterion used to identify the model include Akaike information criterion (AIC) (box and Jenkins; 1976), SchwarzBaysian (BIC) and Hannanquinn (HQC).
The time series y_{t , }where y_{t}=(y_{1t },y_{2t},…,y_{nt} ) denote an (n×1) vector variables, follows a VAR(p) model if it satisfies
Y_{t} = п+ Φ_{1}y_{t1}+…+Φ_{tp}y_{tp }+ e_{t , }t=1,2,3,…t, (1)
In matrix form, this can be expressed as;
Where п is a kdimensional vector, Φ_{1,} Φ_{2,}… ,Φ_{p} are k×k parameter matrices and e_{t }is a sequence of independently and identically distributed error vectors.
Assumptions of the errors
The error term e_{t }is a multivariate normal k × 1 vector of error satisfying the following assumptions;
1. E(e_{t})=0 every error term has mean zero;
2. E(e_{t}eˈ_{t})=s the contemporaneous covariance matrix of error terms is ω (a n× n positivesemidefinite matrix);
3. E(e_{t}eˈ_{tk})=0 for any nonzero k. There is no correlation across time; in particular, no serial correlation in individual error terms.
and Φ_{j }are k×k matrices.
Using the backshift operator b, the VAR (p) model can be written as
(I – Φ_{1} b– ….–Φ_{p}b^{p}) Y_{t} = п+ e_{t,}
Where I is the k×k identity matrix. In a compact form as follows
Φ(b) Y_{t}= п+ e_{t,}
Where Φ (b) =I Φ_{1} b…Φ_{p}b^{p}is a matrix polynomial.
Stability in VAR (p) process is one important characteristic. It generates stationary time series, where the means, variances and covariance are time invariant. Reverse characteristic polynomial can be evaluated to check its stability.
det(I – Φ_{1} b– ….– Φ_{p }b^{p})^{}≠0, b≤1.
If b =1,in the solution, then some variables are integrated
If Y_{t }is weakly stationary, then we have
µ =E(Y_{t}) = (I – Φ_{1} b– ….– Φ_{p }b^{p})^{1}п = [ Φ(1)]
Provided that the determinant exists since determinant of [Φ(1)] is different from zero. Expressing VAR (p) in deviation form from its mean, we define µ =E(Y_{t}). Then in deviation form is given by;
Ỹ_{t}= Y_{t} µ, then the VAR(p) model becomes
Ỹ_{t} = Φ_{1} Ỹ_{t1} +…+ Φ_{p} Ỹ_{tp}_{ }+ e_{t} (2)
Given, as the residual of the model, the residual covariance matrix is defined as;
,
The AIC of a VAR (i) model under the normality assumption is defined
We will select the VAR of order p such that the calculated value is minimum, that is , where p is an integer for a given vector time series. Other multivariate information criterion measures include;
Where is the total number of parameters in each equation, and the VAR model of order p lags such that the criterion information is minimum, is selected.
2.1.2. Estimation of Constants
After obtaining the order of the model, p, of the vector series, we now derives the estimators of the constants.
Consider the consecutive VAR models:
Y_{t}= п+ Φ_{1} y_{t1} +e_{t}
Y_{t =} п + Φ_{1} y_{t1} + Φ_{2} y_{t2 }+e_{t}
…=…
Y_{t} = п + Φ_{1} y_{t1} +…+ Φ_{i}y_{ti }+e_{t} (3)
The most common methods of estimating parameters are the maximum likelihood estimator (MLE) and the ordinary least square estimator (yang & yuan 1991). Here, we will apply the ordinary least squares (OLS) method to estimate the parameters of these models and applied equation by equation.
For equation (3), let _{i}^{(j) }be the OLS estimate of Φ _{i }and _{i}^{(j)} be the OLS estimate of п_{. }Where (i) is used to denote the estimate of a VAR (i) model. The estimates of the coefficients of the VAR model, Φ_{1} Ỹ_{t1} +…+ Φ_{p} Ỹ_{tp}_{ }+ e_{t}, when estimated using unrestricted VAR (p) model are considered to be fixed quantities. These estimates of coefficients do not accurately reflect the underlying relationship because some of the estimated coefficients of the VAR model are nonzero purely by chance when estimated by OLS so restrictions may be imposed to reduce the number of parameters being estimated.
2.1.3. Diagnostic Check
Suppose the orders and constants have been chosen for a VAR model underlying the data, then residuals are checked whether they are normally, identically and independently distributed.
(i) Statistical Tests
The portmanteautest of box and pierce (1970)
It checks whether the estimated residuals behave approximately like realizations from a white noise process. It is defined by;
Where, . The test statistic has an approximate of chisquare with degrees of freedom where is the number of coefficients excluding deterministic terms of a VAR (p) model.
(ii) Arch Tests
The multivariate arch tests, is used to test for heteroskedasticity. The test statistic is defined as
Where and assigns the covariance matrix of the model. The test statistic is a chisquare distribution with degrees of freedom and is implemented in r using vars package.
(iii) Normality Check
The JarqueBera normality tests for multivariate series are implemented and applied to the residuals of a VAR (p) as separate tests for multivariate kurtosis and Skewness. The test statistics is defined as
Where and are computed as follows;
Where and are the third and fourth noncentral moment vectors of the standardized residuals. The is a chi square distribution with 2k degrees of freedom. The multivariate Skewness and kurtosis test are chi square distribution with 2k degrees of freedom.
(iv) Unit Root Tests
(a) DickeyFuller Test & Augument DickeyFuller Test
DF tests the null hypothesis of the unit root against the alternative that there is no unit root in the process. It gives a formal test for checking stationarity of the model. For example, given a differenced equation of an AR (1) as,
Then if the is a random walk, then, but if is stationary, Then coefficient is negative. Using the standard ttest statistic, it will be formed as
and are the least squares estimators for α and ^{2}
For large n, the statistic converges to functional distribution of wiener process
Where w is a standard wiener process.
ADF statistic is an augmented version of the df test. It is used to test for a more complicated and larger set of time series models. The test used in ADF statistic is a negative number and the more negative, the stronger the rejection of the null hypothesis. It can be applied to each variable to check for existence of co integration.
(b) Durbin Watson
This is a test statistic used to detect the presence of autocorrelation in the residuals from a regression analysis. It tests for the null hypothesis that the residuals are serially uncorrelated against the alternative hypothesis that residuals follow a stationary first order auto regression.
Given that is the residual associated with the observation at time t, and then the Durbin Watson test statistic is given by;
Where t is the number of observations and . R is the sample autocorrelation of the residuals.
If d=2, there is no autocorrelation. Then if d<2, there is evidence of positive correlation. For d>2,
Then error terms are negatively correlated which imply an underestimation of level of statistical significance. But if d<1, causes an alarm.
(c) The Qk (m) Statistic
Qk(m) is used to check the adequacy of the fitted model. It can be applied to the residual series to check the assumption that there are no serial or crosscorrelations in the residuals. If the model is fit, the Qk (m) statistic of the residuals is asymptotically a chisquared distribution with mg degrees of freedom, where g is the number of estimated parameters in the AR coefficient matrices.
(d) The T Statistic
tstatistic is used to test the statistical significance of parameters. It is carried out to check if the model is over specified. It is also important to assess whether the stationarity and invertibility conditions are satisfied. If we factorize characteristic polynomials and one of their roots is close to unity, it may be an indication of lack of stationarity and or invertibility.
An inspection of the covariance matrix of the estimated parameters allows us to detect the possible presence of high correlation between the estimates of some parameters which can be a manifestation of the presence of a common factor in the model (Box and Jenkins; 1976).
A vector error correction model (VECM) can be used to incase of presence of a common factor (co integration) in the model. It describes the nature of any non stationarity among different component series. Applying it improves longer term forecasting over an unconstrained model.
2.1.4. Prediction
Prediction is the estimation of values whether future, current or past with respect to the given data. Future values of the original process can be predicted based on the model which fits the given data. It starts with certain assumptions and the estimates projected into the coming future either weeks, months and even years using techniques such as boxJenkins models, exponential smoothing, regression analysis, moving averages and projection. Any error in the assumptions results in a similar or magnified error in forecasting. Therefore, the sensitivity techniques for analysis are used to assign a range of values to uncertain factors or variables.
Once the model has been identified and passed through diagnostic check, we can use it for it for forecasting. For time , the basic problem is then to estimate future values
Of the hsteps ahead forecast made at time t. We need to forecast
In such a way that the mean squared error (MSE) of the prediction is minimum.
For a VAR (p) model, the 1step ahead forecast at the time origin h
Is
The associated forecast error is a_{h }=e_{h}+1
The covariance matrix of the forecast error is Σ. If Y_{t }is weakly stationary, then lstep ahead forecast Y_{1}(l)converges to its mean vector, µ, as the forecast horizon increases.
(i) Advantages of VAR Modeling
1. VARs are more flexible than univariate AR models by allowing the value of a variable to depend on more than just its own lags or combinations of white noise terms
2. The researcher does not need to specify which variables are endogenous or exogenous with VAR model. All variables are endogenous.
3. Forecasting with VARs are often better compared to ‘traditional structural’ models.
(ii) Disadvantages of VARs
1. VARs use little theoretical information about the relationships between the variables to guide the specification of the model.
2. Often not clear how the VAR estimates of coefficient should be interpreted.
3. There are so many parameters to be estimated.
2.2. SARIMA Model
Seasonal autoregressive integrated moving average (SARIMA) processes are designed to model time series with trends, seasonal patterns and short time correlation. They have developed from the standard model of box and Jenkins (1970) and incorporate both seasonal autoregressive and moving average factors into the modeling process.
It is a class of ARIMA models suitable for data exhibiting seasonal variations. Suppose that a time series has a polynomial trend of degree d. Then we can eliminate this trend by considering the process, obtained by d times differencing.
If the filtered process , is an ARMA (p, q) process satisfying the stationarity condition, the original process () is said to be an autoregressive integrated moving average of order p,d,q, denoted by ARIMA(p,d,q). In this case constants exist such that is a white noise.
Let , t=0, 1, 2 … be a nonstationary time series possibly containing seasonality.
Then depends on past values such as as well as Where is a differencing operator, d is the order of nonseasonal differencing is ARMA (p, q) Process.
Then is ARIMA (p, d, q) model but if
Where d is the order of seasonal differencing, the model is referred to as seasonal autoregressive integrated moving average SARIMA written as ARIMA (p,d,q)(P,D,Q)
Is given by
Where and , are polynomials of order p, P, q, Q respectively, z is a backshift operator, s is the seasonal period of the series. is an AR process of order p, is an AR process with order of seasonal component, is an ma process of order q and is an MA process with order of seasonal component.
Fitting SARIMA models to the meagre data using a semiautomated approach based on a combination of the boxJenkins method with smallsample, biascorrected Akaike information criteria (AIC) model selection (Rothschild et al., 1996; Brockwell and Davis, 2002). This approach involved three major steps:
2.2.1. Model Identification
Choosing the parameters p and q, the order q of a moving average MA (q) process can be estimated by means of the empirical autocorrelation function r (k)
The order p of an AR (p) process can be estimated in an analogous way using the partial autocorrelation function. MA (p, q) process we take the pair (p, q), which minimizes some function i. e Akaike's information criterion, Hannan information criterion with bias correction and the Bayesian information criterion discussed earlier under VAR.
2.2.2. Estimation of the Model Coefficients
The coefficients and are estimated using maximum likelihood method or otherwise. Most of the estimation methods are already implemented in existing computer software using iterative techniques.
2.2.3. Diagnostic Check
The fit of the ARMA (p, q) model with the estimated coefficients is checked. These involve scrutinizing the estimated residuals and ensure that they behave approximately like the realizations from a white noise process.
2.2.4. Forecasting
The forecast of future values of the original process is based on the model which is adequate and fits the given data.
2.3. The Concept of ReShaped Time Series
Given n observations in univariate time series data and assuming the series is seasonal with period s.
For illustration, we will use (s=4), where each observation has 4 variables i.e and which in our case will be seasons (s=4 variables). Arranging the data in a serpentine manner for all the n –observations , we will treat each row as a vector of , where i=1,2,…,t. The vector time series is of order s×t.
For a seasonal univariate time series data, we can reshape it as shown in table 1 below;



 



 






























real value 
For , then the VAR form of the reshaped time series is written as;
Y_{t} = п+ Q_{1}y_{t1}+…+Q_{tp}y_{tp }+ e_{t, }t=1,2,3,…t
In matrix form it is represented as follows;
2.4. Performance Measures
Sstep ahead forecast, , will be equivalent to onestep ahead forecast, , for vector.in order to get the best forecasting, the parameters of the models are adjusted to minimize the forecasting error. By defining a forecasting error (loss function), we search for the best parameters in both models that minimizes this function.
To compare performance of the VAR and SARIMA models, we compute a set of indicators for the quality of time series forecast methods. These include;
2.4.1. Mean Squared Error (MSE)
For the forecast in both predictions should be minimum
The results of both forecasts will be compared with the real values.
2.4.2. Root Mean Squared Error (RMSE)
Is the square root of the average of all squared errors, according to Wang and Lim (2005). It ignores any over and under estimation
2.4.3. Mean Percentage Error (MPE)
MPE is a percentage error measurement which allows comparison of under and overestimation. It takes into account whether a forecasting method is biased.
×100
2.4.4. Mean Absolute Percent Error (MAPE)
It is a percentage error measurement, which allows comparison of under and overestimation. It is particularly useful when the units of measurement are relatively large.
2.4.5. Theil’s U Statistic
Theil's Ustatistics see Theil (1958) is used as a measure of forecasting error that is minimized. It is a relative measurement based on comparison of the predicted change with the observed change. The value of u lies between 0 and 1. If u equals to 0, there is a perfect fit, whereas u equals to 1 implies that forecasting of data is very poor.
3. Results
3.1. Introduction
This chapter presents results using secondary data in which the concept of SARIMA model has been known. The data set was used to test the concept of VAR and examine its appropriateness. This was analyzed by use of statistical software for model fitting and testing. Akaike information criterion (AIC) is used to select the best fitting model. Bayesian information criterion (BIC) also was used.
3.2. Seasonal ARIMA Model
3.2.1. Testing the Stationarity
Macro dataset used by stock and Watson in their introduction to econometrics where all unemployment rates for 16 years and above data is used for analysis. We denote the quarterly data on macro dataset as unemploymentrates. Time series plot shows evidence of cyclic and random components. The spikes display evidence of strong seasonality. The time series plot is shown in figure 1 below.
Autocorrelation function and partial correlation function plot of unemploymentrates shown below is evident of presence of a unit root. This is shown by the PACF value close to 1 confirming nonstationarity. The ACF on the other hand shows the autocorrelation between the variables and the lags of itself. It is evident that in all the lags, there is statistically significance of autocorrelation at 1%. ACF and PACF plot is shown in figure 2.
We also test for presence of unit root test using ADF test. The ADF test results show that pvalue is greater than 0.01 hence we fail to reject the null hypothesis of non stationarity of the process and conclude that the data is nonstationary. We take the logarithm transformation of the data in order to reduce problems of heteroscedasticity. Taking the first difference of the data to ensure stationarity assumptions of ARIMA model is satisfied, the data shows stationarity with its plot. Testing for the presence of unit roots using ADF test, the pvalue is less than 0.01. We therefore reject the null hypothesis of nonstationarity and conclude that the data is stationary at critical value of 1%. This is shown in table 2.
Variable  ADF test statistic  Pvalue  Critical value 
Unemploymentrates  1.58286  0.4914  0.01 
Δlnunemploymentrates  4.31345  0.0004159  0.01 
We can also test the stationarity using KPSS test with the null hypothesis of stationarity. The tstatistic is less than the critical value at 1%. We therefore fail to reject the null hypothesis of stationarity at 1%. Therefore, we conclude that the data is stationary. This is shown in table 3 below.
Test statistic  Pvalue  
KPSS test statistic  0.0878173  0.738 
Both plots of differenced data show that there is a little autocorrelation in ACF. The PACF shows that the there is no spike close to one. This confirms the stationarity of the data. The plots are shown in figure 4.
3.2.2. Selection of the Order
The order of the model is checked to identify the appropriate SARIMA model. The model which minimizes the information criterion is chosen. In our data, the ARIMA
(0, 1, 1×0, 1, 0)_{4} model is found to be the best model since its information criterions are minimum.
3.2.3. Estimation of coefficient Parameters
From the fitted model, we can obtain its coefficients. The coefficients of ARIMA (0, 1, 1×0, 1, 0)_{4 }model based on the information criterion is shown below with its coefficients being statistically significant the AIC value is minimum and standard error variance is also small.
Model 30: ARIMA, using observations 1960:2–1999:4 (t = 159)
Dependent variable: (1 − L)(1 − L^{s}) lnunemploymentrates standard errors based on hessian coefficient std. Error z pvalue
Const 0.000773 0.00927 0.0833 0.9336
Θ_{ } 0.928806 0.0809547 1.4732 0.0000
Mean dependent var 0.000312
s.d. Dependent var 0.074352
Mean of innovations 0.000050
s.d. Of innovations 0.060778
Loglikelihood 218.6787
Akaike criterion −431.3574
Schwarz criterion −422.1507
Hannan–quinn −427.6187
Real Imaginary Modulus Frequency
MA
root 1 −1.0767 0.0000 1.0767 0.5000
The results are shown in table 4 below.
Coefficient  Std. Error  Z  Pvalue  
constant  1.000  0.0226214  44.21  0.000 
Thetha1  0.576546  0.0840321  6.861  6.84e012 
S.d dependent var.  0.038212  Aic  629.6445  
Loglikelihood  318.8223  Hannanquinn  624.9712 
From the fitted model we can see clearly that the pvalues of the coefficients are statistically significant. The equation for the fitted model is expressed as;
(1−Z)(1−Z^{4})Xt = (1 + 0.928806B)ε_{t}
3.2.4. Diagnostic Check
To check if the fitted model is significant, one of the ways of checking this is by plotting the frequency distribution of the residuals. Figure 5 below, shows the frequency distribution of the residuals of the fitted model ARIMA (0, 1, 1×0, 1, 0)_{4}. This shows that the residuals are normally distributed since the pvalue is greater than 0.01.
We can also check the normality of the residuals using the QQplot. From the plot, there are no outliers and all points lie on the line making us to conclude that the residuals of the fitted model are normally distributed.
3.2.5. Forecasting
Forecasting for ARIMA (0, 1, 1×0, 1, 0)_{4} model is shown in the table below.
Variable  Actual value  Prediction value  Error 
1998:1  4.6667  4.3456  0.3211 
1998:2  4.4333  4.7373  0.304 
1998:3  4.5000  4.0605  0.4395 
1998:4  4.43333  4.7509  0.3176 
3.3. VAR Modeling
The seasonal univariate time series data is reshaped to form vector form of four variables since the data is quarterly. Data between 1959 and 1998 are used insample estimation and data between 1998 and 1999 are used for the outofsample forecasting purposes. Figures below shows the time series plots of the four variables during the sample period. The variables are denoted Q1, Q2, Q3, and Q4. The figure below displays all the plots.
3.3.1. Testing Stationarity
Autocorrelation function and partial correlation function plot of the variables shown below in figure 7, 8, 9 and 10 is evident of presence of a unit root. This is shown by the PACF value close to 1 confirming nonstationarity. The ACF on the other hand shows the autocorrelation between the variables and the lags of itself.
We can also check the stationarity by testing the unit roots. One of the tests used is the augmented dickeyfuller test applied to the series in order to test for unitroots. Table 6 below shows the ADF results of both the variables and their first differences of the series.
From the table above, the ADF test results indicate that all variables are nonstationary by not rejecting the null hypothesis of unitroot for all the variables at 1% critical value, but they are all stationary after first differencing. We can therefore conclude that all the time series are integrated of at least order one; therefore, we use differenced series in our analysis.
Variable  Deterministic terms  Test value  Pvalue  Conclusion 
Q1  Constant, trend  2.2501  0.4496  Nonstationary 
Differenced q1  Constant  5.57543  3.236e005  Stationary 
Q2  Constant, trend  2.34281  0.402  Nonstationary 
Differenced q2  Constant  5.98549  1.264e005  Stationary 
Q3  Constant, trend  2.57543  0.5616  Nonstationary 
Differenced q3  Constant  5.35931  7.847e005  Stationary 
Q4  Constant, trend  2.75285  0.2152  Nonstationary 
Differenced q4  Constant  4.90132  3.236e005  Stationary 
3.3.2. Model Identification
For VAR model to be identified, the data should be stationary. We then determine the true lag order for the model. Selecting a higher order lag length than the true lag lengths increases the mean square forecast errors of the VAR, and selecting a lower order lag length than the true lag lengths usually causes auto correlated errors. This was pointed out by Lutkepohl. Hence, accuracy of forecasts from VAR models highly depends on selecting the true lag lengths. There are several statistical information criterions for selecting a lag length. This include: Akaike information criterion (AIC), Bayesian information criterion (BIC), and Hannanquinn information criterion (HQC). The table below shows the VAR model lag order selection criteria.
Lag  Aic  Bic  Hqc 
1  1.841568  2.730338*  2.148371 
2  1.493935*  3.093722  2.046181* 
3  1.650475  3.961277  2.448163 
The AIC and HQC is minimized at p=2. But BIC is minimized at p=1. We take BIC since it is appropriate for data with sample less than 120 quarters. Therefore, we have 20 parameters to estimate.
Δq1  Δq2  Δq3  Δq4  
Constant  0.0176441  0.0135887  0.00790380  0.0443565 
Δq1(1)  0.0230364  0.232444  0.153521  0.135078 
Δq2(1)  0.259410  0.501622  0.361898  0.853679 
Δq3(1)  0.556027  0.888681  1.88762  1.90887 
Δq4(1)  1.53304  1.74647  1.77344  1.27800 
The VAR(1) can be expressed as follows:
ΔQ1_{t}=0.0176441+0.0230364*ΔQ1_{t1}0.259410*ΔQ2_{t2}0.556027*ΔQ3_{t3}+1.53304*ΔQ4_{t4}
ΔQ2_{t}=0.0135887+0.0232444*ΔQ1_{t1}0.501622*ΔQ2t20.888681*ΔQ3t3+1.74647*ΔQ4_{t4}
ΔQ3_{t}=0.00790380+0.153521*ΔQ1_{t1}+0.361898*ΔQ2_{t2}1.88762*ΔQ3_{t3}+1.77344*ΔQ4_{t4}
ΔQ4_{t}=0.04435650.135078*ΔQ1_{t1}+0.853679*ΔQ2_{t2}1.90887*ΔQ3_{t3}+1.27800*ΔQ4_{t4}
3.3.3. Model Diagnostic Checking
We test for the autocorrelation of the model, to check whether the residuals of the fitted model are identically distributed. Table 9 below shows that the pvalue of all the four equations are greater than critical value at 1% hence we do not reject the null hypothesis of no autocorrelation, meaning there is consistent estimators.
Ljungbox q  Pvalue  
Equation 1  1.77125  0.183 
Equation 2  0.657489  0.417 
Equation 3  5.18819e005  0.994 
Equation 4  0.140521  0.708 
We also test for the auto effect correlation using arch. The pvalue is found to be greater than 0.01. Meaning we fail to reject the null hypothesis that there is no auto effect. This implies that there is conditional homoscedasticity which allows for a valid inference. This is shown in table 10 below;
Test statistic  Pvalue  
Equation 1  4.10135  0.0428491 
Equation 2  0.0950886  0.757805 
Equation 3  0.746557  0.387568 
Equation 4  0.0531222  0.817717 
The pvalue is greater than 0.01. Hence we fail to reject the null hypothesis that the residuals of the VAR model are normally distributed. Hence we can make an inference.
Residual correlation matrix, c (4 x 4)
1.0000 0.84237 0.70092 0.53609
0.84237 1.0000 0.93270 0.69950
0.70092 0.93270 1.0000 0.87207
0.53609 0.69950 0.87207 1.0000
Doornikhansen test, Chisquare (8) = 9.26726 [0.3203]
Testing all the residuals, the residuals of the VAR model have the pvalue of 0.3203 greater than 0.01, shows that all the residuals are normally distributed with the VAR model fitted being valid.
3.3.4. Prediction and Results
After obtaining the valid model, we now check the forecasting performance measures of each variable and compare with the actual value. The forecasting results of VAR1 model for 1998 are shown in table 11 below.
Variable  Actual value  Prediction value  Error 
Δq1  4.6667  4.5916  0.0751 
Δq2  4.43333  4.5905  0.1572 
Δq3  4.5000  4.6554  0.1554 
Δq4  4.4333  4.7524  0.3191 
3.3.5. Performance Comparison
Using the performance measures mentioned above, the table below shows the summary comparison performance of the two models.
Variable  Actual value  Prediction (SARIMA)  Prediction (var)  Error (SARIMA)  Error (Var) 
Q1  4.6667  4.3456  4.5916  0.3211  0.0751 
Q2  4.4333  4.7373  4.5905  0.304  0.1572 
Q3  4.5000  4.0605  4.6554  0.4395  0.1554 
Q4  4.4333  4.7509  4.7524  0.3176  0.3191 
MsE  0.1224  0.0391  
RmsE  0.3498  0.1977  
MpE  2.646%  3.147%  
MapE  7.662%  3.951%  
Theil's u  0.3498  0.0216 
3.4. Cumulative Results
The purpose for this project is to examine if VAR model is appropriate in forecasting a reshaped univariate seasonal time series. SARIMA model has been used for comparison purposes. Macro dataset used by stock and Watson in their introduction to econometrics where all unemployment rates for 16 years and above data were used. The data was quarterly and was reshaped into four variables and apply VAR for forecasting. This was analyzed in detail and the summary of other data set one was shown. The data contains observations from 1959 to 1999 where observations from 1959 to 1997 were used for model fitting. The year 1998 was used for prediction and comparison purposes.
Taking the detailed sample given above as an example, the performance of the VAR model is found to be better than the SARIMA model with the performance measures being minimum. The mean absolute percentage error for SARIMA is 7.662% and mean absolute percentage error for VAR is 3.951%. The other performance measures also show that VAR is a better model for forecast. Theil's u statistic for VAR is almost close to 0, implying almost a perfect fit.
However, in practice we normally face the challenge of the observations of seasonal univariate time series data that can be reshaped. This becomes a challenge in reshaping the data. Also fitting the reshaped data set into a VAR model is a bit complicated than the SARIMA model since each of the reshaped variables is nonstationary. However, VAR model is still evident to be the best model for all the observations than the SARIMA model.
References