The time series model is automatically selected by using the AIC. The TSUNIMAR call estimates the univariate autoregressive model and computes the AIC. You need to specify the maximum lag or order of the AR process with the MAXLAG= option or put the maximum lag as the sixth argument of the TSUNIMAR call. Here is an example:

   proc iml;
   y = { 2.430 2.506 2.767 2.940 3.169 3.450 3.594 3.774 3.695 3.411
         2.718 1.991 2.265 2.446 2.612 3.359 3.429 3.533 3.261 2.612
         2.179 1.653 1.832 2.328 2.737 3.014 3.328 3.404 2.981 2.557
         2.576 2.352 2.556 2.864 3.214 3.435 3.458 3.326 2.835 2.476
         2.373 2.389 2.742 3.210 3.520 3.828 3.628 2.837 2.406 2.675
         2.554 2.894 3.202 3.224 3.352 3.154 2.878 2.476 2.303 2.360
         2.671 2.867 3.310 3.449 3.646 3.400 2.590 1.863 1.581 1.690
         1.771 2.274 2.576 3.111 3.605 3.543 2.769 2.021 2.185 2.588
         2.880 3.115 3.540 3.845 3.800 3.579 3.264 2.538 2.582 2.907
         3.142 3.433 3.580 3.490 3.475 3.579 2.829 1.909 1.903 2.033
         2.360 2.601 3.054 3.386 3.553 3.468 3.187 2.723 2.686 2.821
         3.000 3.201 3.424 3.531 };
   call tsunimar(arcoef,ev,nar,aic) data=y opt={-1 1} 
                                    print=1 maxlag=20;

You can also invoke the TSUNIMAR subroutine as follows:

   call tsunimar(arcoef,ev,nar,aic,y,20,{-1 1},,1);

The optional arguments can be omitted. In this example, the argument MISSING is omitted, and thus the default option (MISSING=0) is used. The summary table of the minimum AIC method is displayed in Figure 13.4 and Figure 13.5. The final estimates are given in Figure 13.6.

The minimum AIC order is selected as 11. Then the coefficients are estimated as in Figure 13.6. Note that the first 20 observations are used as presample values.

You can estimate the AR(11) model directly by specifying OPT= and using the first 11 observations as presample values. The AR(11) estimates shown in Figure 13.7 are different from the minimum AIC estimates in Figure 13.6 because the samples are slightly different. Here is the code:

   call tsunimar(arcoef,ev,nar,aic,y,11,{-1 0},,1);

The minimum AIC procedure can also be applied to the vector autoregressive (VAR) model by using the TSMULMAR subroutine. See the section Multivariate Time Series Analysis for details. Three variables are used as input. The maximum lag is specified as 10. Here is the code:

   data one;
      input invest income consum @@;
   datalines;
      . . . data lines omitted . . .
      ;
   proc iml;
   use one;
   read all into y var{invest income consum};
   mdel  = 1;
   maice = 2;
   misw  = 0;
   opt   = mdel ||  maice || misw;
   maxlag = 10;
   miss   = 0;
   print  = 1;
   call tsmulmar(arcoef,ev,nar,aic,y,maxlag,opt,miss,print);

The VAR(3) model minimizes the AIC and was selected as an appropriate model (see Figure 13.8). However, AICs of the VAR(4) and VAR(5) models show little difference from VAR(3). You can also choose VAR(4) or VAR(5) as an appropriate model in the context of minimum AIC since this AIC difference is much less than 1.

The TSMULMAR subroutine estimates the instantaneous response model with diagonal error variance. See the section Multivariate Time Series Analysis for details on the instantaneous response model. Therefore, it is possible to select the minimum AIC model independently for each equation. The best model is selected by specifying MAXLAG=5, as in the following code:

   call tsmulmar(arcoef,ev,nar,aic) data=y maxlag=5
        opt={1 1 0} print=1;

You can print the intermediate results of the minimum AIC procedure by using the PRINT=2 option.

Note that the AIC value depends on the MAXLAG=lag option and the number of parameters estimated. The minimum AIC VAR estimation procedure (MAICE=2) uses the following AIC formula:

In this formula, is the order of the -variate VAR process, and intercept=1 if the intercept is specified; otherwise, intercept=0. When you use the MAICE=1 or MAICE=0 option, AIC is computed as the sum of AIC for each response equation. Therefore, there is an AIC difference of since the instantaneous response model contains the additional response variables as regressors.

