Time Series Analysis and Examples

Getting Started

Minimum AIC Model Selection

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 10.4 and Figure 10.5. The final estimates are given in Figure 10.6.

 ORDER  INNOVATION VARIANCE
   M           V(M)              AIC(M)
   0          0.31607294   -108.26753229
   1          0.11481982   -201.45277331
   2          0.04847420   -280.51201122
   3          0.04828185   -278.88576251
   4          0.04656506   -280.28905616
   5          0.04615922   -279.11190502
   6          0.04511943   -279.25356641
   7          0.04312403   -281.50543541
   8          0.04201118   -281.96304075
   9          0.04128036   -281.61262868
  10          0.03829179   -286.67686828
  11          0.03318558   -298.13013264
  12          0.03255171   -297.94298716
  13          0.03247784   -296.15655602
  14          0.03237083   -294.46677874
  15          0.03234790   -292.53337704
  16          0.03187416   -291.92021487
  17          0.03183282   -290.04220196
  18          0.03126946   -289.72064823
  19          0.03087893   -288.90203735
  20          0.02998019   -289.67854830
 



Figure 10.4: Minimum AIC Table - I



                           AIC(M)-AICMIN (truncated at 40.0)
                       0        10        20        30        40
  M   AIC(M)-AICMIN    +---------+---------+---------+---------+
  0    189.862600      |                                       .
  1     96.677359      |                                       .
  2     17.618121      |                 *                     |
  3     19.244370      |                  *                    |
  4     17.841076      |                 *                     |
  5     19.018228      |                  *                    |
  6     18.876566      |                  *                    |
  7     16.624697      |                *                      |
  8     16.167092      |               *                       |
  9     16.517504      |                *                      |
 10     11.453264      |          *                            |
 11             0      *                                       |
 12      0.187145      *                                       |
 13      1.973577      | *                                     |
 14      3.663354      |   *                                   |
 15      5.596756      |     *                                 |
 16      6.209918      |     *                                 |
 17      8.087931      |       *                               |
 18      8.409484      |       *                               |
 19      9.228095      |        *                              |
 20      8.451584      |       *                               |
                       +---------+---------+---------+---------+
 
 



Figure 10.5: Minimum AIC Table - II

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

 ..........................M  A  I  C  E.........................
 .                                                              .
 .                                                              .
 .                                                              .
 .     M     AR Coefficients: AR(M)                             .
 .                                                              .
 .     1     1.181322                                           .
 .     2    -0.551571                                           .
 .     3     0.231372                                           .
 .     4    -0.178040                                           .
 .     5     0.019874                                           .
 .     6    -0.062573                                           .
 .     7     0.028569                                           .
 .     8    -0.050710                                           .
 .     9     0.199896                                           .
 .    10     0.161819                                           .
 .    11    -0.339086                                           .
 .                                                              .
 .                                                              .
 .      AIC =   -298.1301326                                    .
 .      Innovation Variance =    0.033186                       .
 .                                                              .
 .                                                              .
 .            INPUT DATA   START =    21  FINISH =   114        .
 ................................................................
 
 



Figure 10.6: Minimum AIC Estimation

You can estimate the AR(11) model directly by specifying OPT=\{-1 0\} and using the first 11 observations as presample values. The AR(11) estimates shown in Figure 10.7 are different from the minimum AIC estimates in Figure 10.6 because the samples are slightly different. Here is the code:

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


         ..........................M  A  I  C  E.........................
         .                                                              .
         .                                                              .
         .                                                              .
         .     M     AR Coefficients: AR(M)                             .
         .                                                              .
         .     1     1.149416                                           .
         .     2    -0.533719                                           .
         .     3     0.276312                                           .
         .     4    -0.326420                                           .
         .     5     0.169336                                           .
         .     6    -0.164108                                           .
         .     7     0.073123                                           .
         .     8    -0.030428                                           .
         .     9     0.151227                                           .
         .    10     0.192808                                           .
         .    11    -0.340200                                           .
         .                                                              .
         .                                                              .
         .      AIC =   -318.7984105                                    .
         .      Innovation Variance =    0.036563                       .
         .                                                              .
         .                                                              .
         .            INPUT DATA   START =    12  FINISH =   114        .
         ................................................................        
 
 



