### R output of HW#1 ### Problem 1 > p1 <- c(1,-.4,-.6) > s1 <- polyroot(p1) > s1 [1] 1.000000-0i -1.666667+0i > p2 <- c(1,-1.3,.4) > s2 <- polyroot(p2) > s2 [1] 1.25-0i 2.00+0i > (3.5+0.5)/(1-1.3+0.4) [1] 40 #### Problem 2 > source("CornerFun.R") > omega=c(0,3,-1) > delta=c(0.8) > CornerFun(omega,delta,5,5) Impulse responses: 0 3 1.4 1.12 0.896 0.7168 0.57344 0.458752 0.3670016 0.2936013 0.234881 0.1879048 0.1503239 0.1202591 0.09620727 0.07696581 u: 0 1 0.4666667 0.3733333 0.2986667 0.2389333 0.1911467 0.1529173 0.1223339 0.09786709 0.07829367 0.06263494 0.05010795 0.04008636 0.03206909 0.02565527 Corner Table [,1] [,2] [,3] [,4] [,5] [1,] 0.00000 0.00000 0.00000 0.00000 0.00000 [2,] 1.00000 1.00000 1.00000 1.00000 1.00000 [3,] 0.46667 -0.15556 0.05185 -0.01728 0.00576 [4,] 0.37333 0.00000 0.00000 0.00000 0.00000 [5,] 0.29867 0.00000 0.00000 0.00000 0.00000 #### Problem 3 > da <- read.table("m-unemini-6717.txt") > fix(da) ## spreadsheet of the data > da[,1] <- da[,1]/10000 ### re-scale initial claims so that the two series ##################################### are of similar scale > colnames(da) <- c("ini","unem") > ar(da[,1],method="mle") Call: ar(x = da[, 1], method = "mle") Coefficients: 1 2 3 4 5 6 7 1.1009 -0.0792 -0.0364 0.0455 -0.0971 0.1456 -0.1119 Order selected 7 sigma^2 estimated as 3.312 ### Use the command "ar" to identify a marginal model for initial claims > m1 <- arima(da$ini,order=c(7,0,0)) > m1 Call: arima(x = da$ini, order = c(7, 0, 0)) Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 intercept 1.1010 -0.0792 -0.0364 0.0455 -0.0971 0.1456 -0.1119 34.9001 s.e. 0.0405 0.0602 0.0602 0.0602 0.0603 0.0603 0.0407 2.2345 sigma^2 estimated as 3.312: log likelihood = -1216.28, aic = 2450.56 > c1 <- c(NA,0,0,0,NA,NA,NA,NA) #### refine the fitted model > m1a <- arima(da$ini,order=c(7,0,0),fixed=c1) Warning message: In arima(da$ini, order = c(7, 0, 0), fixed = c1) : some AR parameters were fixed: setting transform.pars = FALSE > m1a Call: arima(x = da$ini, order = c(7, 0, 0), fixed = c1) Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 intercept 1.0335 0 0 0 -0.0995 0.1415 -0.1082 34.8844 s.e. 0.0186 0 0 0 0.0462 0.0603 0.0407 2.2271 sigma^2 estimated as 3.335: log likelihood = -1218.34, aic = 2448.67 > tsdiag(m1a,gof=12) ### Not shown, but it looks the model fits the data well. > > ar(da$unem) ### Unemployment rate Call: ar(x = da$unem) Coefficients: 1 2 3 4 5 6 7 1.0393 0.0773 -0.0045 -0.0233 -0.0375 0.0157 -0.0836 Order selected 7 sigma^2 estimated as 0.03913 > m2 <- arima(da$unem,order=c(7,0,0)) > m2 Call: arima(x = da$unem, order = c(7, 0, 0)) Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 intercept 1.0283 0.1135 0.0000 -0.0058 -0.0560 -0.0249 -0.0692 6.0481 s.e. 0.0406 0.0584 0.0585 0.0586 0.0585 0.0584 0.0407 0.4509 sigma^2 estimated as 0.02581: log likelihood = 244.07, aic = -470.14 > c2 <- c(NA,NA,0,0,0,0,NA,NA) > m2a <- arima(da$unem,order=c(7,0,0),fixed=c2) Warning message: In arima(da$unem, order = c(7, 0, 0), fixed = c2) : some AR parameters were fixed: setting transform.pars = FALSE > m2a Call: arima(x = da$unem, order = c(7, 0, 0), fixed = c2) Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 intercept 1.0288 0.0818 0 0 0 0 -0.1253 6.0619 s.e. 0.0405 0.0469 0 0 0 0 0.0133 0.4321 sigma^2 estimated as 0.02595: log likelihood = 242.45, aic = -474.9 > tsdiag(m2a,gof=12) ## Not shown,the fit looks fine. > require(MTS) Loading required package: MTS > VARorder(da) selected order: aic = 12 selected order: bic = 2 selected order: hq = 7 Summary table: p AIC BIC HQ M(p) p-value [1,] 0 4.2793 4.2793 4.2793 0.0000 0.0000 [2,] 1 -2.5570 -2.5278 -2.5457 4010.4707 0.0000 [3,] 2 -2.6119 -2.5534 -2.5892 39.7786 0.0000 [4,] 3 -2.6363 -2.5486 -2.6022 21.9179 0.0002 [5,] 4 -2.6510 -2.5340 -2.6055 16.1952 0.0028 [6,] 5 -2.6619 -2.5157 -2.6050 13.9799 0.0074 [7,] 6 -2.6861 -2.5107 -2.6179 21.6062 0.0002 [8,] 7 -2.7081 -2.5034 -2.6284 20.1871 0.0005 [9,] 8 -2.7039 -2.4700 -2.6129 5.2345 0.2641 [10,] 9 -2.7112 -2.4481 -2.6088 11.7209 0.0196 [11,] 10 -2.7125 -2.4202 -2.5987 8.2832 0.0817 [12,] 11 -2.7017 -2.3801 -2.5765 1.3974 0.8446 [13,] 12 -2.7227 -2.3718 -2.5861 19.3104 0.0007 [14,] 13 -2.7171 -2.3370 -2.5691 4.3085 0.3659 > GrangerTest(da,p=7) ### Use VAR(7) model Number of targeted zero parameters: 7 Chi-square test for Granger Causality and p-value: 19.64567 0.006388102 > > GrangerTest(da,p=7,locInput=c(2)) Number of targeted zero parameters: 7 Chi-square test for Granger Causality and p-value: 138.4576 0 > ##### Problem 4 > gas <- read.csv("w-gasreg.csv") > oil <- read.csv("w-wti.csv") > fix(oil) > x <- cbind(gas[,2],oil[,2]) > y <- log(x) > VARorder(y) selected order: aic = 9 selected order: bic = 3 selected order: hq = 3 Summary table: p AIC BIC HQ M(p) p-value [1,] 0 -6.0290 -6.0290 -6.0290 0.0000 0.0000 [2,] 1 -14.5442 -14.5289 -14.5385 11499.1077 0.0000 [3,] 2 -14.9334 -14.9029 -14.9220 532.3908 0.0000 [4,] 3 -14.9613 -14.9154 -14.9441 45.3161 0.0000 [5,] 4 -14.9665 -14.9054 -14.9436 14.9105 0.0049 [6,] 5 -14.9686 -14.8922 -14.9400 10.7319 0.0297 [7,] 6 -14.9724 -14.8806 -14.9380 12.8172 0.0122 [8,] 7 -14.9707 -14.8637 -14.9306 5.5957 0.2314 [9,] 8 -14.9666 -14.8443 -14.9208 2.3435 0.6729 [10,] 9 -14.9758 -14.8382 -14.9243 20.0801 0.0005 [11,] 10 -14.9744 -14.8216 -14.9172 6.0047 0.1988 [12,] 11 -14.9701 -14.8020 -14.9072 2.0007 0.7356 [13,] 12 -14.9725 -14.7891 -14.9039 11.0324 0.0262 [14,] 13 -14.9671 -14.7684 -14.8927 0.5507 0.9684 > GrangerTest(y,p=3) Number of targeted zero parameters: 3 Chi-square test for Granger Causality and p-value: 56.11119 3.977596e-12 > GrangerTest(y,p=3,locInput=c(2)) Number of targeted zero parameters: 3 Chi-square test for Granger Causality and p-value: 42.05366 3.908289e-09 > #### Problem 5 > C= matrix(c(0.8,-0.3,0.4,0.6),2,2) > S <- matrix(c(2.0,0.5,0.5,1.0),2,2) > C [,1] [,2] [1,] 0.8 0.4 [2,] -0.3 0.6 > S [,1] [,2] [1,] 2.0 0.5 [2,] 0.5 1.0 > m2 <- VARMAsim(200,arlags=c(1),phi=C,sigma=S) > zt <- m2$series > MTSplot(zt) > cc <- ccm(zt,lag=3) [1] "Covariance matrix:" [,1] [,2] [1,] 8.22 -0.93 [2,] -0.93 3.24 CCM at lag: 0 [,1] [,2] [1,] 1.00 -0.18 [2,] -0.18 1.00 Simplified matrix: CCM at lag: 1 + + - + CCM at lag: 2 + + - + CCM at lag: 3 + + - . Hit Enter for p-value plot of individual ccm: > names(cc) [1] "ccm" "pvalue" > cc$ccm [,1] [,2] [,3] [,4] [1,] 1.000000 0.8218888 0.5574345 0.29482696 [2,] -0.180293 -0.5590669 -0.7142569 -0.66169977 [3,] -0.180293 0.1485978 0.3311202 0.36215043 [4,] 1.000000 0.6757953 0.3384625 0.01572472 > mq(zt,lag=5) Ljung-Box Statistics: m Q(m) df p-value [1,] 1 277 4 0 [2,] 2 477 8 0 [3,] 3 610 12 0 [4,] 4 690 16 0 [5,] 5 730 20 0 >