The following code estimates the instantaneous response model. The results are shown in Figure 13.10.

   call tsmulmar(arcoef,ev,nar,aic) data=y 
                           maxlag=3 opt={1 0 0};
   print aic nar;
   print arcoef;

The following code estimates the VAR model. The results are shown in Figure 13.11.

   call tsmulmar(arcoef,ev,nar,aic) data=y maxlag=3
        opt={1 2 0};
   print aic nar;
   print arcoef;

The AIC computed from the instantaneous response model is greater than that obtained from the VAR model estimation by 6. There is a discrepancy between Figure 13.11 and Figure 13.8 because different observations are used for estimation.

The following example shows how to manage nonstationary data by using TIMSAC calls. In practice, time series are considered to be stationary when the expected values of first and second moments of the series do not change over time. This weak or covariance stationarity can be modeled by using the TSMLOCAR, TSMLOMAR, TSDECOMP, and TSTVCAR subroutines.

First, the locally stationary model is estimated. The whole series (1000 observations) is divided into three blocks of size 300 and one block of size 90, and the minimum AIC procedure is applied to each block of the data set. See the section Nonstationary Time Series for more details. Here is the code:

   data one;
      input y @@;
   datalines;
      . . . data lines omitted . . .
      ;

   proc iml;
   use one;
   read all var{y};

   mdel = -1;
   lspan = 300; /* local span of data */
   maice = 1;
   opt = mdel || lspan || maice;
   call tsmlocar(arcoef,ev,nar,aic,first,last)
                   data=y maxlag=10 opt=opt print=2;

Estimation results are displayed with the graphs of power spectrum , where is a rational spectral density function. See the section Spectral Analysis. The estimates for the first block and third block are shown in Figure 13.12 and Figure 13.15, respectively. As the first block and the second block do not have any sizable difference, the pooled model (AIC=45.892) is selected instead of the moving model (AIC=46.957) in Figure 13.13. However, you can notice a slight change in the shape of the spectrum of the third block of the data (observations 611 through 910). See Figure 13.14 and Figure 13.16 for comparison. The moving model is selected since the AIC (106.830) of the moving model is smaller than that of the pooled model (108.867).

Finally, the moving model is selected since there is a structural change in the last block of data (observations 911 through 1000). The final estimates are stored in variables ARCOEF, EV, NAR, AIC, FIRST, and LAST. The final estimates and spectrum are given in Figure 13.17 and Figure 13.18, respectively. The power spectrum of the final model (Figure 13.18) is significantly different from that of the first and second blocks (see Figure 13.14).

The multivariate analysis for locally stationary data is a straightforward extension of the univariate analysis. The bivariate locally stationary VAR models are estimated. The selected model is the VAR(7) process with some zero coefficients over the last block of data. There seems to be a structural difference between observations from 11 to 610 and those from 611 to 896. Here is the code:

   proc iml;
   rudder = {. . . data lines omitted . . .};
   yawing = {. . . data lines omitted . . .};

   y = rudder` || yawing`;
   c = {0.01795 0.02419};
   *-- calibration of data --*/
   y = y # (c @ j(n,1,1));
   mdel = -1;
   lspan = 300; /* local span of data */
   maice = 1;
   call tsmlomar(arcoef,ev,nar,aic,first,last) data=y maxlag=10
           opt=(mdel || lspan || maice) print=1;

The results of the analysis are shown in Figure 13.19.

Consider the time series decomposition

where and are trend and seasonal components, respectively, and is a stationary AR() process. The annual real GNP series is analyzed under second difference stochastic constraints on the trend component and the stationary AR(2) process.

The seasonal component is ignored if you specify SORDER=0. Therefore, the following state space model is estimated:

where

