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(1483509 + ,8036554 + ,4623093 + ,5528662 + ,4221032 + ,8061847 + ,7640066 + ,2935533 + ,8161548 + ,2543967 + ,13163450 + ,3348436 + ,3997440 + ,2322911 + ,2019457 + ,3047748 + ,5728767 + ,2605173 + ,5646743 + ,13121544 + ,3453409 + ,1878333 + ,4247362 + ,23022552 + ,7646203 + ,9016602 + ,3606568 + ,3173510 + ,17568772 + ,10805045 + ,31056269 + ,15623385 + ,6663443 + ,35435745 + ,2823250 + ,5197089 + ,4120632 + ,8832767 + ,3695374 + ,8385805 + ,3777904 + ,5199532 + ,5297275 + ,14847382 + ,5900158 + ,4416718 + ,3926429 + ,4876884 + ,2795297 + ,3385527 + ,3877941 + ,3556729 + ,4982836 + ,2976325 + ,2295026 + ,2218752 + ,4146062 + ,3302091 + ,3864505 + ,5454794 + ,1749836 + ,6684048 + ,2809918 + ,4092664 + ,5070470 + ,9814477 + ,6665318 + ,3912554 + ,6188129 + ,3627991 + ,3308767 + ,3820332 + ,4932979 + ,5567917 + ,5020814 + ,3803273 + ,3999984 + ,4883104 + ,13731747 + ,47531824 + ,8415570 + ,22178158 + ,61211654 + ,18223748 + ,17678085 + ,49299580 + ,25899948 + ,34121754 + ,9859231 + ,29740892 + ,21085212 + ,43003866 + ,59549247 + ,18026465 + ,4680597 + ,5564728 + ,11792347 + ,10371624 + ,3728446 + ,5732978 + ,4067638 + ,2395508 + ,5018801 + ,22068888 + ,7678580 + ,15510095 + ,6471239 + ,14349204 + ,35151574 + ,8210488 + ,5022664 + ,13996871 + ,12822431 + ,14011552 + ,20260980 + ,23718976 + ,45833049 + ,30688420 + ,16576062 + ,14844405 + ,16728286 + ,43477680 + ,57497427 + ,24233726 + ,24921208 + ,9516725 + ,27977239 + ,21632046 + ,22956809 + ,9704324 + ,19871149 + ,5553842 + ,5667858 + ,4348188 + ,10025042 + ,10639796 + ,8639184 + ,10764378 + ,12097733 + ,3988414 + ,4607102 + ,7126895 + ,6009625 + ,21533237 + ,5986771 + ,5455310 + ,1822874 + ,3374062 + ,2920748 + ,2295942 + ,6809829 + ,3318281 + ,13784645 + ,7366577 + ,1628637 + ,4258976 + ,7159779 + ,8098401 + ,6894240 + ,3771246 + ,3249726 + ,3147380 + ,4063037 + ,9621916 + ,5890158 + ,2142901 + ,3145007 + ,1562168 + ,3303103 + ,5886910 + ,3454270 + ,6995348 + ,6487869 + ,12091976 + ,3934625 + ,3999749 + ,3613526 + ,4271706 + ,4253390 + ,5551591 + ,4663041 + ,2104104 + ,5385399 + ,6205877 + ,7529500 + ,17222705 + ,6230913 + ,6508275 + ,4518884 + ,4234991 + ,5625388 + ,5810139 + ,6942187 + ,3711188 + ,4261281 + ,1989945 + ,5033342 + ,7239565 + ,11058795 + ,7384772 + ,3884771 + ,3239201 + ,2316403 + ,4034947 + ,3245271 + ,2387251 + ,2174886 + ,3436080 + ,3738956 + ,1884730 + ,1509144 + ,42728366 + ,3446317 + ,4600683 + ,2953615 + ,3570060 + ,2130208 + ,2442943 + ,4892020 + ,3222192 + ,3121617 + ,3665542 + ,5519432 + ,4113468 + ,1714614 + ,3651985 + ,2419548 + ,2378854 + ,2303949 + ,2555534 + ,1713005 + ,1705960 + ,6115046 + ,3951044 + ,3785568 + ,4670530 + ,2265100 + ,1105643 + ,2814152 + ,3728673 + ,2038949 + ,2402919 + ,2348814 + ,2797822 + ,902505 + ,1331319 + ,4204238 + ,2212485 + ,6797382 + ,4532324 + ,1778808 + ,1890720 + ,5463736 + ,11368931 + ,2040164 + ,4276399 + ,3714445 + ,2068168 + ,1003842 + ,2858535 + ,2355484 + ,2719262 + ,1897741 + ,3945185 + ,3799916 + ,1017654 + ,3052241 + ,3932970 + ,3598151 + ,2296005 + ,2202018 + ,2461777 + ,2452042 + ,2185142 + ,11968502 + ,20395972 + ,21756900 + ,30024300 + ,10811344 + ,1819202 + ,1276885 + ,2946701 + ,3587459 + ,2832691 + ,6674805 + ,3868362 + ,4302909 + ,23265229 + ,22348002 + ,11883953 + ,6634979 + ,2935493 + ,3425669 + ,1171611 + ,6875879 + ,19451908 + ,13885933 + ,7643317 + ,10797966 + ,7297445 + ,8739736 + ,12455537 + ,24291181 + ,4215150 + ,28652176 + ,6851172 + ,3746871 + ,7327861 + ,16829710 + ,13778594 + ,6463717 + ,8956867 + ,21204915 + ,16115855 + ,2536113 + ,16645717 + ,17003730 + ,15969006 + ,31020427 + ,23798897 + ,20770321 + ,44410402 + ,27037491 + ,29627771 + ,18189792 + ,4654610 + ,12307201 + ,15300578 + ,10623864 + ,6880178 + ,29947357 + ,18611399 + ,42432604 + ,20208278 + ,14004392 + ,25737765 + ,16735738 + ,22450825 + ,6880840 + ,8510379 + ,8182481 + ,10948683 + ,4805277 + ,2589229 + ,5658407 + ,12862611 + ,5666188 + ,6875556 + ,7098766 + ,36083309 + ,10200330 + ,7784976) > par10 = 'FALSE' > par9 = '0' > par8 = '0' > par7 = '0' > par6 = '0' > par5 = '12' > par4 = '0' > par3 = '1' > par2 = '0.0' > par1 = '12' > 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 <- as.numeric(par5) #seasonal period > par6 <- 6 #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") Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 -0.5627 -0.4290 -0.3488 -0.2772 -0.2491 -0.1765 s.e. 0.0538 0.0606 0.0631 0.0630 0.0603 0.0537 sigma^2 estimated as 0.4518: log likelihood = -346.7, aic = 707.4 > (forecast <- predict(arima.out,par1)) $pred Time Series: Start = 341 End = 352 Frequency = 1 [1] 16.27451 16.33538 16.42125 16.41985 16.43023 16.26151 16.25746 16.29676 [9] 16.31755 16.33483 16.34380 16.34398 $se Time Series: Start = 341 End = 352 Frequency = 1 [1] 0.6721321 0.7335776 0.7653902 0.7883161 0.8093743 0.8255986 0.8457691 [8] 0.8897820 0.9227970 0.9504768 0.9753370 0.9987302 > (lb <- forecast$pred - 1.96 * forecast$se) Time Series: Start = 341 End = 352 Frequency = 1 [1] 14.95713 14.89757 14.92108 14.87475 14.84386 14.64334 14.59976 14.55279 [9] 14.50887 14.47190 14.43214 14.38647 > (ub <- forecast$pred + 1.96 * forecast$se) Time Series: Start = 341 End = 352 Frequency = 1 [1] 17.59189 17.77319 17.92141 17.96495 18.01660 17.87969 17.91517 18.04074 [9] 18.12623 18.19776 18.25546 18.30149 > 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] 1483509 8036554 4623093 5528662 4221032 8061847 7640066 2935533 [9] 8161548 2543967 13163450 3348436 3997440 2322911 2019457 3047748 [17] 5728767 2605173 5646743 13121544 3453409 1878333 4247362 23022552 [25] 7646203 9016602 3606568 3173510 17568772 10805045 31056269 15623385 [33] 6663443 35435745 2823250 5197089 4120632 8832767 3695374 8385805 [41] 3777904 5199532 5297275 14847382 5900158 4416718 3926429 4876884 [49] 2795297 3385527 3877941 3556729 4982836 2976325 2295026 2218752 [57] 4146062 3302091 3864505 5454794 1749836 6684048 2809918 4092664 [65] 5070470 9814477 6665318 3912554 6188129 3627991 3308767 3820332 [73] 4932979 5567917 5020814 3803273 3999984 4883104 13731747 47531824 [81] 8415570 22178158 61211654 18223748 17678085 49299580 25899948 34121754 [89] 9859231 29740892 21085212 43003866 59549247 18026465 4680597 5564728 [97] 11792347 10371624 3728446 5732978 4067638 2395508 5018801 22068888 [105] 7678580 15510095 6471239 14349204 35151574 8210488 5022664 13996871 [113] 12822431 14011552 20260980 23718976 45833049 30688420 16576062 14844405 [121] 16728286 43477680 57497427 24233726 24921208 9516725 27977239 21632046 [129] 22956809 9704324 19871149 5553842 5667858 4348188 10025042 10639796 [137] 8639184 10764378 12097733 3988414 4607102 7126895 6009625 21533237 [145] 5986771 5455310 1822874 3374062 2920748 2295942 6809829 3318281 [153] 13784645 7366577 1628637 4258976 7159779 8098401 6894240 3771246 [161] 3249726 3147380 4063037 9621916 5890158 2142901 3145007 1562168 [169] 3303103 5886910 3454270 6995348 6487869 12091976 3934625 3999749 [177] 3613526 4271706 4253390 5551591 4663041 2104104 5385399 6205877 [185] 7529500 17222705 6230913 6508275 4518884 4234991 5625388 5810139 [193] 6942187 3711188 4261281 1989945 5033342 7239565 11058795 7384772 [201] 3884771 3239201 2316403 4034947 3245271 2387251 2174886 3436080 [209] 3738956 1884730 1509144 42728366 3446317 4600683 2953615 3570060 [217] 2130208 2442943 4892020 3222192 3121617 3665542 5519432 4113468 [225] 1714614 3651985 2419548 2378854 2303949 2555534 1713005 1705960 [233] 6115046 3951044 3785568 4670530 2265100 1105643 2814152 3728673 [241] 2038949 2402919 2348814 2797822 902505 1331319 4204238 2212485 [249] 6797382 4532324 1778808 1890720 5463736 11368931 2040164 4276399 [257] 3714445 2068168 1003842 2858535 2355484 2719262 1897741 3945185 [265] 3799916 1017654 3052241 3932970 3598151 2296005 2202018 2461777 [273] 2452042 2185142 11968502 20395972 21756900 30024300 10811344 1819202 [281] 1276885 2946701 3587459 2832691 6674805 3868362 4302909 23265229 [289] 22348002 11883953 6634979 2935493 3425669 1171611 6875879 19451908 [297] 13885933 7643317 10797966 7297445 8739736 12455537 24291181 4215150 [305] 28652176 6851172 3746871 7327861 16829710 13778594 6463717 8956867 [313] 21204915 16115855 2536113 16645717 17003730 15969006 31020427 23798897 [321] 20770321 44410402 27037491 29627771 18189792 4654610 12307201 15300578 [329] 10623864 6880178 29947357 18611399 42432604 20208278 14004392 25737765 [337] 16735738 22450825 6880840 8510379 11693132 12426978 13541201 13522342 [345] 13663369 11542122 11495464 11956247 12207338 12420143 12532090 12534264 > (perc.se <- (ub-forecast$pred)/1.96/forecast$pred) Time Series: Start = 341 End = 352 Frequency = 1 [1] 1.394705 1.638506 1.776749 1.881856 1.982652 2.063198 2.166973 2.408174 [9] 2.603265 2.776844 2.940976 3.102898 > postscript(file="/var/www/html/rcomp/tmp/1xajh1291065029.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/24byb1291065029.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/3tuvm1291065029.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/4wcca1291065029.tab") > > try(system("convert tmp/1xajh1291065029.ps tmp/1xajh1291065029.png",intern=TRUE)) character(0) > try(system("convert tmp/24byb1291065029.ps tmp/24byb1291065029.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.765 0.418 1.664