Example: ARMA Model in TSMT
Autoregressive models take the general form:
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.
Preliminary Autocorrelations and Partial Autocorrelations
Suppose we have stationary univariate time series data but are uncertain about the order of the regressive relationship. You may wish to first use the sample autocorrelation function can help identify the general structure of the model. This can be achieved in Gauss using the autocormt and autocovmt functions.
The autocor function computes the autocorrelations for each column of a matrix. It assumes that all data has been demeaned and produces a output matrix containing the autocorrelations for each column of an input matrix. The autocor function requires three inputs: (1) an N x K datamatrix, x, (2) a scalar, f, denoting the first autocorrelation to be computed, and (3) a scalar, l, denoting the last autocorrelation to be computed.
Consider an autoregressive time series randomly generated using the Gauss ARMA data simulation function. For example, consider a time series that follows the purely autoregressive data generation process given by
Realizations of any ARMA data generation process can be randomly created using the function simarmamt ,
//SPECIFY ARMA PARAMETERS b_sim = {0.5,-0.8}; //SPECIFY MA ORDER q_sim = 0; //SPECIFY AR ORDER p_sim = 2; //SPECIFY DETERMINISTIC FEATURES (CONSTANT AND TREND) const_sim=1.5; //NUMBER OF OBSERVATIONS n_sim=200; //NUMBER OF SERIES k_sim=4; //STANDARD DEVIATION OF ERROR TERMS std_sim=1; seed_sim = 5012; y_sim = simarmamt(b_sim,p_sim,q_sim,const_sim,n_sim,k_sim,std_sim,seed_sim); //VISUALIZE DATA x_plot=seqa(1,1,rows(y_sim)); //PLOT Y OBS plotXY(x_plot,y_sim);
The code above will produce a 200X4 matrix of data that follows the AR(2) data generating process specified. Each column of the matrix represents a different random realization of the process. The code graphs the first realization:
Though the data generating process of the data above is known, it is generally more realistic the one will wish to estimate the parameters of an unknown data generating process behind a realization of data. To begin consider finding the ACF and the PACF the autocor and pacf Gauss functions ,
//FIND AUTOCORRELATIONS FOR 0 to 10 lags //DEMEAN DATA plotSetTitle(&myPlot, "Demeaned Data" ); y = y_sim - meanc(y_sim)'; plotXY(myPlot,x_plot,y_sim); //FIND AUTOCORRELATIONS a = autocor(y_sim, 0, 10); a_lab = seqa(0,1,rows(a)); format /rd 12,2; //DISPLAY AUTOCORRELATIONS let a_col_lab = { "Lags" "ACF" }; print /rd $a_col_lab;; print a_lab~a;
This code calculates and lists the sample autocorrelations for lags 0 through 10 for each sample of data. All columns of the list have a precision of two decimal points. Column labels “Lags” and “ACF” are included at the top of the list. This list will change depending on the data that has been randomly generated but should look similar to:
Lags ACF 0.00 1.00 1.00 1.00 1.00 1.00 0.28 0.30 0.28 0.25 2.00 -0.71 -0.64 -0.63 -0.68 3.00 -0.59 -0.61 -0.54 -0.47 4.00 0.31 0.15 0.22 0.31 5.00 0.64 0.58 0.54 0.48 6.00 0.05 0.24 0.11 -0.05 7.00 -0.50 -0.31 -0.37 -0.38 8.00 -0.25 -0.39 -0.30 -0.05 9.00 0.33 -0.01 0.08 0.27 10.00 0.38 0.32 0.28 0.11
As an alternative, a bar graph visually presents the autocorrelation information. For example,
//DISPLAY CORRELOGRAM //Initialize plotControl STRUCTURE struct plotControl myPlot; myPlot = plotGetDefaults("bar"); //ADD A TITLE plotSetTitle(&myPlot, "Correlogram" ); //ADD AXIS LABLES plotSetYLabel(&myPlot, "ACF"); plotSetXLabel(&myPlot, "Lag"); //PLOT ACF x_lag=seqa(0, 1, 11); plotBar(myPlot, x_lag, a[1:rows(a),1]);
This produces the following output:
The pacf function supplements the diagnostic information in the ACF and finds the partial autocorrelation function. The previous example is easily extended to find the PACF for the same randomly generated data. The pacf function requires three inputs: (1) an N x 1 data vector, x, (2) a scalar, k, denoting the maximum number of partial autocorrelation to compute, and (3) a scalar, d, denoting the order of differencing. Note that in GAUSS the pacf function always begins with the minimum of one lag. Hence, if k=10, pacf will find the partial autocorrelations for lags 1 through 10. For example,
//FIND PACF b = pacf(y[1:rows(y),1], 10, 0); format /rd 12,2; b_lab = seqa(0, 1, rows(a)); //DISPLAY AUTOCORRELATIONS let b_col_lab = { "Lags" "ACF" }; print /rd $b_col_lab;; print b_lab~a; //PLOT PACF //INITIALIZE plotControl STRUCTURE struct plotControl myPlot; myPlot = plotGetDefaults("bar"); //ADD A TITLE plotSetTitle(&myPlot, "Partial Autocorrelation Function"); //ADD AXIS LABLES plotSetYLabel(&myPlot, "PACF"); plotSetXLabel(&myPlot, "Lag"); //PLOT ACF plotBar(myPlot, b_lab,b);
This code calculates, lists, and graphs the partial autocorrelations for lags 0 through 10. This list and graph will change depending on the data that has been randomly generated but should look similar to:
Lags PACF 1.00 0.28 2.00 -0.85 3.00 0.05 4.00 0.00 5.00 -0.06 6.00 0.00 7.00 0.04 8.00 0.09 9.00 0.09 10.00 0.07
Estimating the Autoregressive Model
The arimamt function is a convenient tool for estimating the parameters of any ARIMA model, including ARMA models, purely AR models, and purely MA models. The arimamt function takes two arguments, (1) an arimamt control structure, (2) N x 1 data vector, y. Note that this function can only be used to estimate models with a specified number of sequential ordered lags. For example, it will estimate the parameters for the first through fourth lag but will not estimate a model that includes ONLY the first and fourth lag. In addition, note that this arimamt will produce an error if the data is not stationary and invertible.
The arimamt Control Structure
This structure handles the matrices and strings that control the estimation process such as setting specifying whether to include a constant, output appearance, and estimation parameters. Continuing with the example above,
//CREATE AND SPECIFY ALL STRUCTURES struct arimamtControl amc; struct arimamtOut amo; amc = arimamtControlCreate; //SPECIFIES TO INCLUDE CONSTANT (Default) amc.ar = 2; amc.diff = 0; amc.ma = 0; amc.const = 1;
This appropriately passes to arimamt that the model is an AR(2) model and includes a constants. The remaining controls elements, set at the default values, specify preferences for output display and estimation.
The Command File
Finally we put it all together in the command file:
//Aptech Systems, Inc //ARMA MODEL ANALYSIS //Objective: Determining model parameters library tsmt; #include tsmt.sdf //STEP ONE: RANDOMLY GENERATE DATA //SPECIFY ARMA PARAMETERS b_sim = { 0.5, -0.8 }; q_sim = 0; p_sim = 2; const_sim = 1.5; n_sim=200; k_sim=4; std_sim=1; seed_sim = 5012; y_sim = simarmamt(b_sim,p_sim,q_sim,const_sim,n_sim,k_sim,std_sim,seed_sim); //VISUALIZE DATA x_plot=seqa(1, 1, rows(y_sim)); //PLOT Y OBS plotXY(x_plot, y_sim); //CREATE AND SPECIFY ALL STRUCTURES struct arimamtControl amc; struct arimamtOut amo; amc = arimamtControlCreate; amc.ar = 2; amc.diff = 0; amc.ma = 0; amc.const = 1; //DEFINE OUTPUT FILE output file = arimamt1.out reset; //CALL ARIMA ROUTINE amo = arimamt(amc, y_sim[1:rows(y_sim),1]); //FORECAST USING RESULTS f = tsforecastmt(amc, amo.b, z_sim,amo.e,10); {alpha, tstat, zcrit} = vmadfmt(z_sim, -1, 3); output off;
This program produces the following output:
Model: ARIMA(2,0,0) Final Results: Iterations Until Convergence: 18 Log Likelihood: -286.868570 Number of Residuals: 200 AIC : 579.737139 Error Variance : 1.032660354 SBC : 589.632091 Standard Error : 1.016198974 DF: 197 Adj. SSE: 206.257625544 SSE: 203.434089775 Coefficients Std. Errors T-Ratio Approx. Prob. AR1 0.67379713 1.27040000 0.53038 0.59644 AR2 -0.85412220 0.24522685 -3.48299 0.00061 Mean 1.49350894 0.07172501 20.82271 0.00000 Constant: 1.76282603 Total Computation Time: 0.03 (seconds) AR Roots and Moduli: Real : 0.39444 0.39444 Imag.: 1.00758 -1.00758 Mod. : 1.08203 1.08203 Forecasts for ARIMA(2,0,0) Model. 95% Confidence Interval Computed. Period LCL Forecasts UCL Forecast Std. Err. 201 -0.710854 1.280860 3.272573 1.016199 202 0.684241 3.085889 5.487537 1.225353 203 0.217663 2.748079 5.278494 1.291052 204 -2.060362 0.978747 4.017856 1.550594 205 -2.997646 0.075108 3.147862 1.567761 206 -2.297213 0.977464 4.252141 1.670784 207 -1.113494 2.357287 5.828068 1.770839 208 -0.959802 2.516286 5.992373 1.773547 209 -2.204658 1.444881 5.094420 1.862044 210 -3.108973 0.587167 4.283308 1.885821