The parameters of this state space model are , , , and . Here is the code:

   proc iml;
   y = { 116.8 120.1 123.2 130.2 131.4 125.6 124.5 134.3
         135.2 151.8 146.4 139.0 127.8 147.0 165.9 165.5
         179.4 190.0 189.8 190.9 203.6 183.5 169.3 144.2
         141.5 154.3 169.5 193.0 203.2 192.9 209.4 227.2
         263.7 297.8 337.1 361.3 355.2 312.6 309.9 323.7
         324.1 355.3 383.4 395.1 412.8 406.0 438.0 446.1
         452.5 447.3 475.9 487.7 497.2 529.8 551.0 581.1
         617.8 658.1 675.2 706.6 724.7 };
   y = y`; /*-- convert to column vector --*/
   mdel  = 0;
   trade = 0;
   tvreg = 0;
   year  = 0;
   period= 0;
   log   = 0;
   maxit = 100;
   update = .; /* use default update method      */
   line   = .; /* use default line search method */
   sigmax = 0; /* no upper bound for variances   */
   back = 100;
   opt = mdel || trade || year || period || log || maxit ||
         update || line || sigmax || back;
   call tsdecomp(cmp,coef,aic) data=y order=2 sorder=0 nar=2
        npred=5 opt=opt icmp={1 3} print=1;
   y = y[52:61];
   cmp = cmp[52:66,];
   print y cmp;

The estimated parameters are printed when you specify the PRINT= option. In Figure 13.20, the estimated variances are printed under the title of TAU2(I), showing that and . The AR coefficient estimates are and . These estimates are also stored in the output matrix COEF.

The trend and stationary AR components are estimated by using the smoothing method, and out-of-sample forecasts are computed by using a Kalman filter prediction algorithm. The trend and AR components are stored in the matrix CMP since the ICMP={1 3} option is specified. The last 10 observations of the original series Y and the last 15 observations of two components are shown in Figure 13.21. Note that the first column of CMP is the trend component and the second column is the AR component. The last 5 observations of the CMP matrix are out-of-sample forecasts.

Consider the simple time series decomposition

The TSBAYSEA subroutine computes seasonally adjusted series by estimating the seasonal component. The seasonally adjusted series is computed as . The details of the adjustment procedure are given in the section Bayesian Seasonal Adjustment.

The monthly labor force series (1972–1978) are analyzed. You do not need to specify the options vector if you want to use the default options. However, you should change OPT[2] when the data frequency is not monthly (OPT[2]=12). The NPRED= option produces the multistep forecasts for the trend and seasonal components. The stochastic constraints are specified as ORDER=2 and SORDER=1.

In Figure 13.22, the first column shows the trend components; the second column shows the seasonal components; the third column shows the forecasts; the fourth column shows the seasonally adjusted series; the last column shows the value of ABIC. The last 12 rows are the forecasts. The figure is generated by using the following statements:

   proc iml;
   y = { 5447 5412 5215 4697 4344 5426   
         5173 4857 4658 4470 4268 4116    
         4675 4845 4512 4174 3799 4847    
         4550 4208 4165 3763 4056 4058    
         5008 5140 4755 4301 4144 5380    
         5260 4885 5202 5044 5685 6106    
         8180 8309 8359 7820 7623 8569    
         8209 7696 7522 7244 7231 7195    
         8174 8033 7525 6890 6304 7655    
         7577 7322 7026 6833 7095 7022    
         7848 8109 7556 6568 6151 7453    
         6941 6757 6437 6221 6346 5880 }; 
   y = y`;

   call tsbaysea(trend,season,series,adj,abic)
        data=y order=2 sorder=1 npred=12 print=2;
   print trend season series adj abic;

The estimated spectral density function of the irregular series is shown in Figure 13.23 and Figure 13.24.

The forecast values of multivariate time series are computed by using the TSPRED call. In the following example, the multistep-ahead forecasts are produced from the VARMA(2,1) estimates. Since the VARMA model is estimated by using the mean deleted series, you should specify the CONSTANT= option. You need to provide the original series instead of the mean deleted series to get the correct predictions. The forecast variance MSE and the impulse response function IMPULSE are also produced.

The VARMA() model is written