Figure 10.7: AR(11) Estimation

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;  /* instantaneous modeling ? */ 
       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 10.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.

 ORDER  INNOVATION VARIANCE
   M         LOG(|V(M)|)         AIC(M)
   0         25.98001095   2136.36089828
   1         15.70406486   1311.73331883
   2         15.48896746   1312.09533158
   3         15.18567834   1305.22562428
   4         14.96865183   1305.42944974
   5         14.74838535   1305.36759889
   6         14.60269347   1311.42086432
   7         14.54981887   1325.08514729
   8         14.38596333   1329.64899297
   9         14.16383772   1329.43469312
  10         13.85377849   1322.00983656
 
 
                           AIC(M)-AICMIN (truncated at 40.0)
                        0        10        20        30        40
   M   AIC(M)-AICMIN    +---------+---------+---------+---------+
   0    831.135274      |                                       .
   1      6.507695      |      *                                |
   2      6.869707      |      *                                |
   3             0      *                                       |
   4      0.203825      *                                       |
   5      0.141975      *                                       |
   6      6.195240      |     *                                 |
   7     19.859523      |                   *                   |
   8     24.423369      |                       *               |
   9     24.209069      |                       *               |
  10     16.784212      |                *                      |
                        +---------+---------+---------+---------+
 



Figure 10.8: VAR Model Selection

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;
 


 
VARIANCE
256.64375 29.803549 76.846777
29.803549 228.97341 119.60387
76.846777 119.60387 134.21764
 
AR_COEFFICIENT
13.312109 1.5459098 15.963897
0.8257397 0.2514803 0
0.0958916 1.0057088 0
0.0320985 0.3544346 0.4698934
0.044719 -0.201035 0
0.0051931 -0.023346 0
0.1169858 -0.060196 0.0483318
0.1867829 0 0
0.0216907 0 0
-0.117786 0 0.3500366
0.1541108 0 0
0.0178966 0 0
0.0461454 0 -0.191437
-0.389644 0 0
-0.045249 0 0
-0.116671 0 0
 
AIC
1347.6198



Figure 10.9: Model Selection via Instantaneous Response Model

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:

(t-lag) \log(|\hat{\sigma}|) +   2(p x n^2 + n x {intercept})
In this formula, p is the order of the n-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 n(n-1) since the instantaneous response model contains the additional n(n-1)/2 response variables as regressors.

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

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


 
AIC NAR
1403.0762 3
 
ARCOEF
4.8245814 5.3559216 17.066894
0.8855926 0.3401741 -0.014398
0.1684523 1.0502619 0.107064
0.0891034 0.4591573 0.4473672
-0.059195 -0.298777 0.1629818
0.1128625 -0.044039 -0.088186
0.1684932 -0.025847 -0.025671
0.0637227 -0.196504 0.0695746
-0.226559 0.0532467 -0.099808
-0.303697 -0.139022 0.2576405



Figure 10.10: AIC from Instantaneous Response Model

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

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


 
AIC NAR
1397.0762 3
 
ARCOEF
4.8245814 5.3559216 17.066894
0.8855926 0.3401741 -0.014398
0.1684523 1.0502619 0.107064
0.0891034 0.4591573 0.4473672
-0.059195 -0.298777 0.1629818
0.1128625 -0.044039 -0.088186
0.1684932 -0.025847 -0.025671
0.0637227 -0.196504 0.0695746
-0.226559 0.0532467 -0.099808
-0.303697 -0.139022 0.2576405



Figure 10.11: AIC from VAR Model

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 10.11 and Figure 10.8 because different observations are used for estimation.

