### Selected R commands and output for Chapter 6 ### Ouput edited to simplify the file. > require(MTS) > da=read.table("m-hsoldhst6312.txt",header=T) > head(da) year mon hsold hstart 1 1963 1 42 79.0 2 1963 2 35 89.6 3 1963 3 44 124.8 4 1963 4 52 164.2 5 1963 5 58 172.7 6 1963 6 48 154.2 > dim(da) [1] 595 4 > zt=log(da[,3:4]) > tdx=da[,1]+da[,2]/12 > colnames(zt) <- c("Hsold","Hstart") > MTSplot(zt,tdx) ### See Figure 6.1 of the text > dzt=diffM(zt) ### Regular difference > ddzt=diffM(dzt,12) ### Seasonal difference > ccm(ddzt,lags=24) ## See Figure 6.2 of the text [1] "Covariance matrix:" Hsold Hstart Hsold 0.01234 0.00375 Hstart 0.00375 0.01480 CCM at lag: 0 [,1] [,2] [1,] 1.000 0.277 [2,] 0.277 1.000 Simplified matrix: CCM at lag: 1 - . + - CCM at lag: 2 - . . . CCM at lag: 3 . . . . CCM at lag: 4 . . . . CCM at lag: 5 . . . . CCM at lag: 6 . . . . CCM at lag: 7 . . . . CCM at lag: 8 + . - . CCM at lag: 9 . . + + CCM at lag: 10 . . . - CCM at lag: 11 . . . + CCM at lag: 12 - - - - CCM at lag: 13 . + . + CCM at lag: 14 . . . . CCM at lag: 15 . . . . CCM at lag: 16 . . . . CCM at lag: 17 . . . . CCM at lag: 18 . . . . CCM at lag: 19 . . . . CCM at lag: 20 - . . . CCM at lag: 21 . . . - CCM at lag: 22 . . . + CCM at lag: 23 . . . . CCM at lag: 24 . . . . Hit Enter for p-value plot of individual ccm: ## Plot not shown > #### Estimation of Seasonal Time Series models > m3=sVARMA(zt,order=c(0,1,3),sorder=c(0,1,1),s=12) Number of parameters: 16 initial estimates: 0.1985901 -0.08082791 0.05693362 0.02328988 0.01039118 -0.06812058 -0.3552086 0.3327042 -0.2644806 -0.03664861 -0.1818856 -0.1061054 0.8981288 0.04304338 0.02171152 0.8611352 Coefficient(s): Estimate Std. Error t value Pr(>|t|) [1,] 0.24587 0.04621 5.321 1.03e-07 *** [2,] -0.16847 0.04483 -3.758 0.000171 *** [3,] 0.02834 0.05220 0.543 0.587194 [4,] 0.03384 0.04359 0.776 0.437562 [5,] 0.06276 0.04369 1.437 0.150853 [6,] -0.06707 0.04147 -1.617 0.105797 [7,] -0.28993 0.04721 -6.141 8.20e-10 *** [8,] 0.41585 0.04627 8.988 < 2e-16 *** [9,] -0.16783 0.05099 -3.292 0.000996 *** [10,] 0.04792 0.04739 1.011 0.311883 [11,] -0.08497 0.04927 -1.725 0.084564 . [12,] -0.07360 0.04799 -1.534 0.125097 [13,] 0.84990 0.01963 43.296 < 2e-16 *** [14,] 0.03323 0.02406 1.381 0.167355 [15,] 0.04428 0.02444 1.812 0.070012 . [16,] 0.82424 0.02860 28.819 < 2e-16 *** --- Estimates in matrix form: Regular MA coefficient matrix MA( 1 )-matrix [,1] [,2] [1,] 0.246 -0.168 [2,] -0.290 0.416 MA( 2 )-matrix [,1] [,2] [1,] 0.0283 0.0338 [2,] -0.1678 0.0479 MA( 3 )-matrix [,1] [,2] [1,] 0.0628 -0.0671 [2,] -0.0850 -0.0736 Seasonal MA coefficient matrix MA( 12 )-matrix [,1] [,2] [1,] 0.8498980 0.03322685 [2,] 0.0442752 0.82423827 Residuals cov-matrix: resi resi resi 0.006941949 0.002633186 resi 0.002633186 0.007116746 ---- aic= -10.0117 bic= -9.8917 > m4=refsVARMA(m3,thres=1.2) Number of parameters: 13 initial estimates: 0.1985901 -0.08082791 0.01039118 -0.06812058 -0.3552086 0.3327042 -0.2644806 -0.1818856 -0.1061054 0.8981288 0.04304338 0.02171152 0.8611352 Coefficient(s): Estimate Std. Error t value Pr(>|t|) [1,] 0.25255 0.04521 5.586 2.32e-08 *** [2,] -0.16847 0.04357 -3.866 0.00011 *** [3,] 0.07582 0.03910 1.939 0.05247 . [4,] -0.06798 0.03953 -1.720 0.08547 . [5,] -0.29222 0.04659 -6.272 3.56e-10 *** [6,] 0.41585 0.04668 8.909 < 2e-16 *** [7,] -0.16783 0.04498 -3.732 0.00019 *** [8,] -0.08497 0.04840 -1.756 0.07912 . [9,] -0.07345 0.04619 -1.590 0.11182 [10,] 0.84990 0.01959 43.378 < 2e-16 *** [11,] 0.03296 0.02395 1.376 0.16878 [12,] 0.04387 0.02441 1.797 0.07231 . [13,] 0.82472 0.02846 28.980 < 2e-16 *** --- Estimates in matrix form: Regular MA coefficient matrix MA( 1 )-matrix [,1] [,2] [1,] 0.253 -0.168 [2,] -0.292 0.416 MA( 2 )-matrix [,1] [,2] [1,] 0.000 0 [2,] -0.168 0 MA( 3 )-matrix [,1] [,2] [1,] 0.0758 -0.0680 [2,] -0.0850 -0.0734 Seasonal MA coefficient matrix MA( 12 )-matrix [,1] [,2] [1,] 0.84989804 0.03296309 [2,] 0.04387419 0.82471809 Residuals cov-matrix: resi resi resi 0.006964991 0.002655917 resi 0.002655917 0.007139871 ---- aic= -10.0172 bic= -9.9197 #### Switching the MA factors > m5=sVARMA(zt,order=c(0,1,3),sorder=c(0,1,1),s=12,switch=T) Number of parameters: 16 initial estimates: 0.1985901 -0.08082791 0.05693362 0.02328988 0.01039118 -0.06812058 -0.3552086 0.3327042 -0.2644806 -0.03664861 -0.1818856 -0.1061054 0.8981288 0.04304338 0.02171152 0.8611352 Coefficient(s): Estimate Std. Error t value Pr(>|t|) [1,] 0.24578 0.04619 5.321 1.03e-07 *** [2,] -0.16847 0.04489 -3.753 0.000175 *** [3,] 0.02170 0.05252 0.413 0.679443 [4,] 0.03967 0.04409 0.900 0.368237 [5,] 0.05653 0.04369 1.294 0.195703 [6,] -0.06148 0.04175 -1.473 0.140883 [7,] -0.29686 0.04682 -6.340 2.30e-10 *** [8,] 0.41585 0.04593 9.053 < 2e-16 *** [9,] -0.16783 0.05119 -3.278 0.001044 ** [10,] 0.04792 0.04777 1.003 0.315751 [11,] -0.08497 0.04923 -1.726 0.084353 . [12,] -0.06748 0.04781 -1.411 0.158148 [13,] 0.84990 0.01724 49.296 < 2e-16 *** [14,] 0.03715 0.02088 1.779 0.075261 . [15,] 0.02910 0.02153 1.351 0.176536 [16,] 0.84752 0.02453 34.549 < 2e-16 *** --- Estimates in matrix form: Regular MA coefficient matrix MA( 1 )-matrix [,1] [,2] [1,] 0.246 -0.168 [2,] -0.297 0.416 MA( 2 )-matrix [,1] [,2] [1,] 0.0217 0.0397 [2,] -0.1678 0.0479 MA( 3 )-matrix [,1] [,2] [1,] 0.0565 -0.0615 [2,] -0.0850 -0.0675 Seasonal MA coefficient matrix MA( 12 )-matrix [,1] [,2] [1,] 0.84989804 0.03714659 [2,] 0.02909768 0.84752096 Residuals cov-matrix: resi resi resi 0.006938460 0.002605709 resi 0.002605709 0.007118109 ---- aic= -10.0087 bic= -9.8887 > ##### Principal component analysis > da=read.table("m-amdur.txt",header=T) > dur=da[,3:6]/1000 > V0=princomp(dur) > summary(V0) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Standard deviation 196.9977058 30.69968026 12.566213076 3.9317495305 Proportion of Variance 0.9720509 0.02360665 0.003955264 0.0003872027 Cumulative Proportion 0.9720509 0.99565753 0.999612797 1.0000000000 > M0=matrix(V0$loadings[,1:4],4,4) ## eigenvectors > m1=VAR(dur,1) ## VAR(1) fit Constant term: Estimates: 0.008423112 -0.1286131 -8.347717 2.803765 Std.Error: 5.770404 1.382151 4.351705 3.024527 AR coefficient matrix AR( 1 )-matrix [,1] [,2] [,3] [,4] [1,] 0.686 -0.0268 -0.001276 0.357 [2,] 0.116 0.9949 -0.000132 -0.102 [3,] 0.562 -0.0229 0.994850 -0.441 [4,] 0.108 0.0231 -0.002064 0.852 standard error [,1] [,2] [,3] [,4] [1,] 0.0717 0.03719 0.00504 0.0914 [2,] 0.0172 0.00891 0.00121 0.0219 [3,] 0.0541 0.02805 0.00380 0.0689 [4,] 0.0376 0.01949 0.00264 0.0479 Residuals cov-mtx: [,1] [,2] [,3] [,4] [1,] 49.717324 3.2290051 32.681085 17.7787583 [2,] 3.229005 2.8523717 3.226700 0.2811015 [3,] 32.681085 3.2267003 28.275746 4.8528690 [4,] 17.778758 0.2811015 4.852869 13.6587253 det(SSE) = 335.831 AIC = 5.946689 BIC = 6.174678 HQ = 6.03849 > V1=princomp(m1$residuals) > summary(V1) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Standard deviation 8.8492175 3.6873635 1.57201591 0.357258542 Proportion of Variance 0.8286264 0.1438735 0.02614947 0.001350561 Cumulative Proportion 0.8286264 0.9725000 0.99864944 1.000000000 > M1=matrix(V1$loadings[,1:4],4,4) > h4=matrix(c(1,0,-1,-1),4,1) > t(h4)%*%m1$Phi [,1] [,2] [,3] [,4] [1,] 0.01530744 -0.02704302 -0.9940621 -0.05383388 > t(h4)%*%m1$Ph0 [,1] [1,] 5.552374 > m2=VAR(dur,2) ### VAR(2) fit Constant term: Estimates: -0.2213344 3.248112 -6.266722 3.838599 Std.Error: 6.896798 1.278179 5.183505 3.678531 AR coefficient matrix AR( 1 )-matrix [,1] [,2] [,3] [,4] [1,] 1.033 1.012 -0.6384 -0.108 [2,] -0.445 1.549 0.4408 0.537 [3,] 0.307 0.645 1.0054 -0.120 [4,] 0.141 0.452 -0.0723 0.619 standard error [,1] [,2] [,3] [,4] [1,] 0.708 0.2837 0.717 0.682 [2,] 0.131 0.0526 0.133 0.126 [3,] 0.532 0.2132 0.539 0.513 [4,] 0.378 0.1513 0.383 0.364 AR( 2 )-matrix [,1] [,2] [,3] [,4] [1,] 0.2433 -1.028 0.6343 -0.115 [2,] 0.0642 -0.568 -0.4375 -0.166 [3,] 0.2473 -0.663 -0.0100 -0.336 [4,] -0.0161 -0.440 0.0704 0.227 standard error [,1] [,2] [,3] [,4] [1,] 0.0860 0.2814 0.713 0.1473 [2,] 0.0159 0.0522 0.132 0.0273 [3,] 0.0647 0.2115 0.536 0.1107 [4,] 0.0459 0.1501 0.380 0.0785 Residuals cov-mtx: [,1] [,2] [,3] [,4] [1,] 44.157838 1.2597770 28.800155 16.1262464 [2,] 1.259777 1.5166867 1.416396 -0.1285939 [3,] 28.800155 1.4163962 24.943639 4.2286227 [4,] 16.126246 -0.1285939 4.228623 12.5620846 det(SSE) = 97.43829 AIC = 4.839382 BIC = 5.29536 HQ = 5.022983 > V2=princomp(m2$residuals) > summary(V2) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Standard deviation 8.3227499 3.5233108 1.19099873 0.2826413635 Proportion of Variance 0.8327478 0.1492388 0.01705306 0.0009603979 Cumulative Proportion 0.8327478 0.9819865 0.99903960 1.0000000000 > M2=matrix(V2$loadings[,1:4],4,4) > print(round(M0,3)) [,1] [,2] [,3] [,4] [1,] -0.102 0.712 0.342 -0.604 [2,] -0.152 0.315 -0.928 -0.129 [3,] -0.978 -0.182 0.098 0.006 [4,] -0.096 0.600 0.110 0.786 > print(round(M1,3)) [,1] [,2] [,3] [,4] [1,] 0.794 0.161 0.066 0.583 [2,] 0.058 -0.109 -0.990 0.063 [3,] 0.547 -0.592 0.060 -0.588 [4,] 0.260 0.782 -0.106 -0.557 > print(round(M2,3)) [,1] [,2] [,3] [,4] [1,] 0.796 0.150 0.017 0.587 [2,] 0.026 -0.070 -0.997 0.012 [3,] 0.543 -0.600 0.049 -0.585 [4,] 0.267 0.783 -0.055 -0.560 > ### May continue with VAR(3) and VAR(4) fits. #### VARX modeling > da=read.table("m-gasoil.txt",header=T) > head(da) year mon Greg hoil oilp gasp 1 1993 11 1.066 0.502 16.699 2.32 2 1993 12 1.014 0.435 14.510 2.09 3 1994 1 0.998 0.499 15.000 2.41 4 1994 2 1.009 0.557 14.780 2.59 5 1994 3 1.008 0.492 14.660 2.08 6 1994 4 1.027 0.479 16.380 2.04 > yt=da[,3:6] > colnames(yt) <- c("Greg","hoil","oilp","gasp") > tdx=cda[,1]+da[,2]/12 Error: object 'cda' not found > tdx=da[,1]+da[,2]/12 > MTSplot(yt,tdx) Hit return for more plots: ### See Figure 6.6 of the text > zt=yt[,1:2]; xt=yt[,3:4] > VARorder(zt) ## VAR order selection for zt selected order: aic = 11 selected order: bic = 2 selected order: hq = 11 Summary table: p AIC BIC HQ M(p) p-value [1,] 0 -4.0481 -4.0481 -4.0481 0.0000 0.0000 [2,] 1 -8.6871 -8.6265 -8.6626 979.2772 0.0000 [3,] 2 -9.0847 -8.9637 -9.0359 89.8604 0.0000 [4,] 3 -9.1417 -8.9601 -9.0684 18.9893 0.0008 [5,] 4 -9.1252 -8.8830 -9.0274 3.8288 0.4297 [6,] 5 -9.1302 -8.8275 -9.0080 8.1411 0.0865 [7,] 6 -9.1840 -8.8207 -9.0374 17.8015 0.0013 [8,] 7 -9.1574 -8.7336 -8.9864 1.7373 0.7839 [9,] 8 -9.2060 -8.7216 -9.0105 16.4152 0.0025 [10,] 9 -9.2366 -8.6918 -9.0167 12.7826 0.0124 [11,] 10 -9.3124 -8.7070 -9.0681 21.2965 0.0003 [12,] 11 -9.3394 -8.6735 -9.0707 11.8239 0.0187 [13,] 12 -9.3205 -8.5940 -9.0273 3.0860 0.5435 [14,] 13 -9.3131 -8.5261 -8.9955 5.1917 0.2682 > VARXorder(zt,xt,maxp=11,maxm=2) ### VARX order specification selected order(p,s): aic = 11 1 selected order(p,s): bic = 2 1 selected order(p,s): hq = 2 1 > m1=VARX(zt,2,xt,1) ### Estimation constant term: est: 0.1817 -0.017 se: 0.0295 0.0186 AR( 1 ) matrix Greg hoil Greg 1.041 0.169 hoil 0.005 0.844 standard errors [,1] [,2] [1,] 0.059 0.08 [2,] 0.037 0.05 AR( 2 ) matrix Greg hoil Greg -0.322 -0.008 hoil 0.014 0.012 standard errors [,1] [,2] [1,] 0.056 0.068 [2,] 0.035 0.043 Coefficients of exogenous lag- 0 coefficient matrix oilp gasp Greg 0.018 0.010 hoil 0.024 0.026 standard errors [,1] [,2] [1,] 0.001 0.007 [2,] 0.001 0.004 lag- 1 coefficient matrix oilp gasp Greg -0.015 -0.012 hoil -0.019 -0.029 standard errors [,1] [,2] [1,] 0.002 0.007 [2,] 0.001 0.004 Residual Covariance Matrix Greg hoil Greg 0.00671 0.00104 hoil 0.00104 0.00267 =========== Information criteria: AIC: -10.83396 BIC: -10.55981 > cov(m1$residuals) Greg hoil Greg 0.006736072 0.001049648 hoil 0.001049648 0.002678639 > m1a=refVARX(m1,thres=1) ## refinement constant term: est: 0.1828 0 se: 0.028 1 AR( 1 ) matrix [,1] [,2] [1,] 1.044 0.162 [2,] 0.000 0.873 standard errors [,1] [,2] [1,] 0.054 0.056 [2,] 1.000 0.031 AR( 2 ) matrix [,1] [,2] [1,] -0.327 0 [2,] 0.000 0 standard errors [,1] [,2] [1,] 0.04 1 [2,] 1.00 1 Coefficients of exogenous lag- 0 coefficient matrix [,1] [,2] [1,] 0.018 0.010 [2,] 0.023 0.026 standard errors [,1] [,2] [1,] 0.001 0.007 [2,] 0.001 0.004 lag- 1 coefficient matrix [,1] [,2] [1,] -0.015 -0.012 [2,] -0.019 -0.029 standard errors [,1] [,2] [1,] 0.002 0.007 [2,] 0.001 0.004 Residual Covariance Matrix Greg hoil Greg 0.00671 0.00104 hoil 0.00104 0.00269 =========== Information criteria: AIC: -10.87015 BIC: -10.67216 > MTSdiag(m1a) ### See residual CCM of Figure 6.7 of the text [1] "Covariance matrix:" Greg hoil Greg 0.00674 0.00105 hoil 0.00105 0.00270 CCM at lag: 0 [,1] [,2] [1,] 1.000 0.246 [2,] 0.246 1.000 Simplified matrix: CCM at lag: 1 . . . . CCM at lag: 2 . . . . CCM at lag: 3 . . - . CCM at lag: 4 + . . . CCM at lag: 5 . . . . CCM at lag: 6 . . . + CCM at lag: 7 . + . . CCM at lag: 8 . - . . CCM at lag: 9 . . . . CCM at lag: 10 . . . . CCM at lag: 11 . . . . CCM at lag: 12 . . . . CCM at lag: 13 . . . . CCM at lag: 14 . . . . CCM at lag: 15 . . . . CCM at lag: 16 . . . . CCM at lag: 17 . . . . CCM at lag: 18 . . . . CCM at lag: 19 . . . . CCM at lag: 20 . . + . CCM at lag: 21 . . . . CCM at lag: 22 . . . . CCM at lag: 23 . . . . CCM at lag: 24 . . . . Hit Enter for p-value plot of individual ccm: ### Plot not shown Hit Enter to compute MQ-statistics: ### See Figure 6.8 of the text Ljung-Box Statistics: m Q(m) df p-value [1,] 1.00 2.62 4.00 0.62 [2,] 2.00 7.13 8.00 0.52 [3,] 3.00 15.41 12.00 0.22 [4,] 4.00 24.31 16.00 0.08 [5,] 5.00 25.92 20.00 0.17 [6,] 6.00 37.19 24.00 0.04 [7,] 7.00 47.79 28.00 0.01 [8,] 8.00 55.08 32.00 0.01 [9,] 9.00 59.08 36.00 0.01 [10,] 10.00 63.07 40.00 0.01 [11,] 11.00 68.59 44.00 0.01 [12,] 12.00 72.54 48.00 0.01 [13,] 13.00 75.84 52.00 0.02 [14,] 14.00 76.26 56.00 0.04 [15,] 15.00 80.16 60.00 0.04 [16,] 16.00 82.91 64.00 0.06 [17,] 17.00 85.56 68.00 0.07 [18,] 18.00 94.58 72.00 0.04 [19,] 19.00 97.49 76.00 0.05 [20,] 20.00 108.95 80.00 0.02 [21,] 21.00 112.01 84.00 0.02 [22,] 22.00 114.52 88.00 0.03 [23,] 23.00 115.25 92.00 0.05 [24,] 24.00 115.72 96.00 0.08 Hit Enter to obtain residual plots: ## PLot not shown > #### Multivariate regression models with time series errors > m1=Mlm(zt,xt) ### Multivariate Multiple Linear Regression [1] "LSE of parameters" [1] " est s.d. t-ratio prob" [,1] [,2] [,3] [,4] [1,] 0.62686 0.022768 27.53 0.0000 [2,] 0.02864 0.000424 67.59 0.0000 [3,] -0.00675 0.005222 -1.29 0.1977 [4,] -0.02659 0.016308 -1.63 0.1044 [5,] 0.02961 0.000303 97.56 0.0000 [6,] -0.00846 0.003740 -2.26 0.0246 > VARorder(m1$residuals) selected order: aic = 2 selected order: bic = 2 selected order: hq = 2 Summary table: p AIC BIC HQ M(p) p-value [1,] 0 -8.1889 -8.1889 -8.1889 0.0000 0.0000 [2,] 1 -10.0945 -10.0339 -10.0700 406.6203 0.0000 [3,] 2 -10.1920 -10.0710 -10.1432 27.5947 0.0000 [4,] 3 -10.1654 -9.9838 -10.0921 1.7968 0.7731 [5,] 4 -10.1373 -9.8951 -10.0395 1.4782 0.8305 [6,] 5 -10.1491 -9.8464 -10.0270 9.5232 0.0493 [7,] 6 -10.1383 -9.7751 -9.9917 4.9091 0.2967 [8,] 7 -10.1330 -9.7092 -9.9620 5.9411 0.2036 [9,] 8 -10.1488 -9.6644 -9.9533 10.0013 0.0404 [10,] 9 -10.1172 -9.5723 -9.8973 0.7397 0.9463 [11,] 10 -10.1044 -9.4990 -9.8601 4.3268 0.3636 [12,] 11 -10.0776 -9.4116 -9.8088 1.6273 0.8039 [13,] 12 -10.0699 -9.3434 -9.7767 5.1923 0.2681 [14,] 13 -10.0407 -9.2536 -9.7231 1.1487 0.8865 > m3=REGts(zt,2,xt) ### Estimation Number of parameters: 14 initial estimates: 0.6268644 0.02863685 -0.006746305 -0.02658854 0.02960791 -0.008464821 1.007292 -0.07880628 -0.3267944 0.187021 0.07138531 0.9385796 -0.05076304 -0.1429214 Par. Lower-bounds: 0.5972654 0.02808604 -0.01353439 -0.04778846 0.0292134 -0.01332671 0.918251 -0.2195028 -0.4161208 0.04592782 0.01328297 0.8467706 -0.1090514 -0.2349893 Par. Upper-bounds: 0.6564634 0.02918766 4.178094e-05 -0.00538862 0.03000243 -0.003602935 1.096334 0.06189025 -0.2374681 0.3281141 0.1294877 1.030389 0.007525317 -0.05085357 Final Estimates: 0.6495229 0.02808604 -0.004773827 -0.02477847 0.0292134 -0.003602935 1.011848 -0.07054513 -0.3260461 0.1864357 0.07024979 0.9491855 -0.04637316 -0.1383752 Coefficient(s): Estimate Std. Error t value Pr(>|t|) [1,] 0.6495229 0.0486092 13.362 < 2e-16 *** [2,] 0.0280860 0.0009738 28.841 < 2e-16 *** [3,] -0.0047738 0.0079791 -0.598 0.5496 [4,] -0.0247785 0.0431763 -0.574 0.5660 [5,] 0.0292134 0.0008800 33.195 < 2e-16 *** [6,] -0.0036029 0.0064218 -0.561 0.5748 [7,] 1.0118475 0.0686284 14.744 < 2e-16 *** [8,] -0.0705451 0.1132366 -0.623 0.5333 [9,] -0.3260461 0.0679005 -4.802 1.57e-06 *** [10,] 0.1864357 0.1096516 1.700 0.0891 . [11,] 0.0702498 0.0437465 1.606 0.1083 [12,] 0.9491855 0.0712961 13.313 < 2e-16 *** [13,] -0.0463732 0.0442779 -1.047 0.2950 [14,] -0.1383752 0.0719582 -1.923 0.0545 . --- ======= Coefficient matrix for constant + exogenous variable Estimates: [,1] [,2] [,3] [1,] 0.650 0.028 -0.005 [2,] -0.025 0.029 -0.004 Standard errors: [,1] [,2] [,3] [1,] 0.049 0.001 0.008 [2,] 0.043 0.001 0.006 VAR coefficient matrices: AR( 1 ) coefficient: [,1] [,2] [1,] 1.012 -0.071 [2,] 0.070 0.949 standard errors: [,1] [,2] [1,] 0.069 0.113 [2,] 0.044 0.071 AR( 2 ) coefficient: [,1] [,2] [1,] -0.326 0.186 [2,] -0.046 -0.138 standard errors: [,1] [,2] [1,] 0.068 0.110 [2,] 0.044 0.072 Residual Covariance matrix: Greg hoil Greg 0.009158 0.001943 hoil 0.001943 0.003754 ============ Information criteria: AIC: -10.26927 BIC: -10.05605 > m3c=refREGts(m3,thres=1.2) ### Refinement Number of parameters: 9 initial estimates: 0.6495229 0.02808604 0.0292134 1.011848 -0.3260461 0.1864357 0.07024979 0.9491855 -0.1383752 Par. Lower-bounds: 0.5766091 0.02662528 0.02789333 0.9089049 -0.4278969 0.0219583 0.004630105 0.8422414 -0.2463125 Par. Upper-bounds: 0.7224368 0.02954679 0.03053347 1.11479 -0.2241953 0.3509132 0.1358695 1.05613 -0.03043783 Final Estimates: 0.6787531 0.02662528 0.02789333 1.009309 -0.3054354 0.166147 0.04471049 0.9726905 -0.1276744 Coefficient(s): Estimate Std. Error t value Pr(>|t|) [1,] 0.6787531 0.0364469 18.623 < 2e-16 *** [2,] 0.0266253 0.0010140 26.258 < 2e-16 *** [3,] 0.0278933 0.0007738 36.046 < 2e-16 *** [4,] 1.0093088 0.0680000 14.843 < 2e-16 *** [5,] -0.3054354 0.0638083 -4.787 1.69e-06 *** [6,] 0.1661470 0.0703592 2.361 0.0182 * [7,] 0.0447105 0.0295577 1.513 0.1304 [8,] 0.9726905 0.0661847 14.697 < 2e-16 *** [9,] -0.1276744 0.0734336 -1.739 0.0821 . --- ======= Coefficient matrix for constant + exogenous variable Estimates: [,1] [,2] [,3] [1,] 0.679 0.027 0 [2,] 0.000 0.028 0 Standard errors: [,1] [,2] [,3] [1,] 0.036 0.001 1 [2,] 1.000 0.001 1 VAR coefficient matrices: AR( 1 ) coefficient: [,1] [,2] [1,] 1.009 0.000 [2,] 0.045 0.973 standard errors: [,1] [,2] [1,] 0.068 1.000 [2,] 0.030 0.066 AR( 2 ) coefficient: [,1] [,2] [1,] -0.305 0.166 [2,] 0.000 -0.128 standard errors: [,1] [,2] [1,] 0.064 0.070 [2,] 1.000 0.073 Residual Covariance matrix: Greg hoil Greg 0.008935 0.001706 hoil 0.001706 0.003540 ============ Information criteria: AIC: -10.37776 BIC: -10.24069 > MTSdiag(m3c) ### Model checking; see Figure 6.9 of the text [1] "Covariance matrix:" Greg hoil Greg 0.00897 0.00171 hoil 0.00171 0.00355 CCM at lag: 0 [,1] [,2] [1,] 1.000 0.302 [2,] 0.302 1.000 Simplified matrix: CCM at lag: 1 . . . . CCM at lag: 2 . . . . CCM at lag: 3 . . . . CCM at lag: 4 . . . . CCM at lag: 5 . . . . CCM at lag: 6 . . . . CCM at lag: 7 . + - . CCM at lag: 8 . . . . CCM at lag: 9 . . . . CCM at lag: 10 + . . . CCM at lag: 11 . . . . CCM at lag: 12 . . . . CCM at lag: 13 . . . . CCM at lag: 14 . . . . CCM at lag: 15 + . . . CCM at lag: 16 . . . . CCM at lag: 17 . . . . CCM at lag: 18 - . . . CCM at lag: 19 . . . . CCM at lag: 20 . + + . CCM at lag: 21 . . . . CCM at lag: 22 . . - . CCM at lag: 23 . . . . CCM at lag: 24 . . . . Hit Enter for p-value plot of individual ccm: ### Plot not shown Hit Enter to compute MQ-statistics: ### See Figure 6.10 Ljung-Box Statistics: m Q(m) df p-value [1,] 1.000 0.747 4.000 0.95 [2,] 2.000 3.088 8.000 0.93 [3,] 3.000 5.157 12.000 0.95 [4,] 4.000 8.348 16.000 0.94 [5,] 5.000 12.564 20.000 0.90 [6,] 6.000 18.516 24.000 0.78 [7,] 7.000 36.209 28.000 0.14 [8,] 8.000 37.886 32.000 0.22 [9,] 9.000 38.821 36.000 0.34 [10,] 10.000 45.625 40.000 0.25 [11,] 11.000 50.875 44.000 0.22 [12,] 12.000 53.635 48.000 0.27 [13,] 13.000 56.910 52.000 0.30 [14,] 14.000 60.033 56.000 0.33 [15,] 15.000 70.691 60.000 0.16 [16,] 16.000 75.203 64.000 0.16 [17,] 17.000 77.156 68.000 0.21 [18,] 18.000 90.496 72.000 0.07 [19,] 19.000 93.005 76.000 0.09 [20,] 20.000 105.808 80.000 0.03 [21,] 21.000 107.244 84.000 0.04 [22,] 22.000 117.233 88.000 0.02 [23,] 23.000 118.908 92.000 0.03 [24,] 24.000 119.913 96.000 0.05 Hit Enter to obtain residual plots: ## Plot not shown > ##### Missing values > da=read.table("q-gdp-ukcaus.txt",header=T) > gdp=log(da[,3:5]) > VARorder(gdp) selected order: aic = 3 selected order: bic = 2 selected order: hq = 3 Summary table: p AIC BIC HQ M(p) p-value [1,] 0 -19.3392 -19.3392 -19.3392 0.0000 0.0000 [2,] 1 -31.0894 -30.8868 -31.0071 1290.3947 0.0000 [3,] 2 -31.8832 -31.4780 -31.7186 98.8212 0.0000 [4,] 3 -31.9928 -31.3850 -31.7458 25.8717 0.0021 [5,] 4 -31.9925 -31.1822 -31.6633 14.1906 0.1157 [6,] 5 -31.9869 -30.9740 -31.5754 13.2443 0.1519 [7,] 6 -31.9153 -30.6997 -31.4214 6.6597 0.6725 [8,] 7 -31.8341 -30.4159 -31.2579 5.5790 0.7812 [9,] 8 -31.8088 -30.1880 -31.1503 10.2857 0.3279 [10,] 9 -31.8259 -30.0025 -31.0851 13.5164 0.1406 [11,] 10 -31.8339 -29.8080 -31.0108 12.2957 0.1972 [12,] 11 -31.7340 -29.5055 -30.8286 3.3733 0.9476 [13,] 12 -31.7600 -29.3289 -30.7723 12.7515 0.1742 [14,] 13 -31.7203 -29.0866 -30.6503 7.4782 0.5875 > m1=VAR(gdp,3) Constant term: Estimates: -0.05211161 -0.0163411 0.1438155 Std.Error: 0.04891077 0.05049953 0.0534254 AR coefficient matrix AR( 1 )-matrix [,1] [,2] [,3] [1,] 1.321 0.108 0.0617 [2,] 0.367 1.317 0.4519 [3,] 0.532 0.252 1.1169 standard error [,1] [,2] [,3] [1,] 0.0952 0.0973 0.0938 [2,] 0.0983 0.1005 0.0969 [3,] 0.1040 0.1063 0.1025 AR( 2 )-matrix [,1] [,2] [,3] [1,] -0.302 -0.0133 -0.0264 [2,] -0.547 -0.5135 -0.4648 [3,] -0.779 -0.3750 -0.1067 standard error [,1] [,2] [,3] [1,] 0.152 0.141 0.126 [2,] 0.157 0.146 0.130 [3,] 0.166 0.154 0.138 AR( 3 )-matrix [,1] [,2] [,3] [1,] -0.104 -0.0865 0.0273 [2,] 0.223 0.1385 0.0310 [3,] 0.380 0.0157 -0.0310 standard error [,1] [,2] [,3] [1,] 0.0933 0.0900 0.0946 [2,] 0.0963 0.0930 0.0977 [3,] 0.1019 0.0984 0.1034 Residuals cov-mtx: [,1] [,2] [,3] [1,] 2.637286e-05 2.867568e-06 7.732008e-06 [2,] 2.867568e-06 2.811401e-05 1.283314e-05 [3,] 7.732008e-06 1.283314e-05 3.146616e-05 det(SSE) = 1.761671e-14 AIC = -31.24136 BIC = -30.63358 HQ = -30.99444 > piwgt=m1$Phi > sig=m1$Sigma > cnst=m1$Ph0 > m2=Vmiss(gdp,piwgt,sig,50,cnst) Estimate of missing value at time index 50 [,1] [1,] 12.29765 [2,] 13.60911 [3,] 15.92464 > gdp[50,] uk ca us 50 12.29495 13.60897 15.92503 > m2=Vmiss(gdp,piwgt,sig,100,cnst) Estimate of missing value at time index 100 [,1] [1,] 12.67002 [2,] 14.01860 [3,] 16.33219 > gdp[100,] uk ca us 100 12.67209 14.02048 16.33217 > m3=Vpmiss(gdp,piwgt,sig,50,mdx=c(1,0,0),cnst) ### Partial missing Estimate of missing value at time index 50 The missing idicator is: 1 0 0 [,1] [1,] 13.63654 [2,] 15.88154 > m3=Vpmiss(gdp,piwgt,sig,50,mdx=c(1,0,1),cnst) Estimate of missing value at time index 50 The missing idicator is: 1 0 1 [,1] [1,] 13.60977 > #### Factor Models ############ Orthogonal factor models > da=read.table("m-tenstocks.txt",header=T) > rtn=log(da[,2:11]+1) ### Log returns > rtns=scale(rtn,center=F) ### standardize the series (no mean adjustment) > m1=princomp(rtns) > names(m1) [1] "sdev" "loadings" "center" "scale" "n.obs" "scores" "call" > sdev=m1$sdev > M=m1$loadings > summary(m1) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Standard deviation 2.1401291 1.2907427 0.96068476 0.85711091 0.75329507 Proportion of Variance 0.4625933 0.1682669 0.09321401 0.07419821 0.05731254 Cumulative Proportion 0.4625933 0.6308602 0.72407421 0.79827242 0.85558496 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10 Standard deviation 0.63443182 0.56815075 0.53461202 0.47323238 0.44136219 Proportion of Variance 0.04065269 0.03260217 0.02886668 0.02261873 0.01967477 Cumulative Proportion 0.89623765 0.92883982 0.95770650 0.98032523 1.00000000 > screeplot(m1) ### See Figure 6.12 of the text > SD=diag(sdev[1:3]) > L=M[,1:3]%*%SD > print(round(L,3)) [,1] [,2] [,3] TXN -0.788 -0.198 -0.322 MU -0.669 -0.359 -0.285 INTC -0.790 -0.184 -0.335 TSM -0.801 -0.270 -0.158 PFE -0.490 0.643 0.029 MRK -0.402 0.689 -0.226 LLY -0.447 0.697 -0.062 JPM -0.724 -0.020 0.348 MS -0.752 -0.053 0.421 GS -0.745 -0.123 0.497 > LLt=L%*%t(L) > diag(LLt) TXN MU INTC TSM PFE MRK LLY JPM 0.7631747 0.6577358 0.7694556 0.7394939 0.6542743 0.6866215 0.6895597 0.6456710 MS GS 0.7458766 0.8172213 > SigE=1-diag(LLt) > SigE TXN MU INTC TSM PFE MRK LLY JPM 0.2368253 0.3422642 0.2305444 0.2605061 0.3457257 0.3133785 0.3104403 0.3543290 MS GS 0.2541234 0.1827787 > > m2=factanal(rtns,3,scores="regression",rotation="none") > m2 Call: factanal(x = rtns, factors = 3, scores = "regression", rotation = "none") Uniquenesses: TXN MU INTC TSM PFE MRK LLY JPM MS GS 0.299 0.496 0.272 0.331 0.501 0.533 0.430 0.547 0.352 0.096 Loadings: Factor1 Factor2 Factor3 TXN 0.714 -0.132 0.417 MU 0.579 -0.263 0.315 INTC 0.718 -0.130 0.442 TSM 0.725 -0.186 0.330 PFE 0.351 0.580 0.197 MRK 0.282 0.582 0.220 LLY 0.359 0.638 0.184 JPM 0.670 MS 0.789 -0.154 GS 0.879 -0.363 Factor1 Factor2 Factor3 SS loadings 4.061 1.228 0.854 Proportion Var 0.406 0.123 0.085 Cumulative Var 0.406 0.529 0.614 Test of the hypothesis that 3 factors are sufficient. The chi square statistic is 59.43 on 18 degrees of freedom. The p-value is 2.53e-06 > names(m2) [1] "converged" "loadings" "uniquenesses" "correlation" "criteria" [6] "factors" "dof" "method" "scores" "STATISTIC" [11] "PVAL" "n.obs" "call" > L2=matrix(m2$loadings,10,3) > print(round(L2,3)) [,1] [,2] [,3] [1,] 0.714 -0.132 0.417 [2,] 0.579 -0.263 0.315 [3,] 0.718 -0.130 0.442 [4,] 0.725 -0.186 0.330 [5,] 0.351 0.580 0.197 [6,] 0.282 0.582 0.220 [7,] 0.359 0.638 0.184 [8,] 0.670 0.067 0.010 [9,] 0.789 0.052 -0.154 [10,] 0.879 -0.008 -0.363 ###### Forecasting via Diffusion index > da=read.table("m-unempstatesAdj.txt") ## Some version of the data file ## contains state headings. > dim(da) [1] 416 50 > drate=diffM(da) > dim(drate) [1] 415 50 > y=drate[5:415,1] > length(y) [1] 411 > x=cbind(drate[4:414,],drate[3:413,],drate[2:412,],drate[1:411,]) > dim(x) [1] 411 200 > m1=SWfore(y,x,350,10) MSE of out-of-sample forecasts: 0.01510996 > m1=SWfore(y,x,350,20) MSE of out-of-sample forecasts: 0.01274754 > m1=SWfore(y,x,350,30) MSE of out-of-sample forecasts: 0.01136177 > m1=SWfore(y,x,350,50) MSE of out-of-sample forecasts: 0.01044645 > m1=SWfore(y,x,350,200) MSE of out-of-sample forecasts: 0.02750807 > #### Constrained factor models > da=read.table("m-tenstocks.txt",header=T) > rtn=log(da[,2:11]+1) > h1=c(rep(1,4),rep(0,6)) > h2=c(rep(0,4),1,1,1,rep(0,3)) > h3=c(rep(0,7),rep(1,3)) > H=cbind(h1,h2,h3) > m1=hfactor(rtn,H,3) [1] "Data are individually standardized" [1] "First r eigenvalues of the correlation matrix:" [1] 4.6256602 1.6827255 0.9320882 [1] "Variability explained: " [1] 0.7240474 [1] "Loadings:" [,1] [,2] [,3] [1,] -0.368 -0.1532 -0.3331 [2,] -0.313 -0.2792 -0.3000 [3,] -0.368 -0.1419 -0.3465 [4,] -0.374 -0.2090 -0.1649 [5,] -0.229 0.4978 0.0278 [6,] -0.188 0.5333 -0.2354 [7,] -0.209 0.5401 -0.0632 [8,] -0.338 -0.0153 0.3586 [9,] -0.352 -0.0411 0.4421 [10,] -0.348 -0.0946 0.5173 [1] "eigenvalues of constrained part:" [1] 4.576 1.650 0.883 [1] "Omega-Hat" [,1] [,2] [,3] [1,] 0.761 0.2556 0.269 [2,] 0.444 -0.6752 0.101 [3,] 0.738 0.0547 -0.431 [1] "Variation explained by the constrained factors:" [1] 0.7055665 [1] "H*Omega: constrained loadings" [,1] [,2] [,3] [1,] 0.357 0.1997 0.287 [2,] 0.357 0.1997 0.287 [3,] 0.357 0.1997 0.287 [4,] 0.357 0.1997 0.287 [5,] 0.208 -0.5276 0.108 [6,] 0.208 -0.5276 0.108 [7,] 0.208 -0.5276 0.108 [8,] 0.346 0.0427 -0.460 [9,] 0.346 0.0427 -0.460 [10,] 0.346 0.0427 -0.460 [1] "Psi:" TXN MU INTC TSM PFE MRK LLY JPM TXN 0.2830 -0.1252 0.03759 -0.0881 0.04347 -0.0110 0.09870 -0.03267 MU -0.1252 0.2830 -0.19034 -0.0911 -0.00830 -0.0406 -0.20567 -0.03115 INTC 0.0376 -0.1903 0.28296 -0.0209 0.00673 0.0388 0.10487 -0.04752 TSM -0.0881 -0.0911 -0.02089 0.2830 0.04933 -0.0275 -0.01351 0.12385 PFE 0.0435 -0.0083 0.00673 0.0493 0.33689 -0.1574 -0.15304 0.17704 MRK -0.0110 -0.0406 0.03882 -0.0275 -0.15741 0.3369 -0.14916 -0.07708 LLY 0.0987 -0.2057 0.10487 -0.0135 -0.15304 -0.1492 0.33689 -0.00743 JPM -0.0327 -0.0311 -0.04752 0.1239 0.17704 -0.0771 -0.00743 0.26728 MS 0.0334 -0.0582 0.04256 0.0334 0.04181 -0.0196 0.02600 -0.21755 GS 0.0213 -0.0695 0.01492 0.0542 -0.03017 -0.0805 0.00405 -0.15039 MS GS TXN 0.0334 0.02130 MU -0.0582 -0.06949 INTC 0.0426 0.01492 TSM 0.0334 0.05421 PFE 0.0418 -0.03017 MRK -0.0196 -0.08052 LLY 0.0260 0.00405 JPM -0.2176 -0.15039 MS 0.2673 0.01756 GS 0.0176 0.26728 [1] "Diagonal elements of Psi:" TXN MU INTC TSM PFE MRK LLY JPM MS GS 0.283 0.283 0.283 0.283 0.337 0.337 0.337 0.267 0.267 0.267 [1] "eigenvalues of Psi:" [1] 0.7632 0.6350 0.4539 0.3968 0.3741 0.2474 0.2051 0.0262 -0.0670 [10] -0.0902 > ###### Asymptotic Principal Component Analysis > da=read.table("m-sp100y2011.txt",header=T) > sp100=apca(da,3) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Standard deviation 0.2243715 0.1333795 0.09154043 0.05873897 0.05594621 Proportion of Variance 0.5298019 0.1872217 0.08818687 0.03631037 0.03293967 Cumulative Proportion 0.5298019 0.7170236 0.80521047 0.84152084 0.87446052 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10 Standard deviation 0.05387705 0.0495897 0.04381794 0.03738959 0.03489291 Proportion of Variance 0.03054819 0.0258798 0.02020607 0.01471226 0.01281305 Cumulative Proportion 0.90500871 0.9308885 0.95109459 0.96580685 0.97861990 Comp.11 Comp.12 Standard deviation 0.03228930 0.03144792 Proportion of Variance 0.01097224 0.01040787 Cumulative Proportion 0.98959213 1.00000000 > screeplot(sp100) ### See Figure 6.16 > factors=sp100$factors > print(round(factors,4)) [,1] [,2] [,3] [1,] -0.0136 -0.1309 -0.2613 [2,] -0.0258 0.2105 -0.9075 [3,] -0.0206 0.0045 0.1139 [4,] 0.0079 0.1165 -0.0343 [5,] 0.9963 0.0505 -0.0107 [6,] 0.0118 0.0204 0.0208 [7,] -0.0228 -0.0139 0.0384 [8,] -0.0337 0.4671 0.2725 [9,] -0.0208 0.6372 -0.0512 [10,] 0.0167 -0.4830 -0.1069 [11,] -0.0586 0.1457 0.0457 [12,] -0.0117 0.2076 0.0390 > loadings=sp100$loadings > which.max(loadings[,1]) C 71 > which.max(loadings[,3]) FCX 81 > par(mfcol=c(3,1)) > plot(loadings[,1],xlab="Company index",ylab="Loadings-1",type="l") > plot(loadings[,2],xlab="Company index",ylab="Loadings-2",type="l") > plot(loadings[,3],xlab="Company index",ylab="Loadings-3",type="l") #### See Figure 6.17 of the text ##### Model-based clustering > source("MBcluster.R") > da=read.table("m-unempstatesAdj.txt",header=T) > dim(da) [1] 416 50 > da[1:5,1:3] AL AK AZ 1 6.4 7.1 7.1087 2 6.3 7.0 6.9087 3 6.1 7.0 6.6087 4 6.0 7.0 6.4087 5 6.0 7.0 6.2087 > zt=apply(da,2,diff) > dim(zt) [1] 415 50 > mcmc=list(burnin=1000,rep=1000) > m1=MBcluster(zt,4,3,mcmc=mcmc) Use default priors Estimation for Cluster Parameters: Number of Clusters: K= 3 Number of Lags in AR model: p= 4 phi 0 phi 1 phi 2 phi 3 phi 4 sigma Cluster 1 0.0020743 0.31486 0.29799 0.124271 -0.025789 0.11434 Cluster 2 0.0010702 0.61878 0.43858 -0.078071 -0.195115 0.07593 Cluster 3 0.0019771 0.38924 0.26979 0.039070 -0.056382 0.14988 Classification Probabilities: Cluster 1 Cluster 2 Cluster 3 AL 0.000 1 0.000 AK 0.000 1 0.000 AZ 0.964 0 0.036 AR 0.000 1 0.000 CA 0.000 1 0.000 CO 0.000 1 0.000 CT 0.000 1 0.000 DE 0.003 0 0.997 FL 0.964 0 0.036 GA 0.003 0 0.997 HI 0.000 1 0.000 ID 0.000 1 0.000 IL 0.000 1 0.000 IN 0.003 0 0.997 IA 0.964 0 0.036 KS 0.000 1 0.000 KY 0.000 1 0.000 LA 0.000 1 0.000 ME 0.000 1 0.000 MD 0.964 0 0.036 MA 0.000 1 0.000 MI 0.000 1 0.000 MN 0.000 1 0.000 MS 0.008 0 0.992 MO 0.000 1 0.000 MT 0.000 1 0.000 NE 0.000 1 0.000 NV 0.000 1 0.000 NH 0.000 1 0.000 NJ 0.000 1 0.000 NM 0.000 1 0.000 NY 0.000 1 0.000 NC 0.003 0 0.997 ND 0.000 1 0.000 OH 0.003 0 0.997 OK 0.003 0 0.997 OR 0.000 1 0.000 PA 0.964 0 0.036 RI 0.000 1 0.000 SC 0.000 1 0.000 SD 0.000 1 0.000 TN 0.000 1 0.000 TX 0.964 0 0.036 UT 0.000 1 0.000 VT 0.000 1 0.000 VA 0.000 1 0.000 WA 0.000 1 0.000 WV 0.000 1 0.000 WI 0.003 0 0.997 WY 0.000 1 0.000 Classification: Cluster 1 : Number of members: 6 AZ FL IA MD PA TX Cluster 2 : Number of members: 36 AL AK AR CA CO CT HI ID IL KS KY LA ME MA MI MN MO MT NE NV NH NJ NM NY ND OR RI SC SD TN UT VT VA WA WV WY Cluster 3 : Number of members: 8 DE GA IN MS NC OH OK WI Marginal LogLikelihood: 20641.03 > > m2=MBcluster(zt,4,4,mcmc=mcmc) Use default priors Estimation for Cluster Parameters: Number of Clusters: K= 4 Number of Lags in AR model: p= 4 phi 0 phi 1 phi 2 phi 3 phi 4 sigma Cluster 1 0.00133262 0.71090 0.39097 -0.134111 -0.186601 0.082931 Cluster 2 0.00081442 0.51149 0.47222 -0.010046 -0.178475 0.070190 Cluster 3 0.00205418 0.31424 0.29864 0.125113 -0.025592 0.114164 Cluster 4 0.00223332 0.38992 0.26965 0.037280 -0.057303 0.150605 Classification Probabilities: Cluster 1 Cluster 2 Cluster 3 Cluster 4 AL 1.000 0.000 0.000 0.000 AK 0.492 0.508 0.000 0.000 AZ 0.000 0.000 0.989 0.011 AR 0.000 1.000 0.000 0.000 CA 0.000 1.000 0.000 0.000 CO 1.000 0.000 0.000 0.000 CT 0.000 1.000 0.000 0.000 DE 0.000 0.000 0.001 0.999 FL 0.000 0.000 0.989 0.011 GA 0.000 0.000 0.001 0.999 HI 0.995 0.005 0.000 0.000 ID 0.000 1.000 0.000 0.000 IL 1.000 0.000 0.000 0.000 IN 0.000 0.000 0.001 0.999 IA 0.000 0.000 0.989 0.011 KS 0.000 1.000 0.000 0.000 KY 0.890 0.110 0.000 0.000 LA 1.000 0.000 0.000 0.000 ME 1.000 0.000 0.000 0.000 MD 0.000 0.000 0.989 0.011 MA 0.000 1.000 0.000 0.000 MI 1.000 0.000 0.000 0.000 MN 0.000 1.000 0.000 0.000 MS 0.000 0.000 0.000 1.000 MO 0.999 0.001 0.000 0.000 MT 0.000 1.000 0.000 0.000 NE 0.000 1.000 0.000 0.000 NV 1.000 0.000 0.000 0.000 NH 0.000 1.000 0.000 0.000 NJ 0.000 1.000 0.000 0.000 NM 0.000 1.000 0.000 0.000 NY 0.491 0.509 0.000 0.000 NC 0.000 0.000 0.001 0.999 ND 0.000 1.000 0.000 0.000 OH 0.000 0.000 0.001 0.999 OK 0.000 0.000 0.001 0.999 OR 1.000 0.000 0.000 0.000 PA 0.000 0.000 0.989 0.011 RI 0.000 1.000 0.000 0.000 SC 1.000 0.000 0.000 0.000 SD 0.000 1.000 0.000 0.000 TN 0.240 0.760 0.000 0.000 TX 0.000 0.000 0.989 0.011 UT 0.033 0.967 0.000 0.000 VT 1.000 0.000 0.000 0.000 VA 0.024 0.976 0.000 0.000 WA 0.004 0.996 0.000 0.000 WV 1.000 0.000 0.000 0.000 WI 0.000 0.000 0.001 0.999 WY 0.278 0.722 0.000 0.000 Classification: Cluster 1 : Number of members: 14 AL CO HI IL KY LA ME MI MO NV OR SC VT WV Cluster 2 : Number of members: 22 AK AR CA CT ID KS MA MN MT NE NH NJ NM NY ND RI SD TN UT VA WA WY Cluster 3 : Number of members: 6 AZ FL IA MD PA TX Cluster 4 : Number of members: 8 DE GA IN MS NC OH OK WI Marginal LogLikelihood: 20806.63 > > mcmc=list(burnin=2000,rep=1000) > m3=MBcluster(zt,4,4,mcmc=mcmc) Use default priors Estimation for Cluster Parameters: Number of Clusters: K= 4 Number of Lags in AR model: p= 4 phi 0 phi 1 phi 2 phi 3 phi 4 sigma Cluster 1 0.00128788 0.69373 0.40134 -0.124572 -0.190477 0.081734 Cluster 2 0.00214846 0.31495 0.29855 0.124724 -0.025586 0.114190 Cluster 3 0.00230516 0.38906 0.26896 0.040165 -0.055680 0.149503 Cluster 4 0.00080644 0.49235 0.47495 0.004314 -0.171911 0.069182 Classification Probabilities: Cluster 1 Cluster 2 Cluster 3 Cluster 4 AL 1.000 0.000 0.000 0.000 AK 0.979 0.000 0.000 0.021 AZ 0.000 0.949 0.051 0.000 AR 0.000 0.000 0.000 1.000 CA 0.000 0.000 0.000 1.000 CO 1.000 0.000 0.000 0.000 CT 0.000 0.000 0.000 1.000 DE 0.000 0.000 1.000 0.000 FL 0.000 0.949 0.051 0.000 GA 0.000 0.000 1.000 0.000 HI 1.000 0.000 0.000 0.000 ID 0.000 0.000 0.000 1.000 IL 1.000 0.000 0.000 0.000 IN 0.000 0.000 1.000 0.000 IA 0.000 0.949 0.051 0.000 KS 0.000 0.000 0.000 1.000 KY 0.997 0.000 0.000 0.003 LA 1.000 0.000 0.000 0.000 ME 1.000 0.000 0.000 0.000 MD 0.000 0.949 0.051 0.000 MA 0.000 0.000 0.000 1.000 MI 1.000 0.000 0.000 0.000 MN 0.002 0.000 0.000 0.998 MS 0.000 0.007 0.993 0.000 MO 1.000 0.000 0.000 0.000 MT 0.000 0.000 0.000 1.000 NE 0.003 0.000 0.000 0.997 NV 1.000 0.000 0.000 0.000 NH 0.000 0.000 0.000 1.000 NJ 0.000 0.000 0.000 1.000 NM 0.000 0.000 0.000 1.000 NY 0.955 0.000 0.000 0.045 NC 0.000 0.000 1.000 0.000 ND 0.000 0.000 0.000 1.000 OH 0.000 0.000 1.000 0.000 OK 0.000 0.000 1.000 0.000 OR 1.000 0.000 0.000 0.000 PA 0.000 0.949 0.051 0.000 RI 0.000 0.000 0.000 1.000 SC 1.000 0.000 0.000 0.000 SD 0.000 0.000 0.000 1.000 TN 0.856 0.000 0.000 0.144 TX 0.000 0.949 0.051 0.000 UT 0.410 0.000 0.000 0.590 VT 1.000 0.000 0.000 0.000 VA 0.389 0.000 0.000 0.611 WA 0.096 0.000 0.000 0.904 WV 1.000 0.000 0.000 0.000 WI 0.000 0.000 1.000 0.000 WY 0.902 0.000 0.000 0.098 Classification: Cluster 1 : Number of members: 18 AL AK CO HI IL KY LA ME MI MO NV NY OR SC TN VT WV WY Cluster 2 : Number of members: 6 AZ FL IA MD PA TX Cluster 3 : Number of members: 8 DE GA IN MS NC OH OK WI Cluster 4 : Number of members: 18 AR CA CT ID KS MA MN MT NE NH NJ NM ND RI SD UT VA WA Marginal LogLikelihood: 20814.22 > > mcmc=list(burnin=3000,rep=2000) ##### > m4=MBcluster(zt,4,4,mcmc=mcmc) Use default priors Estimation for Cluster Parameters: Number of Clusters: K= 4 Number of Lags in AR model: p= 4 phi 0 phi 1 phi 2 phi 3 phi 4 sigma Cluster 1 0.00127731 0.69889 0.39798 -0.12770645 -0.188757 0.082096 Cluster 2 0.00080087 0.49824 0.47396 -0.00050977 -0.173598 0.069493 Cluster 3 0.00197971 0.31440 0.29790 0.12533160 -0.025974 0.114214 Cluster 4 0.00226912 0.39076 0.26922 0.03767486 -0.057572 0.150859 Classification Probabilities: Cluster 1 Cluster 2 Cluster 3 Cluster 4 AL 0.9995 0.0005 0.0000 0.0000 AK 0.9010 0.0990 0.0000 0.0000 AZ 0.0000 0.0000 0.9945 0.0055 AR 0.0000 1.0000 0.0000 0.0000 CA 0.0000 1.0000 0.0000 0.0000 CO 1.0000 0.0000 0.0000 0.0000 CT 0.0000 1.0000 0.0000 0.0000 DE 0.0000 0.0000 0.0010 0.9990 FL 0.0000 0.0000 0.9945 0.0055 GA 0.0000 0.0000 0.0010 0.9990 HI 0.9995 0.0005 0.0000 0.0000 ID 0.0000 1.0000 0.0000 0.0000 IL 1.0000 0.0000 0.0000 0.0000 IN 0.0000 0.0000 0.0010 0.9990 IA 0.0000 0.0000 0.9945 0.0055 KS 0.0000 1.0000 0.0000 0.0000 KY 0.9890 0.0110 0.0000 0.0000 LA 1.0000 0.0000 0.0000 0.0000 ME 1.0000 0.0000 0.0000 0.0000 MD 0.0000 0.0000 0.9945 0.0055 MA 0.0000 1.0000 0.0000 0.0000 MI 1.0000 0.0000 0.0000 0.0000 MN 0.0005 0.9995 0.0000 0.0000 MS 0.0000 0.0000 0.0015 0.9985 MO 1.0000 0.0000 0.0000 0.0000 MT 0.0000 1.0000 0.0000 0.0000 NE 0.0005 0.9995 0.0000 0.0000 NV 1.0000 0.0000 0.0000 0.0000 NH 0.0000 1.0000 0.0000 0.0000 NJ 0.0000 1.0000 0.0000 0.0000 NM 0.0000 1.0000 0.0000 0.0000 NY 0.8710 0.1290 0.0000 0.0000 NC 0.0000 0.0000 0.0010 0.9990 ND 0.0000 1.0000 0.0000 0.0000 OH 0.0000 0.0000 0.0010 0.9990 OK 0.0000 0.0000 0.0010 0.9990 OR 1.0000 0.0000 0.0000 0.0000 PA 0.0000 0.0000 0.9945 0.0055 RI 0.0000 1.0000 0.0000 0.0000 SC 1.0000 0.0000 0.0000 0.0000 SD 0.0000 1.0000 0.0000 0.0000 TN 0.6730 0.3270 0.0000 0.0000 TX 0.0000 0.0000 0.9945 0.0055 UT 0.2175 0.7825 0.0000 0.0000 VT 1.0000 0.0000 0.0000 0.0000 VA 0.2080 0.7920 0.0000 0.0000 WA 0.0320 0.9680 0.0000 0.0000 WV 1.0000 0.0000 0.0000 0.0000 WI 0.0000 0.0000 0.0010 0.9990 WY 0.7775 0.2225 0.0000 0.0000 Classification: Cluster 1 : Number of members: 18 AL AK CO HI IL KY LA ME MI MO NV NY OR SC TN VT WV WY Cluster 2 : Number of members: 18 AR CA CT ID KS MA MN MT NE NH NJ NM ND RI SD UT VA WA Cluster 3 : Number of members: 6 AZ FL IA MD PA TX Cluster 4 : Number of members: 8 DE GA IN MS NC OH OK WI Marginal LogLikelihood: 20812.45