| 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=
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;
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:
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 10.10.
call tsmulmar(arcoef,ev,nar,aic) data=y maxlag=3
opt={1 0 0};
print aic nar;
print arcoef;
| |||||||||||||||||||||||||||||||||||||
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;
| |||||||||||||||||||||||||||||||||||||
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 
, 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 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
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.
![{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 ] )](images/timeseriesexpls_timeseriesexplseq24.gif)
, 
, 
, 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 10.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.
<<< 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.
|
Figure 10.21: Smoothed and Predicted Values of Two Components
Seasonal Adjustment
Consider the simple time series decomposition
.
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 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;
|
Figure 10.22: Trend and Seasonal Component Estimates and Forecasts
The estimated spectral density function of the irregular
series 
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(
) 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);
|
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() process is written
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:
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).
|
Figure 10.25: Minimum AIC AR Estimation
| |||||||||||||||||||||||||||||||||||||||||||
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.
| |||||||||||||||||||||
Figure 10.27: Polynomial Coefficients
Copyright © 2009 by SAS Institute Inc., Cary, NC, USA. All rights reserved.