Nonstationary Data Analysis

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 (\log 10(f_{yy}(g))), where f_{yy}(g) is a rational spectral density function. See the section "Spectral Analysis". The estimates for the first block and third block are shown in Figure 10.12 and Figure 10.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 10.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 10.14 and Figure 10.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).

 INITIAL LOCAL MODEL: N_CURR =   300
                    NAR_CURR =   8  
                         AIC =   37.583203
 
 
 
 ..........................CURRENT MODEL.........................
 .                                                              .
 .                                                              .
 .                                                              .
 .     M     AR Coefficients: AR(M)                             .
 .                                                              .
 .     1     1.605717                                           .
 .     2    -1.245350                                           .
 .     3     1.014847                                           .
 .     4    -0.931554                                           .
 .     5     0.394230                                           .
 .     6    -0.004344                                           .
 .     7     0.111608                                           .
 .     8    -0.124992                                           .
 .                                                              .
 .                                                              .
 .      AIC =     37.5832030                                    .
 .      Innovation Variance =    1.067455                       .
 .                                                              .
 .                                                              .
 .            INPUT DATA   START =    11  FINISH =   310        .
 ................................................................
 
 



Figure 10.12: Locally Stationary Model for First Block



    --- THE FOLLOWING TWO MODELS ARE COMPARED ---
 
      MOVING MODEL:    (N_PREV =  300, N_CURR = 300)
                      NAR_CURR =  7  
                           AIC =  46.957398
      CONSTANT MODEL: N_POOLED =  600
                    NAR_POOLED =  8  
                           AIC =  45.892350
 
      *****  CONSTANT MODEL ADOPTED  *****
 
 
 
 ..........................CURRENT MODEL.........................
 .                                                              .
 .                                                              .
 .                                                              .
 .     M     AR Coefficients: AR(M)                             .
 .                                                              .
 .     1     1.593890                                           .
 .     2    -1.262379                                           .
 .     3     1.013733                                           .
 .     4    -0.926052                                           .
 .     5     0.314480                                           .
 .     6     0.193973                                           .
 .     7    -0.058043                                           .
 .     8    -0.078508                                           .
 .                                                              .
 .                                                              .
 .      AIC =     45.8923501                                    .
 .      Innovation Variance =    1.047585                       .
 .                                                              .
 .                                                              .
 .            INPUT DATA   START =    11  FINISH =   610        .
 ................................................................
 
 



Figure 10.13: Locally Stationary Model Comparison

         POWER SPECTRAL DENSITY
 20.00+
      |
      |
      |
      |
      |  XXXX
      XXX    XX    XXX
      |        XXXX
      |               X
      |
 10.00+
      |                X
      |
      |                 X
      |
      |                  X               XX
      |                                    X
      |                   X             X
      |
      |                    X           X    X
     0+                     X
      |                      X        X      X
      |                       XX    XX
      |                         XXXX          X
      |
      |                                        X
      |                                         X
      |
      |                                          X
      |                                           X
 -10.0+                                            X
      |                                             XX
      |                                               XX
      |                                                 XX
      |                                                   XXX
      |                                                      XXXXXX
      |
      |
      |
      |
 -20.0+-----------+-----------+-----------+-----------+-----------+
     0.0         0.1         0.2         0.3         0.4         0.5
 
                                FREQUENCY
 



Figure 10.14: Power Spectrum for First and Second Blocks



 --- THE FOLLOWING TWO MODELS ARE COMPARED ---
 
      MOVING MODEL: (N_PREV =  600, N_CURR = 300)
                   NAR_CURR =    7  
                        AIC =  106.829869
      CONSTANT MODEL: N_POOLED =  900
                    NAR_POOLED =    8 
                           AIC =  108.867091
 
 *************************************
 *****                           *****
 *****     NEW MODEL ADOPTED     *****
 *****                           *****
 *************************************
 
 
 ..........................CURRENT MODEL.........................
 .                                                              .
 .                                                              .
 .                                                              .
 .     M     AR Coefficients: AR(M)                             .
 .                                                              .
 .     1     1.648544                                           .
 .     2    -1.201812                                           .
 .     3     0.674933                                           .
 .     4    -0.567576                                           .
 .     5    -0.018924                                           .
 .     6     0.516627                                           .
 .     7    -0.283410                                           .
 .                                                              .
 .                                                              .
 .      AIC =     60.9375188                                    .
 .      Innovation Variance =    1.161592                       .
 .                                                              .
 .                                                              .
 .            INPUT DATA   START =   611  FINISH =   910        .
 ................................................................
 



