R version 2.12.0 (2010-10-15) Copyright (C) 2010 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) 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 + ,109.88 + ,109.93 + ,110.46 + ,109.56 + ,111.49 + ,111.85 + ,111.35 + ,110.95 + ,112.49 + ,113.11 + ,112.54 + ,112.84 + ,111.5 + ,111.52 + ,111.57 + ,112.48 + ,112.31 + ,113.79 + ,114.01 + ,113.64 + ,112.62 + ,113.27 + ,113.51 + ,112.92 + ,113.66 + ,113.14 + ,113.48 + ,113.23 + ,110.56 + ,109.5 + ,109.78 + ,109.49 + ,109.66 + ,109.93 + ,109.82 + ,108.54 + ,108.23 + ,106.19 + ,106.49 + ,107.15 + ,107.74 + ,107.54 + ,107.07 + ,107.54 + ,107.81 + ,108.38 + ,108.42 + ,106.86 + ,106.41 + ,106.46 + ,106.84 + ,107.69 + ,107.04 + ,111.04 + ,111.93 + ,111.98 + ,112.07 + ,112.05 + ,113.14 + ,112.49 + ,113.2 + ,113.52 + ,113.22 + ,113.85 + ,113.68 + ,114.26 + ,114.1 + ,114.8 + ,114.98 + ,115.1 + ,114.21 + ,114.24 + ,113.35 + ,114.23 + ,114.43 + ,114.28 + ,113 + ,113.16 + ,112.59 + ,113.65 + ,113.18 + ,113.21 + ,113.11 + ,112.78 + ,112.57 + ,111.87 + ,111.94 + ,113.18 + ,113.67 + ,115.15 + ,114.41 + ,112.88 + ,112.44 + ,113.48 + ,112.78 + ,112.59 + ,113.31 + ,113.21 + ,112.5 + ,113.72 + ,114.09 + ,113.97 + ,112.5 + ,111.28 + ,111.35 + ,110.92 + ,110.73 + ,109) > par10 = 'FALSE' > par9 = '0' > par8 = '0' > par7 = '0' > par6 = '0' > par5 = '1' > par4 = '0' > par3 = '1' > par2 = '1' > par1 = '15' > par1 <- as.numeric(par1) #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.438: log likelihood = -760.27, aic = 1522.54 > (forecast <- predict(arima.out,par1)) $pred Time Series: Start = 477 End = 491 Frequency = 1 [1] 112.44 112.44 112.44 112.44 112.44 112.44 112.44 112.44 112.44 112.44 [11] 112.44 112.44 112.44 112.44 112.44 $se Time Series: Start = 477 End = 491 Frequency = 1 [1] 1.199169 1.695881 2.077022 2.398338 2.681423 2.937352 3.172703 3.391762 [9] 3.597507 3.792105 3.977194 4.154043 4.323665 4.486880 4.644362 > (lb <- forecast$pred - 1.96 * forecast$se) Time Series: Start = 477 End = 491 Frequency = 1 [1] 110.0896 109.1161 108.3690 107.7393 107.1844 106.6828 106.2215 105.7921 [9] 105.3889 105.0075 104.6447 104.2981 103.9656 103.6457 103.3371 > (ub <- forecast$pred + 1.96 * forecast$se) Time Series: Start = 477 End = 491 Frequency = 1 [1] 114.7904 115.7639 116.5110 117.1407 117.6956 118.1972 118.6585 119.0879 [9] 119.4911 119.8725 120.2353 120.5819 120.9144 121.2343 121.5429 > 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 125.81 124.95 [341] 123.66 122.66 119.34 117.84 120.97 117.38 118.06 116.99 115.55 114.17 [351] 115.32 112.49 111.93 112.08 111.63 109.53 111.35 110.79 113.06 112.62 [361] 110.65 112.36 113.74 111.73 109.86 109.32 109.99 109.84 111.13 112.43 [371] 111.77 112.15 112.89 112.12 113.10 111.09 110.76 109.59 109.99 110.25 [381] 108.31 108.79 108.14 109.88 109.93 110.46 109.56 111.49 111.85 111.35 [391] 110.95 112.49 113.11 112.54 112.84 111.50 111.52 111.57 112.48 112.31 [401] 113.79 114.01 113.64 112.62 113.27 113.51 112.92 113.66 113.14 113.48 [411] 113.23 110.56 109.50 109.78 109.49 109.66 109.93 109.82 108.54 108.23 [421] 106.19 106.49 107.15 107.74 107.54 107.07 107.54 107.81 108.38 108.42 [431] 106.86 106.41 106.46 106.84 107.69 107.04 111.04 111.93 111.98 112.07 [441] 112.05 113.14 112.49 113.20 113.52 113.22 113.85 113.68 114.26 114.10 [451] 114.80 114.98 115.10 114.21 114.24 113.35 114.23 114.43 114.28 113.00 [461] 113.16 112.59 113.65 113.18 113.21 113.11 112.78 112.57 111.87 111.94 [471] 113.18 113.67 115.15 114.41 112.88 112.44 112.44 112.44 112.44 112.44 [481] 112.44 112.44 112.44 112.44 112.44 112.44 112.44 112.44 112.44 112.44 [491] 112.44 > (perc.se <- (ub-forecast$pred)/1.96/forecast$pred) Time Series: Start = 477 End = 491 Frequency = 1 [1] 0.01066497 0.01508254 0.01847227 0.02132994 0.02384759 0.02612373 [7] 0.02821685 0.03016508 0.03199490 0.03372559 0.03537170 0.03694453 [13] 0.03845309 0.03990466 0.04130524 > postscript(file="/var/www/rcomp/tmp/1agky1292756193.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/rcomp/tmp/2hhh91292756193.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/36id31292756193.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/rcomp/tmp/4s1cr1292756193.tab") > > try(system("convert tmp/1agky1292756193.ps tmp/1agky1292756193.png",intern=TRUE)) character(0) > try(system("convert tmp/2hhh91292756193.ps tmp/2hhh91292756193.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.82 0.36 1.16