• Previous Page
• |
• Next Page

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

Seasonal Adjustment

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.

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= 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;

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.

• Previous Page
• |
• Next Page
• |
• Top of Page