Figure 10.15: Locally Stationary Model for Third Block



  POWER SPECTRAL DENSITY
 20.00+             X
      |            X
      |           X
      |              X
      |        XXX
      |   XXXXX
      | XX
      XX              X
      |
      |
 10.00+                X
      |
      |
      |                 X
      |
      |                  X
      |                                   X
      |                   X              X
      |                                    X
      |                    X            X
     0+                     X          X    X
      |                      X
      |                       X      XX      X
      |                        XXXXXX
      |                                       X
      |
      |                                        X
      |
      |                                         X
      |                                          X
 -10.0+                                           X
      |                                            XX
      |                                              XX       XXXXX
      |                                                XXXXXXX
      |
      |
      |
      |
      |
      |
 -20.0+-----------+-----------+-----------+-----------+-----------+
     0.0         0.1         0.2         0.3         0.4         0.5
 
                                FREQUENCY
 



Figure 10.16: Power Spectrum for Third Block

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 10.17 and Figure 10.18, respectively. The power spectrum of the final model (Figure 10.18) is significantly different from that of the first and second blocks (see Figure 10.14).

   --- THE FOLLOWING TWO MODELS ARE COMPARED ---
 
        MOVING MODEL:    (N_PREV =  300, N_CURR =  90)
                        NAR_CURR =    6  
                             AIC =  139.579012
        CONSTANT MODEL: N_POOLED = 390
                      NAR_POOLED =   9  
                             AIC =  167.783711
 
 *************************************
 *****                           *****
 *****     NEW MODEL ADOPTED     *****
 *****                           *****
 *************************************
 
 
 
 ..........................CURRENT MODEL.........................
 .                                                              .
 .                                                              .
 .                                                              .
 .     M     AR Coefficients: AR(M)                             .
 .                                                              .
 .     1     1.181022                                           .
 .     2    -0.321178                                           .
 .     3    -0.113001                                           .
 .     4    -0.137846                                           .
 .     5    -0.141799                                           .
 .     6     0.260728                                           .
 .                                                              .
 .                                                              .
 .      AIC =     78.6414932                                    .
 .      Innovation Variance =    2.050818                       .
 .                                                              .
 .                                                              .
 .            INPUT DATA   START =   911  FINISH =  1000        .
 ................................................................
 
 



Figure 10.17: Locally Stationary Model for Last Block



    POWER SPECTRAL DENSITY
 30.00+
      |
      |
      |
      |
      |           X
      |
      |            X
      |
      |
 20.00+          X
      |
      |
      |         X   X
      |
      |        X
      XXX     X
      |  XXXXX       X
      |
      |
 10.00+               X
      |
      |                X
      |
      |                 X
      |
      |                  X
      |                   X
      |                    X
      |                     XX
     0+                       XX      XXXXXX
      |                         XXXXXX      XX
      |                                       XX
      |                                         XX               XX
      |                                           XX         XXXX
      |                                             XXXXXXXXX
      |
      |
      |
      |
 -10.0+-----------+-----------+-----------+-----------+-----------+
     0.0         0.1         0.2         0.3         0.4         0.5
 
                                FREQUENCY
 



