VAR Model

Example: VAR Model in TSMT

Vector autoregressive models generalize the single variable autoregressive model to allow for more than one endogenous variable. The standard (or reduced) form of a VAR(1) model with two separate time series is:

yt,1 = α1 + φ11 yt-1,11,2 yt-1,2 + ut,1
yt,2 = α2 + φ21 yt-1,122 yt-1,2 + ut,2

Assuming:

  1. y is an N by p matrix, where N is the length of the time series and p is the order of the VAR.
  2. alpha is a 1 by p matrix, representing the intercept term.
  3. e is an N by p matrix containing the error term.

We can express a VAR(p) model in matrix form, using GAUSS code as:

y[t,1:p] = alpha + y[t-1,.] * phi' + e[t,.];
The time series module provides a number of routines for performing pre-estimation data analysis, model parameter estimation, and post-estimation diagnosis of autoregressive time series.

Load and plot time series

For this example, we will use a data set that comes with TSMT, named minkmt.asc. This file contains data for the prices of two commodities from 1850 to 1910. The first column contains the dates and the second and third columns contain the data.

//Get path to GAUSSHOME directory
path = getGAUSSHome();

//Add file name to path
fname = path $+ "examples/minkmt.asc";

//Load data from the second and third columns
//of a space separated text file
y = csvReadM(fname, 1, 2|3, " ");

//Load first date
first_date = csvReadM(fname, 1|1, 1|1, " ");

//Plot data
plotTS(first_date, 1, y);

The above code should produce a plot that looks similar to this:

Time series plot of 2 commodities

The first thing that we notice is that the, at least the first series, seems to have a trend and that neither series has a mean near zero. We can use the TSMT function, vmdetrendmt to remove the trend and demean the data. The vmdetrendmt function takes two inputs:

  1. An NxK matrix containing one or more time series to be acted upon.
  2. A scalar, indicator variable with the following options:
    0Return the original data untouched.
    1Detrend the data, based upon a linear model.
    2+Detrend the data based upon a polynomial model.

Now we will detrend and demean our data, assuming only a linear trend.

//Load TSMT library functions
library tsmt;

//Detrend and demean data
y_detrend = vmdetrendmt(y, 1);

//Plot transformed data
plotTS(first_date, 1, y_detrend);

The graph of the detrended data should look more like this:

Plot of two time series after demeaning and detrending.

After detrending, the data appears to have a mean of zero and does not have a trend. However, the variance of the data seems to increase with time. Let's try differencing the data instead of detrending:

//Calculate the first difference
y_diff = vmdiffmt(y, 1);

//Plot transformed data
plotTS(first_date, 1, y_diff);

The plot of the differenced data should look like this:

Plot of two, once differenced, time series

This plot does not appear to have a trend, the mean is near zero and the variance does not appear to increase with time.

Estimating the Vector Autoregressive Model

The varmaxmt function is a convenient tool for estimating the parameters of VAR models with or without exogenous variables.  The varmaxmt function takes two required arguments:
  1. A varmamtControl structure.
  2. An N x K data vector, y, where each column of y is a different time series.
Note that varmaxmt will produce an error if the data is not stationary and invertible.

The varmamtControl Structure

This structure handles the matrices and strings that control the estimation process such as specifying whether to include a constant, output appearance, and estimation parameters.  Continuing with the example above,
//Declare 'vctl' to be a varmamtControl structure
//and fill with default values
struct varmamtControl vctl;
vctl = varmamtControlCreate();

//Declare 'vout' to be an varmamtOut structure,
//to hold the results from 'varmaxmt'
struct varmamtOut vout;

//Specify the AR order to 2
vctl.ar = 2;

//Specify order of differencing
vctl.diff = 1;

//Specify MA order
vctl.ma = 0;

//No constant in model
vctl.const = 0;
This appropriately passes to varmaxmt that the model is an VAR(2) model and does not include a constant.  The remaining control elements (the ones that we did not set above), left at the default values, specify preferences for output display and estimation. With the control structure options set, we can now perform the estimation:
//Perform estimation with 'varmaxmt'
//Note that we are passing in the original, 'y'
//Since we told 'varmaxmt' to difference 'y',
//by setting vctl.diff = 1
vout = varmaxmt(vctl, y);
which will return output similar to:
===============================================
 VARMAX Version 2.1.3                          
===============================================


Phi


Plane [1,.,.] 

  -0.12309    0.57672 
  -0.67611    0.34425 

Plane [2,.,.] 

   0.26190   -0.12894 
  -0.26000  -0.044378 

standard errors


Plane [1,.,.] 

    1.0000     1.0000 
    1.0000     1.0000 

Plane [2,.,.] 

    1.0000     1.0000 
    1.0000     1.0000 

t-values


Plane [1,.,.] 

  -0.12309    0.57672 
  -0.67611    0.34425 

Plane [2,.,.] 

   0.26190   -0.12894 
  -0.26000  -0.044378 

probabilities


Plane [1,.,.] 

   0.90253    0.56672 
   0.50209    0.73210 

Plane [2,.,.] 

   0.79447    0.89792 
   0.79593    0.96478 


