Example: Error Correction Model in TSMT
The difference equation for a stationary VAR(P) model can be written as (Hamilton, 1994):Δyt = ζ1 Δyt-1 + ζ2 Δyt-2 + ... + ζp-1 Δyt-p+1 + ζ0 yt-1 + εt
where:ζ0 ≡ -(Ιn - Φ1 - Φ2 - ... Φp)
The error-correction form of a cointegrated VAR(P) system can be represented by:Δyt = ζ1 Δyt-1 + ζ2 Δyt-2 + ... + ζp-1 Δyt-p+1 + α - BA'yt-1 + εt
For a model with h cointegrating relationships, we can define:
zt ≡ A'yt-1 (a stationary hx1 vector).
we can rewrite the equation as:
Δyt = ζ1 Δyt-1 + ζ2 Δyt-2 + ... + ζp-1 Δyt-p+1 + α - Bzt-1 + εt
The time series module provides a number of routines for performing pre-estimation data analysis, model parameter estimation, and post-estimation diagnosis of cointegrated time series.Load and plot time series
For this example, we will use a data set that comes with TSMT, named rpie.xlsx. This file contains data for seasonally adjusted real personal income and expenditures, quarterly. 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/rpie.xlsx"; //Load data from the second and third columns //of a space separated text file y = xlsReadM(fname, "B2:C277"); //Load first date first_date = xlsReadM(fname, "A2:A2"); //Declare 'myPlot' to be a plotControl structure //and fill with default settings for XY plots struct plotControl myPlot; myPlot = plotGetDefaults("xy"); //Set text and location for legend plotSetLegend(&myPlot, "Consumption" $| "Income", "top left inside"); //Place x-tics every 10 years (4 quarters/year * 10 years = 40) //and place the first tic at 1950 plotSetXTicInterval(&myPlot, 40, 1950); //Only print 4 digit year at x-tic labels plotSetXTicLabel(&myPlot, "YYYY"); //Plot data using settings applied to 'myPlot' plotTS(myPlot, first_date, 1, y);
The above code should produce a plot that looks similar to this:
The two series definitely seem to have a relationship. Their change over time seems to be superlinear. In this tutorial, we will model the logged percent change of these variables. We can make this data transformation with the following code:
//Make TSMT library functions available library tsmt; //Transform data by taking first difference of each //of our variables after taking their natural log y_chg = vmdiffmt(ln(y), 1);
Let's now plot our transformed data.
//Continuing with 'myPlot' structure made above //Add a plot title and set its: text, font and font size plotSetTitle(&myPlot, "Log percent change", "Helvetica", 18); //Plot transformed data plotTS(myPlot, first_date, 4, y_chg);
Which should produce a graph similar to this:
This graph looks like what we would expect. The variance is not constant and any wild swings are corrected back towards normal quickly. The only problem here is the scale of the data. It appears that we have only one or, maybe two, data points with an absolute value greater than 0.05. The function we will use to estimate the parameters of this error-correction model is a Full Information Maximum Likelihood (FIML) method. The numerical optimization procedures can be sensitive to scaling problems with the data. In order to improve the scale of our data for the estimation, we will multiply it by 100 so that the range of our data will be from about -5 to 5, rather than -0.05 to 0.05.
Estimating the Error Correction Model
The ecmmt function is a convenient tool for estimating the parameters of ECM models with or without exogenous variables. The ecmmt function takes two required arguments:- A varmamtControl structure.
- An N x K data vector, y, where each column of y is a different time series.
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 'ctl' to be a varmamtControl structure //and fill with default values struct varmamtControl ctl; ctl = varmamtControlCreate(); //Declare 'out' to be an varmamtOut structure, //to hold the results from 'ecmmt' struct varmamtOut out; //No constraints on the optimization problem--for simplicity ctl.setConstraints = 0; //Specify the AR order to 2 ctl.ar = 2; //No constant in model ctl.const = 0;
//Scale data y_chg = 100 * y_chg; //Perform estimation with 'ecmmt' out = ecmmt(ctl, y_chg);
Δyt = ζ1 Δyt-1 + ζ2 Δyt-2 + ... + ζp-1 Δyt-p+1 - Πyt-1 + εt
===================================================== ECM Version 2.1.5 ===================================================== Residual Covariance Matrix 150.73574 73.83026 73.83026 249.96929 Zeta Plane [1,.,.] -0.29990 0.10971 -0.24282 0.12436 Plane [2,.,.] -0.14146 0.08794 -0.15634 0.18373 Pi 1.33370 -0.15140 0.00141 0.90544 Augmented Dickey-Fuller UNIT ROOT Test for Y1 Critical Values ADF Stat 1% 5% 10% 90% 95% 99% No Intercept -3.6411 -2.5897 -1.9567 -1.6179 0.9025 1.3221 2.0538 Intercept -7.0128 -3.4404 -2.8697 -2.5829 -0.4516 -0.0789 0.6319 Intercept and Time Trend -7.2323 -3.9542 -3.4305 -3.1392 -1.2333 -0.9390 -0.3755 Augmented Dickey-Fuller UNIT ROOT Test for Y2 Critical Values ADF Stat 1% 5% 10% 90% 95% 99% No Intercept -4.6936 -2.5897 -1.9567 -1.6179 0.9025 1.3221 2.0538 Intercept -8.9575 -3.4404 -2.8697 -2.5829 -0.4516 -0.0789 0.6319 Intercept and Time Trend -9.5108 -3.9542 -3.4305 -3.1392 -1.2333 -0.9390 -0.3755 Phillips-Perron UNIT ROOT Test for Y1 PPt 1% 5% No Intercept -9.7015 -2.6349 -1.9678 Intercept -15.5683 -3.4456 -2.8418 Intercept and Time Trend -15.7258 -3.9914 -3.4154 Phillips-Perron UNIT ROOT Test for Y2 PPt 1% 5% No Intercept -12.4127 -2.6349 -1.9678 Intercept -17.2373 -3.4456 -2.8418 Intercept and Time Trend -17.7911 -3.9914 -3.4154 Augmented Dickey-Fuller COINTEGRATION Test for Y1 Y2 Critical Values ADF Stat 1% 5% 10% 90% 95% 99% No Intercept -10.4949 -3.3567 -2.7752 -2.4659 -0.2541 0.1961 1.0722 Intercept -7.3543 -3.9243 -3.3804 -3.0821 -1.0076 -0.6342 0.0946 Intercept and Time Trend -7.4509 -4.4817 -3.8340 -3.5454 -1.6020 -1.3184 -0.7343 Johansen's Trace and Maximum Eigenvalue Statistics. r = # of CI Equations Critical Values r Trace Max. Eig 1% 5% 10% 90% No Intercept 0 170.8718 153.6202 1 17.2516 17.2516 1.0524 1.7046 2.1927 9.3918 Intercept 0 213.8515 156.7800 1 57.0716 57.0716 2.2515 3.3599 4.0975 12.8635 Intercept and Time Trend 0 221.3177 160.5904 1 60.7273 60.7273 4.0389 5.3796 6.1879 16.1762 Dep. Variable(s) : D(Y1 ) D(Y2 ) No. of Observations : 275 275 Degrees of Freedom : 260 260 Mean of Y : 0.8163 0.8157 Std. Dev. of Y : 0.8202 1.0107 Y Sum of Squares : 184.3448 279.8773 SSE : 150.7957 250.0293 MSE : 0.5637 0.9347 sqrt(MSE) : 0.7508 0.9668 R-Squared : 0.1820 0.1066 Adjusted R-Squared : 0.1379 0.0585 Model Selection (Information) Criteria ...................................... Likelihood Function : -1933.9647 Akaike AIC : 3837.9293 Schwarz BIC : 3952.1809 Likelihood Ratio : 3867.9293 Characteristic Equation(s) for Stationarity and Invertibility AR Roots and Moduli: Real : 1.95 1.23 1.23 Imag.: 0.00 1.63 -1.63 Mod. : 1.95 2.04 2.04 MULTIVARIATE ACF LAG01 LAG02 LAG03 0.008238 -0.04828 0.02109 0.009994 0.05089 0.05083 0.01015 -0.0364 0.002671 -0.009055 0.04582 -0.03054 LAG04 LAG05 LAG06 -0.06684 0.09705 -0.1203 -0.09942 -0.07708 0.01643 -0.06681 -0.1699 -0.09808 -0.0361 0.01039 0.06907 LAG07 LAG08 LAG09 0.05776 0.04065 -0.1715 -0.07736 0.1374 0.08797 0.1593 0.05391 -0.03936 -0.0111 0.08368 0.01126 LAG10 LAG11 LAG12 0.08925 0.06825 0.03497 0.08071 -0.02322 -0.006081 -0.006279 0.1147 0.03267 0.05381 0.0337 0.005717 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 3.8602 0.4253 4 24.4700 0.0019 5 30.5424 0.0023 6 34.2922 0.0050 7 41.6454 0.0031 8 50.2697 0.0013 9 56.8432 0.0010 10 64.2697 0.0006 11 66.3653 0.0015 12 67.1565 0.0046
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 view original data ** ** ** *******************************************/ //Get path to GAUSSHOME directory path = getGAUSSHome(); //Add file name to path fname = path $+ "examples/rpie.xlsx"; //Load data from the second and third columns //of a space separated text file y = xlsReadM(fname, "B2:C277"); //Load first date first_date = xlsReadM(fname, "A2:A2"); //Declare 'myPlot' to be a plotControl structure //and fill with default settings for XY plots struct plotControl myPlot; myPlot = plotGetDefaults("xy"); //Set text and location for legend plotSetLegend(&myPlot, "Consumption" $| "Income", "top left inside"); //Place x-tics every 10 years (4 quarters/year * 10 years = 40) //and place the first tic at 1950 plotSetXTicInterval(&myPlot, 40, 1950); //Only print 4 digit year at x-tic labels plotSetXTicLabel(&myPlot, "YYYY"); //Plot data using settings applied to 'myPlot' plotTS(myPlot, first_date, 4, y); /******************************************* ** ** ** Transform data ** ** ** *******************************************/ //Transform both columns (variables) in one step y = 100 * vmdiffmt(ln(y), 1); /******************************************* ** ** ** Plot transformed data ** ** ** *******************************************/ //Create new graph tab so we do not //overwrite our last plot plotOpenWindow(); //Continuing with 'myPlot' structure made above //Add a plot title and set its: text, font and font size plotSetTitle(&myPlot, "Log percent change", "Helvetica", 18); //Plot transformed data plotTS(myPlot, first_date, 4, y); /******************************************* ** ** ** Estimate ** ** ** *******************************************/ //Declare 'ctl' to be a varmamtControl structure struct varmamtControl ctl; //Fill 'ctl' with default values ctl = varmamtControlCreate(); //Specify a VAR(2) model ctl.ar = 2; //Do not set constraints for internal optimization model ctl.setConstraints = 0; //Declare 'out' to be a varmamtOut structure //to hold the return values from 'out' struct varmamtOut out; //Estimate the parameters of the model //and print diagnostic information out = ecmmt(ctl, y);