Figure 10.18: Power Spectrum for Last Block

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 10.19.

 --- THE FOLLOWING TWO MODELS ARE COMPARED ---
 
 MOVING MODEL:    (N_PREV =  600, N_CURR = 286)
                 NAR_CURR =    7  
                      AIC = -823.845234
 CONSTANT MODEL: N_POOLED =  886
               NAR_POOLED =   10  
                      AIC = -716.818588
 
 *************************************
 *****                           *****
 *****     NEW MODEL ADOPTED     *****
 *****                           *****
 *************************************
 
 
 
 ..........................CURRENT MODEL.........................
 .                                                              .
 .                                                              .
 .                                                              .
 .     M     AR Coefficients                                    .
 .                                                              .
 .     1     0.932904    -0.130964                              .
 .          -0.024401     0.599483                              .
 .     2     0.163141     0.266876                              .
 .          -0.135605     0.377923                              .
 .     3    -0.322283     0.178194                              .
 .           0.188603    -0.081245                              .
 .     4     0.166094    -0.304755                              .
 .          -0.084626    -0.180638                              .
 .     5            0            0                              .
 .                  0    -0.036958                              .
 .     6            0            0                              .
 .                  0     0.034578                              .
 .     7            0            0                              .
 .                  0     0.268414                              .
 .                                                              .
 .                                                              .
 .      AIC =   -114.6911872                                    .
 .                                                              .
 .           Innovation Variance                                .
 .                                                              .
 .           1.069929     0.145558                              .
 .           0.145558     0.563985                              .
 .                                                              .
 .                                                              .
 .            INPUT DATA   START =   611  FINISH =   896        .
 ................................................................
 



Figure 10.19: Locally Stationary VAR Model Analysis

Consider the time series decomposition

y_t = t_t + s_t + u_t + \epsilon_t
where t_t and s_t are trend and seasonal components, respectively, and u_t is a stationary AR(p) process. The annual real GNP series is analyzed under second difference stochastic constraints on the trend component and the stationary AR(2) process.
t_t & = & 2t_{t-1} - t_{t-2} + w_{1t} \    u_t & = & \alpha_1 u_{t-1} + \alpha_2 u_{t-2} + w_{2t}
The seasonal component is ignored if you specify SORDER=0. Therefore, the following state space model is estimated:
y_t & = & {h}{z}_t + \epsilon_t \    {z}_t & = & {f}{z}_{t-1} + {w}_t
where
{h}& = & [ 1 & 0 & 1 & 0    ] \    {f}& = & [ 2 & -1 & 0 & 0 \    1 & 0 & 0 & 0 \   ...   ... & 0 & 0 \    0 & 0 & 0 & 0 \    0 & 0 & \sigma_2^2 & 0 \    0 & 0 & 0 & 0    ]    )
The parameters of this state space model are \sigma^2_1, \sigma^2_2, \alpha_1, and \alpha_2. 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 10.20, the estimated variances are printed under the title of TAU2(I), showing that \hat{\sigma}_1^2=2.915 and \hat{\sigma}_2^2=113.9577. The AR coefficient estimates are \hat{\alpha}_1=1.397 and \hat{\alpha}_2=-0.595. These estimates are also stored in the output matrix COEF.

                      <<<  Final Estimates  >>>
 
                   ---  PARAMETER VECTOR  ---
 
    1.607426E-01 6.283837E+00 8.761628E-01 -5.94879E-01
 
                      ---  GRADIENT  ---
 
    3.385021E-04 5.760929E-06 3.029534E-04 -1.18396E-04
 
 
    LIKELIHOOD = -249.937193     SIG2 =       18.135052
     AIC        =  509.874385
 
 
    I   TAU2(I)       AR(I)      PARCOR(I)
    1     2.915075     1.397374     0.876163
    2   113.957712    -0.594879    -0.594879
 



Figure 10.20: Nonstationary Time Series and State Space Modeling

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 10.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.

 
Y CMP  
487.7 514.01142 -26.94343
497.2 532.62744 -32.48674
529.8 552.02403 -24.46593
551 571.90122 -20.15113
581.1 592.31944 -10.58647
617.8 613.21855 5.2504378
658.1 634.43665 20.799213
675.2 655.70431 22.161602
706.6 677.21249 27.927985
724.7 698.72363 25.95797
  720.23477 19.659209
  741.74591 12.029407
  763.25705 5.1147239
  784.76819 -0.008863
  806.27933 -3.055027



Figure 10.21: Smoothed and Predicted Values of Two Components

Seasonal Adjustment

Consider the simple time series decomposition

y_t = t_t + s_t + \epsilon_t
The TSBAYSEA subroutine computes seasonally adjusted series by estimating the seasonal component. The seasonally adjusted series is computed as y_t^* = y_t - \hat{s}_t. 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.

