Time-series Analysis with VAR & VECM: Statistical approach with complete Python code (2022)

A case study with stock market data

Modeling Time-series Stochastic Data

https://sarit-maitra.medium.com/membership

VECTOR auto-regressive (VAR) integrated model comprises multiple time series and is quite a useful tool for forecasting. It can be considered an extension of the auto-regressive (AR part of ARIMA) model. VAR model involves multiple independent variables and therefore has more than one equations. Each equation uses as its explanatory variables lags of all the variables and likely a deterministic trend. Time series models for VAR are usually based on applying VAR to stationary series with first differences to original series and because of that, there is always a possibility of loss of information about the relationship among integrated series.

Therefore, differencing the series to make them stationary is one solution, but at the cost of ignoring possibly important (“long run”) relationships between the levels. A better solution is to test whether the levels regressions are trustworthy (“cointegration”.) The usual approach is to use Johansen’s method for testing whether or not cointegration exists. If the answer is “yes” then a vector error correction model (VECM), which combines levels and differences, can be estimated instead of a VAR in levels. So, we shall check if VECM is been able to outperform VAR for the series we have.

This an extension of my previously published article.

Loading all datasets ( Gold, Silver & Crude Oil)

After necessary cleaning & pre-processing (filling the missing values with previous ones), we finally have the three time series for necessary analysis.

A quick test is to check if the data is random. Random data will not exhibit a structure in the lag plot.

The time series plot clearly indicates some kind of relationships among the series. The linear shape of the lag plot suggests that an AR model is a better choice. We also don’t see any outlier in the data. Data here showing linear pattern, indicating the presence of positive auto-correlation.

“Time series for economic data is generally stochastic or has a trend that is not stationary, meaning that the data has a root unit”

# plots the autocorrelation plots at 75 lags
for i in dataset:
 plot_acf(dataset[i], lags = 50)
 plt.title(‘ACF for %s’ % i) 
 plt.show()

Stationarity check

def augmented_dickey_fuller_statistics(time_series):
 result = adfuller(time_series.values)
 print('ADF Statistic: %f' % result[0])
 print('p-value: %f' % result[1])
 print('Critical Values:')
for key, value in result[4].items():
 print('\t%s: %.3f' % (key, value))

Through the above function, we can run Augmented Dickey Fuller (ADF) test on all columns, which clearly shows the original series are non-stationary and contain unit root.