Residual Covariance Matrix

  0.066376   0.026428 
  0.026428   0.065819 

standard errors

    1.0000     1.0000 
    1.0000     1.0000 

t-values

  0.066376   0.026428 
  0.026428   0.065819 

probabilities

   0.94734    0.97902 
   0.97902    0.94778 


Beta0

 -0.010973 
  0.021210 

Augmented Dickey-Fuller UNIT ROOT Test for    Y1
                                                          Critical Values
                               ADF Stat        1%       5%       10%      90%      95%     99%
   No Intercept                -5.2278     -2.5328  -1.9498  -1.6266   0.9152   1.3168   2.1179
   Intercept                   -5.1774     -3.5663  -2.9370  -2.6152  -0.4393  -0.0499   0.6942
   Intercept and Time Trend    -5.2981     -4.0892  -3.4615  -3.1709  -1.2584  -0.9195  -0.2986

Augmented Dickey-Fuller UNIT ROOT Test for    Y2
                                                          Critical Values
                               ADF Stat        1%       5%       10%      90%      95%     99%
   No Intercept                -4.7720     -2.5328  -1.9498  -1.6266   0.9152   1.3168   2.1179
   Intercept                   -4.7352     -3.5663  -2.9370  -2.6152  -0.4393  -0.0499   0.6942
   Intercept and Time Trend    -4.6155     -4.0892  -3.4615  -3.1709  -1.2584  -0.9195  -0.2986


Phillips-Perron UNIT ROOT Test for    Y1
                                  PPt      1%      5%
   No Intercept                -7.2621  -2.5328  -1.9498
   Intercept                   -6.9186  -3.5663  -2.9370
   Intercept and Time Trend    -6.9497  -4.0892  -3.4615

Phillips-Perron UNIT ROOT Test for    Y2
                                  PPt      1%      5%
   No Intercept                -6.3499  -2.5328  -1.9498
   Intercept                   -5.9127  -3.5663  -2.9370
   Intercept and Time Trend    -5.8505  -4.0892  -3.4615


Augmented Dickey-Fuller COINTEGRATION Test for    Y1    Y2

                                                          Critical Values
                               ADF Stat        1%       5%       10%      90%      95%     99%
   No Intercept                -5.0448     -3.4003  -2.8198  -2.4901  -0.2841   0.1628   0.9912
   Intercept                   -5.0375     -4.0246  -3.4040  -3.0890  -0.9988  -0.6383   0.0929
   Intercept and Time Trend    -5.1691     -4.5041  -3.9157  -3.6062  -1.6464  -1.3413  -0.6750



Johansen's  Trace  and  Maximum  Eigenvalue  Statistics. r = # of CI Equations
                                                                      Critical Values
                               r    Trace    Max. Eig          1%       5%       10%      90%
   No Intercept                0   65.5906   46.3981
                               1   19.1925   19.1925         1.0524   1.7046   2.1927   9.3918

   Intercept                   0   65.6102   46.4140
                               1   19.1962   19.1962         2.2515   3.3599   4.0975  12.8635

   Intercept and Time Trend    0   66.2392   46.3986
                               1   19.8406   19.8406         4.0389   5.3796   6.1879  16.1762

Dep. Variable(s)    :         Y1                Y2       
No. of Observations :             61                61
Degrees of Freedom  :             50                50
Mean of Y           :         0.0018            0.0279
Std. Dev. of Y      :         0.3014            0.3506
Y Sum of Squares    :         5.4509            7.3772
SSE                 :         4.0489            4.0150
MSE                 :         0.0730            0.0723
sqrt(MSE)           :         0.2701            0.2690
R-Squared           :         0.2572            0.4558
Adjusted R-Squared  :         0.1086            0.3469


Model Selection (Information) Criteria
......................................
Likelihood Function :        -2.0830
Akaike AIC          :       -17.8339
Schwarz BIC         :        49.3857
Likelihood Ratio    :         4.1661