t_t & = & 2t_{t-1} - t_{t-2} + w_{1t} \   s_t & = & -s_{t-1} -  ...  - s_{t-11} + w_{2t}

In Figure 10.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;
 


 
OBS TREND SEASON SERIES ADJ ABIC
1 4843.2502 576.86675 5420.1169 4870.1332 874.04585
2 4848.6664 612.79607 5461.4624 4799.2039  
3 4871.2876 324.02004 5195.3077 4890.98  
4 4896.6633 -198.7601 4697.9032 4895.7601  
5 4922.9458 -572.5562 4350.3896 4916.5562  
. . . . .  
71 6551.6017 -266.2162 6285.3855 6612.2162  
72 6388.9012 -440.3472 5948.5539 6320.3472  
73 6226.2006 650.7707 6876.9713    
74 6063.5001 800.93733 6864.4374    
75 5900.7995 396.19866 6296.9982    
76 5738.099 -340.2852 5397.8137    
77 5575.3984 -719.1146 4856.2838    
78 5412.6979 553.19764 5965.8955    
79 5249.9973 202.06582 5452.0631    
80 5087.2968 -54.44768 5032.8491    
81 4924.5962 -295.2747 4629.3215    
82 4761.8957 -487.6621 4274.2336    
83 4599.1951 -266.1917 4333.0034    
84 4436.4946 -440.3354 3996.1591    



Figure 10.22: Trend and Seasonal Component Estimates and Forecasts

The estimated spectral density function of the irregular series \hat{\epsilon}_t is shown in Figure 10.23.

     I   Rational   0.0       10.0      20.0      30.0      40.0      50.0      60.0
           Spectrum  +---------+---------+---------+---------+---------+---------+
     0 1.366798E+00  |*                                                           ===>X
     1 1.571261E+00  |*
     2 2.414836E+00  |  *
     3 5.151906E+00  |      *
     4 1.634887E+01  |           *
     5 8.085674E+01  |                  *
     6 3.805530E+02  |                        *
     7 8.082536E+02  |                            *
     8 6.366350E+02  |                           *
     9 3.479435E+02  |                        *
    10 3.872650E+02  |                        *                                   ===>X
    11 1.264805E+03  |                              *
    12 1.726138E+04  |                                         *
    13 1.559041E+03  |                              *
    14 1.276516E+03  |                              *
    15 3.861089E+03  |                                  *
    16 9.593184E+03  |                                      *
    17 3.662145E+03  |                                  *
    18 5.499783E+03  |                                    *
    19 4.443303E+03  |                                   *
    20 1.238135E+03  |                             *                              ===>X
    21 8.392131E+02  |                            *
    22 1.258933E+03  |                              *
    23 2.932003E+03  |                                 *
    24 1.857923E+03  |                               *
    25 1.171437E+03  |                             *
    26 1.611958E+03  |                               *
    27 4.822498E+03  |                                   *
    28 4.464961E+03  |                                   *
    29 1.951547E+03  |                               *
    30 1.653182E+03  |                               *                            ===>X
    31 2.308152E+03  |                                *
    32 5.475758E+03  |                                    *
    33 2.349584E+04  |                                          *
    34 5.266969E+03  |                                    *
    35 2.058667E+03  |                                *
    36 2.215595E+03  |                                *
    37 8.181540E+03  |                                      *
    38 3.077329E+03  |                                 *
    39 7.577961E+02  |                           *
    40 5.057636E+02  |                          *                                 ===>X
    41 7.312090E+02  |                           *
    42 3.131377E+03  |                                 *                          ===>T
    43 8.173276E+03  |                                      *
    44 1.958359E+03  |                               *
    45 2.216458E+03  |                                *
    46 4.215465E+03  |                                   *
    47 9.659340E+02  |                            *
    48 3.758466E+02  |                        *
    49 2.849326E+02  |                       *
    50 3.617848E+02  |                        *                                   ===>X
    51 7.659839E+02  |                           *
    52 3.191969E+03  |                                  *
    53 1.768107E+04  |                                         *
    54 5.281385E+03  |                                    *
    55 2.959704E+03  |                                 *
    56 3.783522E+03  |                                  *
    57 1.896625E+04  |                                         *
    58 1.041753E+04  |                                       *
    59 2.038940E+03  |                                *
    60 1.347568E+03  |                              *                             ===>X
 
 X: If peaks (troughs) appear at these frequencies, try lower (higher) values
    of rigid and watch ABIC
 T: If a peaks appears here try trading-day adjustment
 