print('Augmented Dickey-Fuller Test: Gold Price Time Series')
augmented_dickey_fuller_statistics(X_train['Gold'])
print('Augmented Dickey-Fuller Test: Silver Price Time Series')
augmented_dickey_fuller_statistics(X_train['Silver'])print('Augmente
d Dickey-Fuller Test: Oil Price Time Series')
augmented_dickey_fuller_statistics(X_train['Oil'])

VAR models can also be used for analyzing the relation between the variables involved using Granger Causality tests. Granger causality specifies that a variable y1t is causal for a variable y2t if the information in y1t is helpful for improving the forecasts of y2t.

Granger Causality tests try to determine if one variable(x1) can be used as a predictor of another variable(x2) where the past values of that another variable may or may not help. This means that x1 explains beyond the past values of x2. Two important assumptions here are -

• both x1 and x2 are stationary
• there exists a linear relation between their current and past values.

This means that if x1 and x2 are non-stationary, we have to make them stationary before testing for Granger Causality.

Split the Series into Training and Testing Data

We will fit the VAR model on X_train to forecast the next 10 observations. These forecasts will be compared against the actuals present in test data (X_test). We shall use multiple forecast accuracy metrics.

Series transformation

A difference transform is a simple way for removing a systematic structure from the time series. We will remove trend by subtracting the previous value from each value in the series which is the first order differencing. To keep it simple, we will do first order differencing or seasonal differencing.

If we have an integrated order to n time series and if we take first order to difference and time, we will be left with series integrated order of zero.

X_train_log = np.log(X_train)
X_train_log_diff =(X_train_log).diff().dropna()
X_train_log_diff.describe()

looking at the plot, we could figure out that data set looks like normalized

Auto-Correlation Function Analysis of Transformed Series

Below function plots the auto-correlation plots for the difference in each stock’s price from the price the previous trading day at 75 lags.

fig, ax = plt.subplots(1,2, figsize=(10,5)) 
ax[0] = plot_acf(X_train_log_diff['Gold'], ax=ax[0])
ax[1] = plot_pacf(X_train_log_diff['Gold'], ax=ax[1])

We have shown ACF & PACF of transformed Gold series; likewise, other series can be plotted.

ADF Test - transformed series

print('Augmented Dickey-Fuller Test: Gold Price Time Series')
augmented_dickey_fuller_statistics(X_train_log_diff['Gold'])
print('Augmented Dickey-Fuller Test: Silver Price Time Series')
augmented_dickey_fuller_statistics(X_train_log_diff['Silver'])
print('Augmented Dickey-Fuller Test: Oil Price Time Series')
augmented_dickey_fuller_statistics(X_train_log_diff['Oil'])
Augmented Dickey-Fuller Test on "Oil" 

The Granger causality test is conducted to determine whether one time series is useful in forecasting another. A time series X is said to Granger-cause Y if it can be shown, usually through a series of t-tests and F-tests on lagged values of X (and with lagged values of Y also included), that those X values provide statistically significant information about future values of Y.

Multivariate analysis

Here, for multivariate Granger causality analysis performed by fitting a VAR to the time series. Considering below is a d-dimensional multivariate time series —

Granger causality is performed by fitting a VAR model with L time lags as follows:

where ε ( t ) is a white Gaussian random vector, and A τ is a matrix for every τ. A time series X i is called a Granger cause of another time series X j, if at least one of the elements A τ ( j , i ) for τ = 1 , … , L is significantly larger than zero.

print(grangercausalitytests(X_train_log_diff[['Gold','Silver']], maxlag=15, addconst=True, verbose=True))
print(grangercausalitytests(X_train_log_diff[['Gold','Oil']], maxlag=15, addconst=True, verbose=True))
print(grangercausalitytests(X_train_log_diff[['Oil','Silver']], maxlag=15, addconst=True, verbose=True))

Below output shown for Gold & Oil which differs the test hypothesis till lag 4.

A VAR(p) process in its basic form is:

Here yt​​ represents a set of variables collected in a vector, c denotes a vector of constants, a is a matrix of autoregressive coefficients and e​t​​ is white noise. Since the parameters of a​​ are unknown, we have to estimate these parameters. Each variable in the model has one equation. The current (time t) observation of each variable depends on its own lagged values as well as on the lagged values of each other variable in the VAR.

I have implemented Akaike’s Information Criteria (AIC) through the VAR (p) to determine the lag order value. In the fit function, I have passed a maximum number of lags and the order criterion to use for order selection.

#Initiate VAR model
model = VAR(endog=X_train_log_diff)
res = model.select_order(15)
res.summary()

#Fit to a VAR model
model_fit = model.fit(maxlags=3)
#Print a summary of the model results
model_fit.summary()

Forecast VAR Model

The forecasts are generated on the training data used by the model.

# Get the lag order
lag_order = model_fit.k_ar
print(lag_order)# Input data for forecasting
input_data = X_train_log_diff.values[-lag_order:]
print(input_data)# forecasting
pred = model_fit.forecast(y=input_data, steps=nobs)
pred = (pd.DataFrame(pred, index=X_test.index, columns=X_test.columns + '_pred'))
print(pred)

Invert transformation

So, to bring it back up to its original scale, we need to de-difference to the original input data. Our data is 1st logarithm transformed and then differenced. So, to inverse, we have to first use cumulative sum to de-differentiate and then use exponential. Natural logarithm is the inverse of the exp().

# inverting transformation
def invert_transformation(X_train, pred_df):
 forecast = pred.copy()
 columns = X_train.columns
for col in columns:
 forecast[str(col)+'_pred'] = X_train[col].iloc[-1] + forecast[str(col) +'_pred'].cumsum()
 return forecastoutput = invert_transformation(X_train, pred)
print(output)output_original = np.exp(output)
print(output_original)

VAR Forecast evaluation

#Calculate forecast bias
forecast_errors = [X_test['Oil'][i]- output_original['Oil_pred'][i] for i in range(len(X_test['Oil']))]
bias = sum(forecast_errors) * 1.0/len(X_test['Oil'])
print('Bias: %f' % bias)#Calculate mean absolute error
mae = mean_absolute_error(X_test['Oil'],output_original['Oil_pred'])
print('MAE: %f' % mae)#Calculate mean squared error and root mean squared error
mse = mean_squared_error(X_test['Oil'], output_original['Oil_pred'])
print('MSE: %f' % mse)
rmse = sqrt(mse)
print('RMSE: %f' % rmse)

“Least squares parameter estimation of dynamic regression models is known to exhibit substantial bias in small samples when the data is fairly persistent”

See Also

SUSTENTACIONES VIRTUALES DE LA UNIVERSIDAD RICARDO PALMALet's talk about Rosa Var Attre, the impossible romance of The Witcher 3 - Kill Screen - PreviouslyJapanese Boxwood Plant Care Guide (Buxus microphylla) | Home for the Harvest

VECM imposes additional restriction due to the existence of non-stationary but co-integrated data forms. It utilizes the co-integration restriction information into its specifications. After the cointegration is known then the next test process is done by using error correction method. Through VECM we can interpret long term and short term equations. We need to determine the number of co-integrating relationships. The advantage of VECM over VAR is that the resulting VAR from VECM representation has more efficient coefficient estimates.

Johansen Cointegration

In order to fit a VECM model, we need to determine the number of co-integrating relationships using a VEC rank test.

vec_rank1 = vecm.select_coint_rank(X_train, det_order = 1, k_ar_diff = 1, method = 'trace', signif=0.01)
print(vec_rank.summary())

We find the λtrace statistics in the third column, together with the corresponding critical values. The test statistic of 38.25 is lower than the critical value (41.08) and so the null of at most one co-integrating vector cannot be rejected.

Let us employ an alternative statistic, the maximum-eigenvalue statistic (λmax).

vec_rank2 = vecm.select_coint_rank(X_train, det_order = 1, k_ar_diff = 1, method = 'maxeig', signif=0.01)
print(vec_rank2.summary())

The test output reports the results for the λmax statistics which does not differ much from trace statistic; the critical value (29.28) is still higher than test statistic.

We will still go ahead and estimate VECM, since it can still valuable for short-run dynamics in absence of co-integration. Let’s estimates the VECM on the prices with 9 lags, 1 co-integrating relationship, and a constant within the co-integration relationship. I have used ‘cili’ a combination of “ci” — constant within the co-integration relation and “li” — linear trend within the co-integration relation

vecm = VECM(endog = X_train, k_ar_diff = 9, coint_rank = 3, deterministic = ‘ci’)
vecm_fit = vecm.fit()
vecm_fit.predict(steps=10)

forecast, lower, upper = vecm_fit.predict(10, 0.05)
print(“lower bounds of confidence intervals:”)
print(lower.round(3))
print(“\npoint forecasts:”)
print(forecast.round(3))
print(“\nupper bounds of confidence intervals:”)
print(upper.round(3))

VECM Forecast evaluation

Though we had an indication that, VAR would be best for our data set for price prediction; however, we have shown VECM for experimentation and illustration purpose. This is a simple procedure to explain VAR. However, there are other procedures like impulse response analysis and variance decomposition also can be introduced to experiment if we are able to see how a shock to one variable affects other variable in subsequent periods.

Time series for economic data is generally stochastic or has a trend that is not stationary, meaning that the data has a root unit. To be able to estimate a model using the data-

Steps for VAR-

1. Test stationarity of data and degree of integration
2. Determination of lag length
3. Test the granger causality
4. Estimation of VAR
5. Variance decomposition

Forecasting Steps for VECM-

1. Determination of lag length
2. Test the granger causality
3. Cointegration degree test
4. Estimation of VECM
5. Variance decomposition

I can be reached here.

Notice: The programs described here are experimental and should be used with caution. All such use at your own risk.

References:

(1) Rao, B. (2007). Cointegration: for the Applied Economist, Springer.

(2) Ashley, R. A., & Verbrugge, R. J. (2009). To difference or not to difference: a Monte Carlo investigation of inference in vector autoregression models. International Journal of Data Analysis Techniques and Strategies, 1(3), 242–274.

(3) Lütkepohl, H. (2011). Vector autoregressive models. In International Encyclopedia of Statistical Science (pp. 1645–1647). Springer Berlin Heidelberg.

(4) Kuo, C. Y. (2016). Does the vector error correction model perform better than others in forecasting stock price? An application of residual income svaluation theory. Economic Modelling, 52, 772–789.