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.0046The 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);