Figure 10.23: Spectrum of Irregular Component

Miscellaneous Time Series Analysis Tools

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=-1 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(p,q) model is written

{y}_t + \sum_{i=1}^p {a}_i{y}_{t-i} =   \epsilon_t + \sum_{i=1}^q {m}_i\epsilon_{t-i}

Then the COEF matrix is constructed by stacking matrices {a}_1, ... ,{a}_p, {m}_1, ... ,{m}_q. 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);
 


 
OBSERVED
Y1
Y2 PREDICTED
P1
P2
264 690 269.950 700.750
235 690 256.764 691.925
239 688 239.996 693.467
239 690 242.320 690.951
275 694 247.169 693.214
277 702 279.024 696.157
274 702 284.041 700.449
334 702 286.890 701.580
334 700 321.798 699.851
306 702 330.355 702.383
308 702 297.239 700.421
309 694 302.651 701.928
295 708 294.570 696.261
271 702 283.254 703.936
277 702 269.600 703.110
221 708 270.349 701.557
223 700 231.523 705.438
227 700 233.856 701.785
215 702 234.883 700.185
223 694 229.156 701.837
241 698 235.054 697.060
250 694 249.288 698.181
270 700 257.644 696.665
303 702 272.549 699.281
311 700 301.947 701.667
307 702 306.422 700.708
322 708 304.120 701.204
335 708 311.590 704.654
335 710 320.570 706.389
334 704 315.127 706.439
309 704 311.083 703.735
262 700 288.159 702.801
228 700 251.352 700.805
191 694 226.749 700.247
188 702 199.775 696.570
215 694 202.305 700.242
215 710 222.951 696.451
249 710 226.553 704.483
291 710 259.927 707.610
296 708 291.446 707.861
    293.899 707.430
    293.477 706.933
    292.564 706.190
    290.313 705.384
    286.559 704.618



Figure 10.24: Multivariate ARMA Prediction

The first 40 forecasts in Figure 10.24 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(p) process is written

y_t = \sum_{i=1}^p \alpha_i y_{t-i} + \epsilon_t
you can specify the following polynomial equation:
z^p - \sum_{i=1}^p \alpha_i z^{p-i} = 0

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

z^q + \sum_{i=1}^q \theta_i z^{q-i} = 0
where \theta_i are the MA coefficients. For example, the best AR model is selected and estimated by the TSUNIMAR subroutine (see Figure 10.25). 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 10.26, 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).

 
LAG AR_COEF
1 1.3003068
2 -0.72328
3 0.2421928
4 -0.378757
5 0.1377273
 
AIC INNOVATION_VARINACE
-318.6138 0.0490554



Figure 10.25: Minimum AIC AR Estimation

 
Roots of AR Characteristic Polynomial
I Real Imaginary MOD(Z) ATAN(I/R) Degree
1 -0.29755 0.55991 0.6341 2.0593 117.9869
2 -0.29755 -0.55991 0.6341 -2.0593 -117.9869
3 0.40529 0 0.4053 0 0
4 0.74505 0.53866 0.9194 0.6260 35.8660
5 0.74505 -0.53866 0.9194 -0.6260 -35.8660

Z**5-AR(1)*Z**4-AR(2)*Z**3-AR(3)*Z**2-AR(4)*Z**1-AR(5)=0



Figure 10.26: Roots of AR Characteristic Polynomial Equation

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 10.27.

 
Polynomial Coefficents
I AR(real) AR(imag)
1 1.30031 0
2 -0.72328 5.55112E-17
3 0.24219 1.61885E-16
4 -0.37876 0
5 0.13773 -4.1674E-18



Figure 10.27: Polynomial Coefficients

Previous Page | Next Page | Top of Page