Then the COEF matrix is constructed by stacking matrices . Here is the code:

   proc iml;
   c = { 264 235 239 239 275 277 274 334 334 306
         308 309 295 271 277 221 223 227 215 223
         241 250 270 303 311 307 322 335 335 334
         309 262 228 191 188 215 215 249 291 296 };
   f = { 690 690 688 690 694 702 702 702 700 702
         702 694 708 702 702 708 700 700 702 694
         698 694 700 702 700 702 708 708 710 704
         704 700 700 694 702 694 710 710 710 708 };
   t = { 1152 1288 1288 1288 1368 1456 1656 1496 1744 1464
         1560 1376 1336 1336 1296 1296 1280 1264 1280 1272
         1344 1328 1352 1480 1472 1600 1512 1456 1368 1280
         1224 1112 1112 1048 1176 1064 1168 1280 1336 1248 };
   p = { 254.14 253.12 251.85 250.41 249.09 249.19 249.52 250.19
         248.74 248.41 249.95 250.64 250.87 250.94 250.96 251.33
         251.18 251.05 251.00 250.99 250.79 250.44 250.12 250.19
         249.77 250.27 250.74 250.90 252.21 253.68 254.47 254.80
         254.92 254.96 254.96 254.96 254.96 254.54 253.21 252.08 };

   y = c` || f` || t` || p`;
   ar = {  .82028   -.97167    .079386   -5.4382,
          -.39983    .94448    .027938   -1.7477,
          -.42278  -2.3314    1.4682    -70.996,
           .031038  -.019231  -.0004904   1.3677,
          -.029811   .89262   -.047579    4.7873,
           .31476    .0061959 -.012221    1.4921,
           .3813    2.7182    -.52993    67.711,
          -.020818   .01764    .00037981  -.38154 };
   ma = {  .083035 -1.0509     .055898   -3.9778,
          -.40452    .36876    .026369    -.81146,
           .062379 -2.6506     .80784   -76.952,
           .03273   -.031555  -.00019776  -.025205 };
   coef = ar // ma;
   ev = { 188.55   6.8082    42.385   .042942,
          6.8082   32.169    37.995  -.062341,
          42.385   37.995    5138.8  -.10757,
          .042942  -.062341  -.10757  .34313 };

   nar = 2; nma = 1;
   call tspred(forecast,impulse,mse,y,coef,nar,nma,ev,
                  5,nrow(y),-1);

The first 40 forecasts in Figure 13.25 are one-step predictions. The last observation is the five-step forecast values of variables C and F. You can construct the confidence interval for these forecasts by using the mean square error matrix, MSE. See the section Multivariate Time Series Analysis for more details about impulse response functions and the mean square error matrix.

The TSROOT call computes the polynomial roots of the AR and MA equations. When the AR() process is written

you can specify the following polynomial equation:

When all roots of the preceding equation are inside the unit circle, the AR() process is stationary. The MA() process is invertible if the following polynomial equation has all roots inside the unit circle:

where are the MA coefficients. For example, the best AR model is selected and estimated by the TSUNIMAR subroutine (see Figure 13.26). You can obtain the roots of the preceding equation by calling the TSROOT subroutine. Since the TSROOT subroutine can handle the complex AR or MA coefficients, note that you should add zero imaginary coefficients for the second column of the MATIN matrix for real coefficients. Here is the code:

   proc iml;
   y = { 2.430 2.506 2.767 2.940 3.169 3.450 3.594 3.774 3.695 3.411
         2.718 1.991 2.265 2.446 2.612 3.359 3.429 3.533 3.261 2.612
         2.179 1.653 1.832 2.328 2.737 3.014 3.328 3.404 2.981 2.557
         2.576 2.352 2.556 2.864 3.214 3.435 3.458 3.326 2.835 2.476
         2.373 2.389 2.742 3.210 3.520 3.828 3.628 2.837 2.406 2.675
         2.554 2.894 3.202 3.224 3.352 3.154 2.878 2.476 2.303 2.360
         2.671 2.867 3.310 3.449 3.646 3.400 2.590 1.863 1.581 1.690
         1.771 2.274 2.576 3.111 3.605 3.543 2.769 2.021 2.185 2.588
         2.880 3.115 3.540 3.845 3.800 3.579 3.264 2.538 2.582 2.907
         3.142 3.433 3.580 3.490 3.475 3.579 2.829 1.909 1.903 2.033
         2.360 2.601 3.054 3.386 3.553 3.468 3.187 2.723 2.686 2.821
         3.000 3.201 3.424 3.531 };

   call tsunimar(ar,v,nar,aic) data=y maxlag=5
        opt=({-1 1}) print=1;
   /*-- set up complex coefficient matrix --*/
   ar_cx = ar || j(nrow(ar),1,0);
   call tsroot(root) matin=ar_cx nar=nar nma=0 print=1;

In Figure 13.27, the roots and their lengths from the origin are shown. The roots are also stored in the matrix ROOT. All roots are within the unit circle, while the MOD values of the fourth and fifth roots appear to be sizable (0.9194).

The TSROOT subroutine can also recover the polynomial coefficients if the roots are given as an input. You should specify the QCOEF=1 option when you want to compute the polynomial coefficients instead of polynomial roots. You can compare the result with the preceding output of the TSUNIMAR call. Here is the code:

   call tsroot(ar_cx) matin=root nar=nar qcoef=1
                       nma=0 print=1;

The results are shown in Figure 13.28.