Characteristic Equation(s) for Stationarity and Invertibility
   AR Roots and Moduli:
   Real :         4.2121928       -1.9895302 
   Imag.:         0.0000000        0.0000000 
   Mod. :         4.2121928        1.9895302 



                                 MULTIVARIATE ACF

               LAG01                    LAG02                    LAG03
       0.03143    0.06429       -0.1224    -0.0856       -0.0715    -0.1648
      0.002152   -0.03561       -0.1163    -0.2181         0.133     0.1341


               LAG04                    LAG05                    LAG06
       -0.1864    -0.1613       -0.1483    0.09278       -0.3151    0.02871
      -0.04398   -0.04798       0.01561   -0.04506       -0.1768   -0.03803


               LAG07                    LAG08                    LAG09
      0.008986   -0.02272      -0.04262    0.01081        0.2493     0.3194
      -0.04619    -0.1493      -0.09121    0.06298        0.1083    0.08513


               LAG10                    LAG11                    LAG12
        0.4103    0.04377       0.01628    -0.1909      -0.01679   0.005056
        0.1636    -0.1654        0.1557     0.1854      -0.09051    0.07984


    ACF INDICATORS: SIGNIFICANCE = 0.95 
  (using Bartlett's large sample standard errors)

LAG01     LAG02     LAG03     LAG04     LAG05     LAG06     
-+-+      -+-+      -+-+      -+-+      -+-+      -+-+      
-+-+      -+-+      -+-+      -+-+      -+-+      -+-+      


LAG07     LAG08     LAG09     LAG10     LAG11     LAG12     
-+-+      -+-+      -+-+      -+-+      -+-+      -+-+      
-+-+      -+-+      -+-+      -+-+      -+-+      -+-+      


Multivariate Goodness of Fit Test
Lag       Qs     P-Value
3      9.5468    0.0488
4     12.4737    0.1313
5     17.3960    0.1353
6     26.3349    0.0495
7     28.0138    0.1091
8     29.4429    0.2039
9     38.2221    0.0943
10    55.2182    0.0066
11    64.9744    0.0022
12    66.9787    0.0048

Forecasting

We will produce forecasts from our model, using vmforecastmt. It takes five inputs:

  1. A varmamtControl structure which, as we saw above, contains the AR order and other model settings.
  2. A varmamtOut structure, possibly returned by varmaxmt, which will contain the parameter estimates for our VAR model.
  3. A time series vector from which to produce forecasts. The length of this vector must be (p + 1), where p is the specified order of the AR process. However, vmforecast will read from the end of the vector. This allows you to pass in the entire time series if you want forecasts from the end of your data.
  4. A vector or matrix, containing any exogenous variables in the model over the forecast horizon.
  5. A scalar, the number of periods to forecast.

vmforecastmt returns a matrix in which the first column contains the forecast observation number (from 1 to the number of forecast periods), the remaining columns contain the forecast for each of the endogenous variables. Continuing on with our example, we will produce a 10 period forecast:

//Produce 10 period forecast
//Since we have no exogenous variables,
//we pass a scalar zero for the 4th input
f = vmforecastmt(vctl, vout, y, 0, 10);

The forecast matrix, f should look similar to this:

 1            0.437            0.222 
 2           0.0291           -0.286 
 3            0.151           -0.338 
 4          -0.0829           -0.242 
 5           -0.169           -0.213 
 6          -0.0461          -0.0513 
 7          -0.0927           0.0731 
 8          -0.0407           0.0669 
 9           0.0481            0.102 
10          0.00988           0.0714

After calculating our forecast, we can plot the original differenced data and our forecasts together like this:

//Declare 'myPlot' to be a plotControl structure
struct plotControl myPlot;

//Fill 'myPlot' with default values
myPlot = plotGetDefaults("xy");

//Set line colors
plotSetLineColor(&myPlot, "orange"$|"blue");

plotSetTitle(&myPlot, "Differenced data and VAR(2) model forecast", "times", 18);

//Plot differenced 'y' data
plotTS(myPlot, 1850, 1, y);

//Set line style dots for forecast period
plotSetLineStyle(&myPlot, 3);

//Add forecasts to plot
plotAddTS(myPlot, 1910, 1, f[.,2:3]);

Which should produce the following graph.

Graph of two differenced time series and 10 period forecasts.

The Command File

Finally we put it all together in the command file below, which will reproduce all output shown above.
//Load TSMT library functions
library tsmt;

/*******************************************
**                                        **
**       Load and process data            **
**                                        **
*******************************************/

//Create file name with full path
fname = getGAUSSHome() $+ "examples/minkmt.asc";

//Load all rows of the second and third columns of
//a space separated file
y = csvReadM(fname, 1, 2, " ");

//Difference the data
y = vmdiffmt(y, 1);

/*******************************************
**                                        **
**             Estimate                   **
**                                        **
*******************************************/

//Declare 'vctl' to be a varmamtControl structure
struct varmamtControl vctl;

//Fill 'vctl' with default values
vctl = varmamtControlCreate();

//Specify a VAR(2) model
vctl.ar = 2;

//Declare 'vout' to be a varmamtOut structure
//to hold the return values from 'vout'
struct varmamtOut vout;

//Estimate the parameters of the VAR(2) model
//and print diagnostic information
vout = varmaxmt(vctl, y); 

/*******************************************
**                                        **
**              Forecast                  **
**                                        **
*******************************************/

//Create a 10 period forecast
f = vmforecastmt(vctl, vout, y, 0, 10);

/*******************************************
**                                        **
**   Plot differenced data and forecast   **
**                                        **
*******************************************/

//Declare 'myPlot' to be a plotControl structure
struct plotControl myPlot;

//Fill 'myPlot' with default values
myPlot = plotGetDefaults("xy");

//Set line colors
plotSetLineColor(&myPlot, "orange"$|"blue");

plotSetTitle(&myPlot, "Differenced data and VAR(2) model forecast", "times", 18);

//Plot differenced 'y' data
plotTS(myPlot, 1850, 1, y);
plotSetLineStyle(&myPlot, 3);

//Add forecasts to plot
plotAddTS(myPlot, 1910, 1, f[.,2:3]);

Have a Specific Question?

Get a real answer from a real person

Need Support?

Get help from our friendly experts.