## Selected R commands and output used in Chapter 5 >require(glarma) > data("Asthma",package="glarma") > y <- Asthma[,1] > X <- as.matrix(Asthma[,-1]) > z <- sqrt(y) > m1 <- lm(z~-1+.,data=data.frame(X)) > summary(m1) Call: lm(formula = z ~ -1 + ., data = data.frame(X)) Coefficients: Estimate Std. Error t value Pr(>|t|) Intercept 1.19440 0.05177 23.072 < 2e-16 *** Sunday 0.12873 0.04855 2.652 0.008095 ** Monday 0.18631 0.04776 3.901 0.000100 *** CosAnnual -0.16244 0.03022 -5.376 8.88e-08 *** SinAnnual 0.13169 0.03272 4.025 5.99e-05 *** H7 0.10311 0.04667 2.210 0.027292 * NO2max -0.08138 0.02761 -2.947 0.003256 ** T1.1990 0.20539 0.05374 3.822 0.000138 *** T2.1990 0.11551 0.05321 2.171 0.030124 * T1.1991 0.05616 0.05481 1.025 0.305715 T2.1991 0.15497 0.05457 2.840 0.004579 ** T1.1992 0.21848 0.05383 4.059 5.19e-05 *** T2.1992 0.28513 0.05315 5.364 9.46e-08 *** T1.1993 0.37351 0.05393 6.926 6.48e-12 *** T2.1993 0.06352 0.05461 1.163 0.244928 --- Residual standard error: 0.6296 on 1446 degrees of freedom Multiple R-squared: 0.7977, Adjusted R-squared: 0.7956 F-statistic: 380 on 15 and 1446 DF, p-value: < 2.2e-16 #### > m2 <-arima(z,seasonal=list(order=c(0,0,1),period=7),xreg=X,include.mean=F) > m2 Call:arima(x=z,seasonal=list(order=c(0,0,1),period=7),xreg=X,include.mean=F) Coefficients: sma1 Intercept Sunday Monday CosAnnual SinAnnual H7 NO2max 0.0648 1.1899 0.1296 0.1870 -0.1607 0.1314 0.1009 -0.0789 s.e. 0.0269 0.0524 0.0512 0.0505 0.0315 0.0344 0.0470 0.0277 T1.1990 T2.1990 T1.1991 T2.1991 T1.1992 T2.1992 T1.1993 T2.1993 0.2076 0.116 0.0561 0.1565 0.2189 0.2866 0.3731 0.0621 s.e. 0.0565 0.056 0.0575 0.0573 0.0566 0.0559 0.0567 0.0573 sigma^2 estimated as 0.3908: log likelihood = -1386.74, aic = 2807.48 > Box.test(m2$residuals,lag=400,type="Ljung") Box-Ljung test data: m2$residuals X-squared = 332.02, df = 400, p-value = 0.9943 #### m4 <- glarma(y,X,thetaLags=7,type="NegBin",method="NR",alphaInit=0, maxit=100,grad=1e-6,residuals="Pearson") > summary(m4) Call: glarma(y = y,X=X,type="NegBin", method="NR", residuals = "Pearson", thetaLags = 7, alphaInit = 0, maxit = 100, grad = 1e-06) Negative Binomial Parameter: Estimate Std.Error z-ratio Pr(>|z|) alpha 37.19 25.44 1.462 0.144 GLARMA Coefficients: Estimate Std.Error z-ratio Pr(>|z|) theta_7 0.04392 0.01936 2.269 0.0233 * Linear Model Coefficients: Estimate Std.Error z-ratio Pr(>|z|) Intercept 0.58397 0.06331 9.225 < 2e-16 *** Sunday 0.19455 0.05760 3.377 0.000732 *** Monday 0.22999 0.05642 4.076 4.57e-05 *** CosAnnual -0.21450 0.03965 -5.410 6.31e-08 *** SinAnnual 0.17728 0.04153 4.269 1.96e-05 *** H7 0.16843 0.05634 2.990 0.002794 ** NO2max -0.10404 0.03392 -3.067 0.002163 ** T1.1990 0.19903 0.05845 3.405 0.000662 *** T2.1990 0.13087 0.05897 2.219 0.026477 * T1.1991 0.08587 0.06746 1.273 0.203058 T2.1991 0.17082 0.05951 2.871 0.004096 ** T1.1992 0.25276 0.05669 4.459 8.24e-06 *** T2.1992 0.30572 0.05103 5.991 2.08e-09 *** T1.1993 0.43607 0.05233 8.334 < 2e-16 *** T2.1993 0.11412 0.06269 1.821 0.068679 . Null deviance: 1989.9 on 1460 degrees of freedom Residual deviance: 1442.6 on 1444 degrees of freedom AIC: 4873.511 Number of Newton Raphson iterations: 6 LRT and Wald Test: Alternative hypothesis: model is a GLARMA process Null hypothesis: model is a GLM with the same regression structure Statistic p-value LR Test 7.047 0.00794 ** Wald Test 5.147 0.02329 * --- > par(mfcol=c(2,2)) > plot(m4,which=c(1,2,3,5),titles=c("","","","PIT for glarma(NegBin)")) #### > da <- read.table(``taq-JBES2011.txt'',header=T) > y <- da$GLT > require(NTS) > X <- hfDummy(5,Fopen=4,Tend=4,days=61,skipmin=15) > X <- as.matrix(X) > colnames(X) <- paste(``X'',1:8,sep=``'') > m1 <- ACMx(y,X=X) ## Output omitted. Use default order (1,1). > names(m1) > acf(m1$sresi,lag=230) > m2 <- ACMx(y,X=X,cond.dist=``nb'') > m3 <- ACMx(y,X=X,cond.dist=``dp'') #### > require(fGarch) > m1 <- garchFit(~garch(1,1),data=spy,trace=F) > summary(m1) Title: GARCH Modelling Call: garchFit(formula = ~garch(1, 1), data = spy, trace=F) Conditional Distribution: norm Error Analysis: Estimate Std. Error t value Pr(>|t|) mu 6.842e-04 1.681e-04 4.070 4.70e-05 *** omega 3.494e-06 5.469e-07 6.388 1.68e-10 *** alpha1 1.275e-01 1.403e-02 9.086 < 2e-16 *** beta1 8.464e-01 1.463e-02 57.856 < 2e-16 *** --- Standardised Residuals Tests: Statistic p-Value Jarque-Bera Test R Chi^2 557.7797 0 Ljung-Box Test R Q(10) 16.93149 0.07589384 Ljung-Box Test R Q(20) 26.75075 0.1424233 Ljung-Box Test R^2 Q(10) 14.30636 0.1594707 Ljung-Box Test R^2 Q(20) 20.29979 0.4393226 > plot(m1) > m2 <- garchFit(~garch(1,1),data=spy,cond.dist="std",trace=F) > summary(m2) Title: GARCH Modelling Call: garchFit(formula=~garch(1,1),data=spy,cond.dist="std",trace=F) Conditional Distribution: std Error Analysis: Estimate Std. Error t value Pr(>|t|) mu 9.031e-04 1.509e-04 5.986 2.16e-09 *** omega 2.477e-06 6.248e-07 3.965 7.34e-05 *** alpha1 1.456e-01 1.915e-02 7.601 2.95e-14 *** beta1 8.490e-01 1.681e-02 50.499 < 2e-16 *** shape 5.160e+00 5.813e-01 8.876 < 2e-16 *** --- Standardised Residuals Tests: Statistic p-Value Ljung-Box Test R Q(10) 17.45736 0.06483688 Ljung-Box Test R Q(20) 27.11551 0.1320453 Ljung-Box Test R^2 Q(10) 10.25488 0.4184234 Ljung-Box Test R^2 Q(20) 19.70427 0.4765604 LM Arch Test R TR^2 11.08226 0.5218827 Information Criterion Statistics: AIC BIC SIC HQIC -6.430537 -6.418958 -6.430545 -6.426335 > m3 <- garchFit(~garch(1,1),data=spy,cond.dist="sstd",trace=F) > summary(m3) Title: GARCH Modelling Call:garchFit(formula=~garch(1,1),data=spy,cond.dist="sstd",trace=F) Conditional Distribution: sstd Error Analysis: Estimate Std. Error t value Pr(>|t|) mu 6.441e-04 1.603e-04 4.017 5.90e-05 *** omega 2.298e-06 5.858e-07 3.922 8.77e-05 *** alpha1 1.405e-01 1.790e-02 7.848 4.22e-15 *** beta1 8.514e-01 1.628e-02 52.310 < 2e-16 *** skew 8.882e-01 2.352e-02 37.760 < 2e-16 *** shape 5.637e+00 6.844e-01 8.236 2.22e-16 *** --- Standardised Residuals Tests: Statistic p-Value Ljung-Box Test R Q(10) 17.46512 0.06468508 Ljung-Box Test R Q(20) 26.92805 0.1373005 Ljung-Box Test R^2 Q(10) 9.8221 0.4562374 Ljung-Box Test R^2 Q(20) 19.38767 0.4967704 Information Criterion Statistics: AIC BIC SIC HQIC -6.437849 -6.423954 -6.437861 -6.432807 > predict(m3,5) meanForecast meanError standardDeviation 1 0.0006440697 0.005880318 0.005880318 2 0.0006440697 0.006049519 0.006049519 3 0.0006440697 0.006212803 0.006212803 4 0.0006440697 0.006370635 0.006370635 5 0.0006440697 0.006523421 0.006523421 #### > m4 <- garchFit(~garch(1,1),data=spy,leverage=T,trace=F,cond.dist="sstd") > summary(m4) Title: GARCH Modelling Call: garchFit(formula=~garch(1,1),data spy,cond.dist="sstd",leverage=T,trace=F) Conditional Distribution: sstd Error Analysis: Estimate Std. Error t value Pr(>|t|) mu 2.457e-04 1.571e-04 1.563 0.118 omega 2.890e-06 5.379e-07 5.372 7.80e-08 *** alpha1 6.951e-02 1.646e-02 4.223 2.41e-05 *** gamma1 1.000e+00 2.018e-01 4.956 7.21e-07 *** beta1 8.492e-01 1.446e-02 58.713 < 2e-16 *** skew 8.463e-01 2.330e-02 36.321 < 2e-16 *** shape 6.083e+00 7.768e-01 7.831 4.88e-15 *** --- Standardised Residuals Tests: Statistic p-Value Ljung-Box Test R Q(10) 14.5514 0.1492931 Ljung-Box Test R Q(20) 23.74802 0.2535727 Ljung-Box Test R^2 Q(10) 13.71561 0.1863633 Ljung-Box Test R^2 Q(20) 21.74093 0.3547326 LM Arch Test R TR^2 16.33507 0.1763667 Information Criterion Statistics: AIC BIC SIC HQIC -6.484754 -6.468543 -6.484769 -6.478871 #### > require(rugarch) >spec4 <- ugarchspec(variance.model=list(model="eGARCH"),mean.model= list(armaOrder=c(0,0)),distribution.model="sstd") > m6 <- ugarchfit(data=spy,spec=spec4) > m6 *---------------------------------* * GARCH Model Fit * *---------------------------------* Conditional Variance Dynamics ----------------------------------- GARCH Model : eGARCH(1,1) Mean Model : ARFIMA(0,0,0) Distribution : sstd Optimal Parameters ------------------------------------ Estimate Std. Error t value Pr(>|t|) mu 0.000197 0.000295 0.66883 0.503602 omega -0.262853 0.054392 -4.83257 0.000001 alpha1 -0.214794 0.025768 -8.33573 0.000000 beta1 0.971598 0.004886 198.83937 0.000000 gamma1 0.141275 0.007465 18.92445 0.000000 skew 0.832645 0.013163 63.25713 0.000000 shape 6.198791 1.636482 3.78788 0.000152 Information Criteria ------------------------------------ Akaike -6.4939, Bayes -6.4777 Shibata -6.4939, Hannan-Quinn -6.4880 Weighted Ljung-Box Test on Standardized Residuals ------------------------------------ statistic p-value Lag[1] 5.963 0.01461 Lag[2*(p+q)+(p+q)-1][2] 6.437 0.01683 Lag[4*(p+q)+(p+q)-1][5] 7.626 0.03647 d.o.f=0 H0 : No serial correlation Weighted Ljung-Box Test on Standardized Squared Residuals ------------------------------------ statistic p-value Lag[1] 4.246 0.03934 Lag[2*(p+q)+(p+q)-1][5] 5.234 0.13547 Lag[4*(p+q)+(p+q)-1][9] 6.351 0.25986 d.o.f=2 #### require(fds) data(sundaydemand) electric <- sundaydemand$y dim(electric) ## 48 by 508 ele=as.matrix(electric) test=matrix(ele,48*508,1) par(mfrow=c(1,1)) plot(test[1:(24*150)],type='l') plot(test,type='l') data(sundaytempairport) temp=sundaytempairport$y ### fill in the missing values. temp[1,90] <- (temp[1,89]+temp[1,91])/2 temp[12,6] <- (temp[12,5]+temp[12,7])/2 temp[14,269] <- (temp[14,268]+temp[14,270])/2 temp[15,79] <- (temp[15,78]+temp[15,80])/2 temp[17,11] <- (temp[17,10]+temp[17,12])/2 temp[17,46] <- (temp[17,45]+temp[17,47])/2 temp[22,37] <- (temp[22,36]+temp[22,38])/2 temp[28,92] <- (temp[28,91]+temp[28,93])/2 temp[29,92] <- (temp[29,91]+temp[29,93])/2 temp[30,92] <- (temp[30,91]+temp[30,93])/2 temp[33,9] <- (temp[33,8]+temp[33,10])/2 n.row=dim(temp)[1] n.col=dim(temp)[2] dep=matrix(electric,n.row*n.col,1) temp=as.matrix(temp) indep=matrix(temp,n.row*n.col,1) indep2=indep for(i in 1:(n.row*n.col)){ indep2[i]=max(0,indep[i]-25) } lm3=lm(dep~indep+indep2) summary(lm3) data=t(matrix(lm3$residuals,n.row,n.col)) t=dim(data)[1] mean=apply(data,2,mean) plot(mean,type='l') ############### ### Average removed data_re=data- matrix(rep(mean,each=t),t, dim(data)[2]) #pdf('plot.pdf') par(mfrow=c(1,3)) plot(data[1,],type='l',ylim=c(min(data)*1.1,max(data)*1.1), main='Original data',xlab='time',ylab='temperature',xaxt="n") axis(1,at=seq(0,48,by=6),labels=c("00:00","3:00","6:00","9:00","12:00", "15:00","18:00","21:00","24:00")) for(i in 2:508){ points(data[i,],type='l',col=grey(i/508)) } plot(mean,type='l',main='Temperature pattern', xlab='time',ylab='temperature',xaxt="n") axis(1,at=seq(0,48,by=6),labels=c("00:00","3:00","6:00","9:00","12:00", "15:00","18:00","21:00","24:00")) plot(data_re[1,],type='l',ylim=c(min(data_re)*1.1,max(data_re)*1.1), main='Average removed',xlab='time',ylab='temperature',xaxt="n") axis(1,at=seq(0,48,by=6),labels=c("00:00","3:00","6:00","9:00","12:00", "15:00","18:00","21:00","24:00")) for(i in 2:508){ points(data_re[i,],type='l',col=grey(i/508)) } #dev.off() par(mfrow=c(1,1)) #pdf('mean.pdf') plot(mean,type='l',main='Electricity pattern without temperature effect', xlab='time', ylab='Electricity demand',xaxt="n") axis(1,at=seq(0,48,by=6),labels=c("00:00","3:00","6:00","9:00","12:00", "15:00","18:00","21:00","24:00")) #dev.off() require(NTS) grid=470 df_b=5 p.max=8 F_test(data_re,p.max,df_b,grid) ###CFAR(4) #### Estimation p=4 est=est_cfar(data_re,p,6,grid) ### end of addition #pdf("phi.pdf") par(mfrow=c(2,2)) plot(seq(-1,1,by=1/470),est$phi_func[1,],type='l',main='phi1', xlab='x',ylab='',ylim=c(min(est$phi_func)*1.1,max(est$phi_func)*1.1)) plot(seq(-1,1,by=1/470),est$phi_func[2,],type='l',main='phi2', xlab='x',ylab='',ylim=c(min(est$phi_func)*1.1,max(est$phi_func)*1.1)) plot(seq(-1,1,by=1/470),est$phi_func[3,],type='l',main='phi3', xlab='x',ylab='',ylim=c(min(est$phi_func)*1.1,max(est$phi_func)*1.1)) plot(seq(-1,1,by=1/470),est$phi_func[4,],type='l',main='phi4', xlab='x',ylab='',ylim=c(min(est$phi_func)*1.1,max(est$phi_func)*1.1)) #dev.off() f=data_re n=dim(f)[2]-1 t=dim(f)[1] grid=470 index=seq(1,grid+1,by=grid/n) pred=p_cfar(est,f,1) pred_ele=pred[index]+mean #pdf("prediction.pdf") plot(pred_ele,type='l',main="Prediction of Electricity Demand", xlab='Time', ylab='Electricity demand',xaxt="n") axis(1,at=seq(0,48,by=6),labels=c("00:00","3:00","6:00","9:00","12:00", "15:00","18:00","21:00","24:00")) #dev.off() ####