R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > x <- c(126.64 + ,126.81 + ,125.84 + ,126.77 + ,124.34 + ,124.4 + ,120.48 + ,118.54 + ,117.66 + ,116.97 + ,120.11 + ,119.16 + ,116.9 + ,116.11 + ,114.98 + ,113.65 + ,115.82 + ,117.59 + ,118.57 + ,118.07 + ,114.98 + ,114.04 + ,115.02 + ,114.28 + ,115.04 + ,116.7 + ,119.21 + ,118.39 + ,116.5 + ,115.46 + ,117.59 + ,117.33 + ,116.2 + ,116.83 + ,118.99 + ,118.62 + ,121.09 + ,122.4 + ,123.76 + ,125.33 + ,123.23 + ,122.52 + ,123.64 + ,124.67 + ,124.71 + ,122.53 + ,124.4 + ,125.45 + ,125.35 + ,124.3 + ,127.03 + ,128.51 + ,128.1 + ,128.94 + ,129.67 + ,129.87 + ,131.12 + ,132.68 + ,132.24 + ,133.63 + ,129.91 + ,127.93 + ,131.17 + ,130.86 + ,133.48 + ,134.08 + ,136.02 + ,132.8 + ,132.37 + ,133.05 + ,132.57 + ,130.7 + ,130.5 + ,129.67 + ,127.8 + ,126.82 + ,126.85 + ,128.28 + ,128.3 + ,126.82 + ,125.08 + ,128.53 + ,130.34 + ,131.52 + ,132.59 + ,131.17 + ,132.72 + ,133.36 + ,132.82 + ,132.9 + ,130.9 + ,129.41 + ,128.67 + ,129.28 + ,130.91 + ,131.06 + ,130.84 + ,131.41 + ,133.22 + ,132.06 + ,132.48 + ,134.38 + ,135.22 + ,134.89 + ,136.09 + ,136.33 + ,136.32 + ,137.48 + ,136.53 + ,136.8 + ,138.03 + ,137.39 + ,137.55 + ,136.08 + ,134.78 + ,133.28 + ,133.57 + ,134.84 + ,133.02 + ,133.49 + ,133.77 + ,134.34 + ,134.5 + ,134.03 + ,135.51 + ,136.53 + ,135.95 + ,134.32 + ,132.44 + ,133.61 + ,131.02 + ,130.05 + ,128.21 + ,129.03 + ,130.34 + ,131.57 + ,132.63 + ,132.06 + ,134.44 + ,134.1 + ,132.49 + ,134.23 + ,134.92 + ,135.61 + ,134.53 + ,133.86 + ,133.89 + ,135.33 + ,135.86 + ,136.22 + ,137.38 + ,137.31 + ,136.89 + ,138.01 + ,136.72 + ,135.77 + ,137.52 + ,135.61 + ,132.94 + ,134.12 + ,132.55 + ,134.11 + ,134.19 + ,135.57 + ,135.05 + ,134.32 + ,133.61 + ,134.75 + ,133.1 + ,133.26 + ,131.63 + ,132.47 + ,132.45 + ,133.33 + ,133.57 + ,134.13 + ,133.92 + ,132.62 + ,132.3 + ,133.26 + ,132.6 + ,134.38 + ,134.17 + ,135.46 + ,135.09 + ,134.96 + ,133.85 + ,132.59 + ,131.15 + ,130.91 + ,131.07 + ,130.78 + ,129.95 + ,131.41 + ,131.21 + ,130.68 + ,130.46 + ,131.12 + ,132.99 + ,133.02 + ,133.39 + ,134.07 + ,135.6 + ,135.66 + ,135.53 + ,135.82 + ,136.9 + ,137.97 + ,138.09 + ,136.91 + ,134.76 + ,135.13 + ,134.66 + ,132.95 + ,132.25 + ,134.3 + ,134.3 + ,134.76 + ,134.81 + ,134.51 + ,135.11 + ,134.32 + ,133.51 + ,134.02 + ,132.76 + ,133.39 + ,132.05 + ,131.87 + ,133.03 + ,132.57 + ,132.1 + ,130.7 + ,129.2 + ,129.77 + ,131.02 + ,131.55 + ,133.17 + ,133.08 + ,133.24 + ,130.74 + ,129.91 + ,130.03 + ,131.13 + ,129.55 + ,130.22 + ,130.61 + ,129.27 + ,129.68 + ,130.1 + ,130.83 + ,130.95 + ,131.73 + ,131.86 + ,132.44 + ,132.35 + ,133.16 + ,133.62 + ,132.54 + ,132.69 + ,133.5 + ,133.36 + ,134.23 + ,132.41 + ,133.02 + ,132.88 + ,130.76 + ,130.33 + ,129.79 + ,128.65 + ,129.14 + ,127.35 + ,127.74 + ,126.31 + ,125.95 + ,126.36 + ,126.15 + ,125.6 + ,126.2 + ,126.73 + ,125.68 + ,122.49 + ,122.07 + ,123.4 + ,123.01 + ,123.03 + ,122.33 + ,122.42 + ,122.68 + ,124.69 + ,123.3 + ,124.17 + ,124.38 + ,123.19 + ,122.16 + ,120.66 + ,120.92 + ,120.67 + ,120.68 + ,121.1 + ,120.86 + ,121.48 + ,123.48 + ,121.72 + ,123.16 + ,123.84 + ,124.57 + ,124.3 + ,124.22 + ,124.43 + ,123.33 + ,122.86 + ,121.25 + ,122.16 + ,122.62 + ,123.44 + ,124 + ,124.75 + ,124.8 + ,125.93 + ,126.28 + ,126.04 + ,125.04 + ,123.76 + ,125.34 + ,126.99 + ,126.34 + ,127.42 + ,126.18 + ,125.3 + ,123.5 + ,125.32 + ,124.65 + ,124.03 + ,125.11 + ,125.46 + ,124.7 + ,124.48 + ,124.76 + ,125.81 + ,124.95 + ,123.66 + ,122.66 + ,119.34 + ,117.84 + ,120.97 + ,117.38 + ,118.06 + ,116.99 + ,115.55 + ,114.17 + ,115.32 + ,112.49 + ,111.93 + ,112.08 + ,111.63 + ,109.53 + ,111.35 + ,110.79 + ,113.06 + ,112.62 + ,110.65 + ,112.36 + ,113.74 + ,111.73 + ,109.86 + ,109.32 + ,109.99 + ,109.84 + ,111.13 + ,112.43 + ,111.77 + ,112.15 + ,112.89 + ,112.12 + ,113.1 + ,111.09 + ,110.76 + ,109.59 + ,109.99 + ,110.25 + ,108.31 + ,108.79 + ,108.14) > par10 = 'FALSE' > par9 = '0' > par8 = '0' > par7 = '0' > par6 = '0' > par5 = '1' > par4 = '0' > par3 = '1' > par2 = '1' > par1 = '15' > par1 <- 45 #cut off periods > par2 <- as.numeric(par2) #lambda > par3 <- as.numeric(par3) #degree of non-seasonal differencing > par4 <- as.numeric(par4) #degree of seasonal differencing > par5 <- 5 #seasonal period > par6 <- as.numeric(par6) #p > par7 <- as.numeric(par7) #q > par8 <- as.numeric(par8) #P > par9 <- as.numeric(par9) #Q > if (par10 == 'TRUE') par10 <- TRUE > if (par10 == 'FALSE') par10 <- FALSE > if (par2 == 0) x <- log(x) > if (par2 != 0) x <- x^par2 > lx <- length(x) > first <- lx - 2*par1 > nx <- lx - par1 > nx1 <- nx + 1 > fx <- lx - nx > if (fx < 1) { + fx <- par5 + nx1 <- lx + fx - 1 + first <- lx - 2*fx + } > first <- 1 > if (fx < 3) fx <- round(lx/10,0) > (arima.out <- arima(x[1:nx], order=c(par6,par3,par7), seasonal=list(order=c(par8,par4,par9), period=par5), include.mean=par10, method='ML')) Call: arima(x = x[1:nx], order = c(par6, par3, par7), seasonal = list(order = c(par8, par4, par9), period = par5), include.mean = par10, method = "ML") sigma^2 estimated as 1.505: log likelihood = -547.03, aic = 1096.06 > (forecast <- predict(arima.out,par1)) $pred Time Series: Start = 339 End = 383 Frequency = 1 [1] 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 [11] 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 [21] 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 [31] 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 [41] 124.76 124.76 124.76 124.76 124.76 $se Time Series: Start = 339 End = 383 Frequency = 1 [1] 1.226664 1.734765 2.124645 2.453329 2.742905 3.004702 3.245449 3.469531 [9] 3.679993 3.879054 4.068386 4.249290 4.422802 4.589758 4.750851 4.906658 [17] 5.057667 5.204296 5.346906 5.485810 5.621283 5.753566 5.882876 6.009404 [25] 6.133322 6.254786 6.373935 6.490898 6.605790 6.718718 6.829779 6.939062 [33] 7.046651 7.152621 7.257045 7.359987 7.461508 7.561667 7.660517 7.758107 [41] 7.854485 7.949694 8.043777 8.136771 8.228715 > (lb <- forecast$pred - 1.96 * forecast$se) Time Series: Start = 339 End = 383 Frequency = 1 [1] 122.3557 121.3599 120.5957 119.9515 119.3839 118.8708 118.3989 117.9597 [9] 117.5472 117.1571 116.7860 116.4314 116.0913 115.7641 115.4483 115.1430 [17] 114.8470 114.5596 114.2801 114.0078 113.7423 113.4830 113.2296 112.9816 [25] 112.7387 112.5006 112.2671 112.0378 111.8127 111.5913 111.3736 111.1594 [33] 110.9486 110.7409 110.5362 110.3344 110.1354 109.9391 109.7454 109.5541 [41] 109.3652 109.1786 108.9942 108.8119 108.6317 > (ub <- forecast$pred + 1.96 * forecast$se) Time Series: Start = 339 End = 383 Frequency = 1 [1] 127.1643 128.1601 128.9243 129.5685 130.1361 130.6492 131.1211 131.5603 [9] 131.9728 132.3629 132.7340 133.0886 133.4287 133.7559 134.0717 134.3770 [17] 134.6730 134.9604 135.2399 135.5122 135.7777 136.0370 136.2904 136.5384 [25] 136.7813 137.0194 137.2529 137.4822 137.7073 137.9287 138.1464 138.3606 [33] 138.5714 138.7791 138.9838 139.1856 139.3846 139.5809 139.7746 139.9659 [41] 140.1548 140.3414 140.5258 140.7081 140.8883 > if (par2 == 0) { + x <- exp(x) + forecast$pred <- exp(forecast$pred) + lb <- exp(lb) + ub <- exp(ub) + } > if (par2 != 0) { + x <- x^(1/par2) + forecast$pred <- forecast$pred^(1/par2) + lb <- lb^(1/par2) + ub <- ub^(1/par2) + } > if (par2 < 0) { + olb <- lb + lb <- ub + ub <- olb + } > (actandfor <- c(x[1:nx], forecast$pred)) [1] 126.64 126.81 125.84 126.77 124.34 124.40 120.48 118.54 117.66 116.97 [11] 120.11 119.16 116.90 116.11 114.98 113.65 115.82 117.59 118.57 118.07 [21] 114.98 114.04 115.02 114.28 115.04 116.70 119.21 118.39 116.50 115.46 [31] 117.59 117.33 116.20 116.83 118.99 118.62 121.09 122.40 123.76 125.33 [41] 123.23 122.52 123.64 124.67 124.71 122.53 124.40 125.45 125.35 124.30 [51] 127.03 128.51 128.10 128.94 129.67 129.87 131.12 132.68 132.24 133.63 [61] 129.91 127.93 131.17 130.86 133.48 134.08 136.02 132.80 132.37 133.05 [71] 132.57 130.70 130.50 129.67 127.80 126.82 126.85 128.28 128.30 126.82 [81] 125.08 128.53 130.34 131.52 132.59 131.17 132.72 133.36 132.82 132.90 [91] 130.90 129.41 128.67 129.28 130.91 131.06 130.84 131.41 133.22 132.06 [101] 132.48 134.38 135.22 134.89 136.09 136.33 136.32 137.48 136.53 136.80 [111] 138.03 137.39 137.55 136.08 134.78 133.28 133.57 134.84 133.02 133.49 [121] 133.77 134.34 134.50 134.03 135.51 136.53 135.95 134.32 132.44 133.61 [131] 131.02 130.05 128.21 129.03 130.34 131.57 132.63 132.06 134.44 134.10 [141] 132.49 134.23 134.92 135.61 134.53 133.86 133.89 135.33 135.86 136.22 [151] 137.38 137.31 136.89 138.01 136.72 135.77 137.52 135.61 132.94 134.12 [161] 132.55 134.11 134.19 135.57 135.05 134.32 133.61 134.75 133.10 133.26 [171] 131.63 132.47 132.45 133.33 133.57 134.13 133.92 132.62 132.30 133.26 [181] 132.60 134.38 134.17 135.46 135.09 134.96 133.85 132.59 131.15 130.91 [191] 131.07 130.78 129.95 131.41 131.21 130.68 130.46 131.12 132.99 133.02 [201] 133.39 134.07 135.60 135.66 135.53 135.82 136.90 137.97 138.09 136.91 [211] 134.76 135.13 134.66 132.95 132.25 134.30 134.30 134.76 134.81 134.51 [221] 135.11 134.32 133.51 134.02 132.76 133.39 132.05 131.87 133.03 132.57 [231] 132.10 130.70 129.20 129.77 131.02 131.55 133.17 133.08 133.24 130.74 [241] 129.91 130.03 131.13 129.55 130.22 130.61 129.27 129.68 130.10 130.83 [251] 130.95 131.73 131.86 132.44 132.35 133.16 133.62 132.54 132.69 133.50 [261] 133.36 134.23 132.41 133.02 132.88 130.76 130.33 129.79 128.65 129.14 [271] 127.35 127.74 126.31 125.95 126.36 126.15 125.60 126.20 126.73 125.68 [281] 122.49 122.07 123.40 123.01 123.03 122.33 122.42 122.68 124.69 123.30 [291] 124.17 124.38 123.19 122.16 120.66 120.92 120.67 120.68 121.10 120.86 [301] 121.48 123.48 121.72 123.16 123.84 124.57 124.30 124.22 124.43 123.33 [311] 122.86 121.25 122.16 122.62 123.44 124.00 124.75 124.80 125.93 126.28 [321] 126.04 125.04 123.76 125.34 126.99 126.34 127.42 126.18 125.30 123.50 [331] 125.32 124.65 124.03 125.11 125.46 124.70 124.48 124.76 124.76 124.76 [341] 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 [351] 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 [361] 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 [371] 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 124.76 [381] 124.76 124.76 124.76 > (perc.se <- (ub-forecast$pred)/1.96/forecast$pred) Time Series: Start = 339 End = 383 Frequency = 1 [1] 0.009832193 0.013904821 0.017029858 0.019664387 0.021985453 0.024083857 [7] 0.026013538 0.027809642 0.029496580 0.031092125 0.032609696 0.034059717 [13] 0.035450477 0.036788699 0.038079921 0.039328773 0.040539171 0.041714463 [19] 0.042857537 0.043970905 0.045056770 0.046117074 0.047153543 0.048167713 [25] 0.049160966 0.050134545 0.051089575 0.052027077 0.052947981 0.053853141 [31] 0.054743335 0.055619284 0.056481650 0.057331046 0.058168040 0.058993160 [37] 0.059806897 0.060609710 0.061402027 0.062184250 0.062956755 0.063719895 [43] 0.064474003 0.065219392 0.065956358 > postscript(file="/var/www/html/rcomp/tmp/1hdij1292756619.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > opar <- par(mar=c(4,4,2,2),las=1) > ylim <- c( min(x[first:nx],lb), max(x[first:nx],ub)) > plot(x,ylim=ylim,type='n',xlim=c(first,lx)) > usr <- par('usr') > rect(usr[1],usr[3],nx+1,usr[4],border=NA,col='lemonchiffon') > rect(nx1,usr[3],usr[2],usr[4],border=NA,col='lavender') > abline(h= (-3:3)*2 , col ='gray', lty =3) > polygon( c(nx1:lx,lx:nx1), c(lb,rev(ub)), col = 'orange', lty=2,border=NA) > lines(nx1:lx, lb , lty=2) > lines(nx1:lx, ub , lty=2) > lines(x, lwd=2) > lines(nx1:lx, forecast$pred , lwd=2 , col ='white') > box() > par(opar) > dev.off() null device 1 > prob.dec <- array(NA, dim=fx) > prob.sdec <- array(NA, dim=fx) > prob.ldec <- array(NA, dim=fx) > prob.pval <- array(NA, dim=fx) > perf.pe <- array(0, dim=fx) > perf.mape <- array(0, dim=fx) > perf.mape1 <- array(0, dim=fx) > perf.se <- array(0, dim=fx) > perf.mse <- array(0, dim=fx) > perf.mse1 <- array(0, dim=fx) > perf.rmse <- array(0, dim=fx) > for (i in 1:fx) { + locSD <- (ub[i] - forecast$pred[i]) / 1.96 + perf.pe[i] = (x[nx+i] - forecast$pred[i]) / forecast$pred[i] + perf.se[i] = (x[nx+i] - forecast$pred[i])^2 + prob.dec[i] = pnorm((x[nx+i-1] - forecast$pred[i]) / locSD) + prob.sdec[i] = pnorm((x[nx+i-par5] - forecast$pred[i]) / locSD) + prob.ldec[i] = pnorm((x[nx] - forecast$pred[i]) / locSD) + prob.pval[i] = pnorm(abs(x[nx+i] - forecast$pred[i]) / locSD) + } > perf.mape[1] = abs(perf.pe[1]) > perf.mse[1] = abs(perf.se[1]) > for (i in 2:fx) { + perf.mape[i] = perf.mape[i-1] + abs(perf.pe[i]) + perf.mape1[i] = perf.mape[i] / i + perf.mse[i] = perf.mse[i-1] + perf.se[i] + perf.mse1[i] = perf.mse[i] / i + } > perf.rmse = sqrt(perf.mse1) > postscript(file="/var/www/html/rcomp/tmp/2oexc1292756619.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(forecast$pred, pch=19, type='b',main='ARIMA Extrapolation Forecast', ylab='Forecast and 95% CI', xlab='time',ylim=c(min(lb),max(ub))) > dum <- forecast$pred > dum[1:par1] <- x[(nx+1):lx] > lines(dum, lty=1) > lines(ub,lty=3) > lines(lb,lty=3) > dev.off() null device 1 > > #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/rcomp/createtable") > > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Univariate ARIMA Extrapolation Forecast',9,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'time',1,header=TRUE) > a<-table.element(a,'Y[t]',1,header=TRUE) > a<-table.element(a,'F[t]',1,header=TRUE) > a<-table.element(a,'95% LB',1,header=TRUE) > a<-table.element(a,'95% UB',1,header=TRUE) > a<-table.element(a,'p-value
(H0: Y[t] = F[t])',1,header=TRUE) > a<-table.element(a,'P(F[t]>Y[t-1])',1,header=TRUE) > a<-table.element(a,'P(F[t]>Y[t-s])',1,header=TRUE) > mylab <- paste('P(F[t]>Y[',nx,sep='') > mylab <- paste(mylab,'])',sep='') > a<-table.element(a,mylab,1,header=TRUE) > a<-table.row.end(a) > for (i in (nx-par5):nx) { + a<-table.row.start(a) + a<-table.element(a,i,header=TRUE) + a<-table.element(a,x[i]) + a<-table.element(a,'-') + a<-table.element(a,'-') + a<-table.element(a,'-') + a<-table.element(a,'-') + a<-table.element(a,'-') + a<-table.element(a,'-') + a<-table.element(a,'-') + a<-table.row.end(a) + } > for (i in 1:fx) { + a<-table.row.start(a) + a<-table.element(a,nx+i,header=TRUE) + a<-table.element(a,round(x[nx+i],4)) + a<-table.element(a,round(forecast$pred[i],4)) + a<-table.element(a,round(lb[i],4)) + a<-table.element(a,round(ub[i],4)) + a<-table.element(a,round((1-prob.pval[i]),4)) + a<-table.element(a,round((1-prob.dec[i]),4)) + a<-table.element(a,round((1-prob.sdec[i]),4)) + a<-table.element(a,round((1-prob.ldec[i]),4)) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/3dfd61292756619.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Univariate ARIMA Extrapolation Forecast Performance',7,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'time',1,header=TRUE) > a<-table.element(a,'% S.E.',1,header=TRUE) > a<-table.element(a,'PE',1,header=TRUE) > a<-table.element(a,'MAPE',1,header=TRUE) > a<-table.element(a,'Sq.E',1,header=TRUE) > a<-table.element(a,'MSE',1,header=TRUE) > a<-table.element(a,'RMSE',1,header=TRUE) > a<-table.row.end(a) > for (i in 1:fx) { + a<-table.row.start(a) + a<-table.element(a,nx+i,header=TRUE) + a<-table.element(a,round(perc.se[i],4)) + a<-table.element(a,round(perf.pe[i],4)) + a<-table.element(a,round(perf.mape1[i],4)) + a<-table.element(a,round(perf.se[i],4)) + a<-table.element(a,round(perf.mse1[i],4)) + a<-table.element(a,round(perf.rmse[i],4)) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/4gfbu1292756619.tab") > > try(system("convert tmp/1hdij1292756619.ps tmp/1hdij1292756619.png",intern=TRUE)) character(0) > try(system("convert tmp/2oexc1292756619.ps tmp/2oexc1292756619.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.791 0.362 1.727