R version 2.8.0 (2008-10-20) Copyright (C) 2008 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(236.23 + ,239.63 + ,239.38 + ,240.44 + ,238.91 + ,235.21 + ,231.98 + ,226.43 + ,232.52 + ,235.84 + ,232.5 + ,238.61 + ,241.71 + ,242.42 + ,241.39 + ,239.78 + ,234.79 + ,237.41 + ,232.07 + ,230.78 + ,229.89 + ,233.2 + ,233.89 + ,230.92 + ,231.11 + ,231.66 + ,231.37 + ,232.38 + ,232.87 + ,233.42 + ,226.36 + ,224.56 + ,234.4 + ,235.94 + ,233.1 + ,232.92 + ,232.63 + ,234.98 + ,230.24 + ,231.83 + ,223.13 + ,227 + ,219.67 + ,217.99 + ,219.86 + ,216.47 + ,218.05 + ,223.38 + ,212.21 + ,216.49 + ,221.29 + ,225.45 + ,218.37 + ,216.42 + ,221.54 + ,216.41 + ,212.87 + ,217.03 + ,222.65 + ,230.4 + ,227.91 + ,228.72 + ,228.66 + ,228.27 + ,225.9 + ,216.84 + ,222.73 + ,225.56 + ,228.76 + ,226.32 + ,230.96 + ,226.52 + ,227.84 + ,220.77 + ,217.66 + ,224.68 + ,223.89 + ,226.87 + ,228.37 + ,225.27 + ,234.83 + ,238.08 + ,228.55 + ,228.71 + ,228.96 + ,227.75 + ,227.72 + ,229.29 + ,223.93 + ,226.79 + ,229.52 + ,226.39 + ,226.79 + ,230.86 + ,232.33 + ,234.26 + ,232.2 + ,228.83 + ,232.41 + ,228.21 + ,227.66 + ,229.35 + ,228.75 + ,218.74 + ,222.16 + ,226.21 + ,225.36 + ,223.28 + ,227.44 + ,226.12 + ,224.6 + ,225.23 + ,217.21 + ,216.98 + ,215.83 + ,219.92 + ,208.53 + ,206.45 + ,204.67 + ,202.97 + ,212.09 + ,210.09 + ,203.55 + ,204.38 + ,198.59 + ,199.47 + ,192.62 + ,193.4 + ,198.76 + ,191.77 + ,191.93 + ,203.56 + ,205 + ,205.9 + ,204.4 + ,199.64 + ,207.09 + ,212.96 + ,210.42 + ,210.88 + ,203.7 + ,204.15 + ,208.41 + ,211.19 + ,209.8 + ,212.54 + ,213.59 + ,217.29 + ,221.8 + ,213.19 + ,211.2 + ,211.32 + ,204.09 + ,202.08 + ,192.94 + ,198.55 + ,195.59 + ,199.02 + ,199.26 + ,196.71 + ,197.76 + ,190.02 + ,191.58 + ,187.2 + ,183.79 + ,185.24 + ,183.33 + ,197.2 + ,195.4 + ,189.02 + ,186.52 + ,193.54 + ,188.9 + ,192.89 + ,194.87 + ,192.97 + ,184.52 + ,182.49 + ,192.25 + ,192.31 + ,195.73 + ,195.14 + ,200.15 + ,195.93 + ,190.38 + ,190.98 + ,189.72 + ,192.69 + ,189.93 + ,188.39 + ,193.52 + ,192.69 + ,187.88 + ,182.19 + ,174.16 + ,173.66 + ,177.02 + ,178.39 + ,171.35 + ,165.72 + ,165.72 + ,165.72 + ,165.72 + ,167 + ,170.23 + ,158.67 + ,174.02 + ,168.34 + ,172.97 + ,175.05 + ,173.2 + ,167.52 + ,165.88 + ,161.84 + ,179.08 + ,175.53 + ,175.38 + ,172.44 + ,173.62 + ,171.85 + ,170.71 + ,173 + ,167.85 + ,163.6 + ,168.33 + ,172.2 + ,171.49 + ,178.36 + ,175.39 + ,169.05 + ,172.36 + ,177.3 + ,181.56 + ,177.02 + ,175.79 + ,175.76 + ,186.22 + ,182.15 + ,182.73 + ,180.77 + ,180.9 + ,181.17 + ,182.86 + ,183.93 + ,184.99 + ,187.57 + ,185.13 + ,184.88 + ,179.89 + ,172.42 + ,175.82 + ,171.72 + ,171.05 + ,176.08 + ,178.14 + ,183.39 + ,182.12 + ,177.1 + ,174.23 + ,176.24 + ,178.5 + ,168.19 + ,166.42 + ,174.26 + ,171.34 + ,172.6 + ,181.73 + ,181.08 + ,174.91 + ,171.93 + ,172.25 + ,167.29 + ,171.25 + ,169.23 + ,171.24 + ,160.77 + ,160.41 + ,160.89 + ,161.43 + ,154.26 + ,158.31 + ,159.76 + ,157.71 + ,160.66 + ,156.89 + ,157.18 + ,152.54 + ,150.51 + ,146.63 + ,155.64 + ,157.35 + ,158.43 + ,160.37 + ,163.1 + ,163.77 + ,157.35 + ,155.81 + ,153.36 + ,162.89 + ,168.29 + ,162.18 + ,161.04 + ,159.74 + ,170.55 + ,171.62 + ,168.08 + ,180.79 + ,181.64 + ,180.55 + ,180.9 + ,184.55 + ,189.13 + ,183.56 + ,181.04 + ,185.02 + ,182.67 + ,186.73 + ,181.84 + ,186.4 + ,187.64 + ,191.9 + ,188.16 + ,186.55 + ,189.47 + ,188.8 + ,181.66 + ,184.74 + ,195.04 + ,182.38 + ,187.96 + ,194.58 + ,191.64 + ,177.05 + ,186.46 + ,172.76 + ,178.95 + ,193.84 + ,176.37 + ,167.35 + ,168.37 + ,174.88 + ,175.77 + ,155.81 + ,175.03 + ,179.96 + ,171.59 + ,170.93 + ,165.22 + ,149.6 + ,149.63 + ,137.71 + ,134.98 + ,164.02 + ,154.66 + ,133.36 + ,143.2 + ,148.79 + ,162.56 + ,154.19 + ,145.09 + ,148.44 + ,146.24 + ,141.91 + ,165.92 + ,156.85 + ,165.01 + ,159.89 + ,167.07 + ,174.15 + ,151.02 + ,148.16 + ,153.91 + ,150.94 + ,151.97 + ,140.18 + ,153.15 + ,147.44 + ,140.43 + ,144.43 + ,130.63 + ,117.87 + ,133.24 + ,140.61 + ,136.35 + ,146.25 + ,146.97 + ,151.8 + ,131.72 + ,140.79 + ,145.98 + ,137.51 + ,146.06 + ,151.38 + ,144.46 + ,143.33) > par10 = 'FALSE' > par9 = '0' > par8 = '2' > par7 = '1' > par6 = '0' > par5 = '1' > par4 = '0' > par3 = '1' > par2 = '1' > par1 = '20' > #'GNU S' R Code compiled by R2WASP v. 1.0.44 () > #Author: Prof. Dr. P. Wessa > #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/ > #Source of accompanying publication: Office for Research, Development, and Education > #Technical description: Write here your technical program description (don't use hard returns!) > 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 <- 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") Coefficients: ma1 sar1 sar2 -0.3129 0.0314 -0.1974 s.e. 0.2025 0.1976 0.0751 sigma^2 estimated as 31.32: log likelihood = -1174.87, aic = 2357.75 > (forecast <- predict(arima.out,fx)) $pred Time Series: Start = 376 End = 395 Frequency = 1 [1] 143.3895 145.8176 145.2604 144.7636 144.8580 144.9590 144.9436 144.9231 [9] 144.9255 144.9297 144.9293 144.9285 144.9285 144.9287 144.9287 144.9287 [17] 144.9287 144.9287 144.9287 144.9287 $se Time Series: Start = 376 End = 395 Frequency = 1 [1] 5.596642 6.891601 7.464252 8.098530 8.774909 9.386056 9.943772 [8] 10.474250 10.982051 11.467038 11.931786 12.379147 12.810995 13.228749 [15] 13.633688 14.026943 14.409473 14.782108 15.145576 15.500524 > (lb <- forecast$pred - 1.96 * forecast$se) Time Series: Start = 376 End = 395 Frequency = 1 [1] 132.4201 132.3100 130.6305 128.8905 127.6592 126.5623 125.4538 124.3936 [9] 123.4007 122.4543 121.5430 120.6654 119.8190 119.0003 118.2067 117.4359 [17] 116.6861 115.9557 115.2433 114.5476 > (ub <- forecast$pred + 1.96 * forecast$se) Time Series: Start = 376 End = 395 Frequency = 1 [1] 154.3589 159.3251 159.8903 160.6368 162.0568 163.3557 164.4334 165.4527 [9] 166.4504 167.4051 168.3156 169.1916 170.0381 170.8570 171.6507 172.4215 [17] 173.1712 173.9016 174.6140 175.3097 > 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] 236.2300 239.6300 239.3800 240.4400 238.9100 235.2100 231.9800 226.4300 [9] 232.5200 235.8400 232.5000 238.6100 241.7100 242.4200 241.3900 239.7800 [17] 234.7900 237.4100 232.0700 230.7800 229.8900 233.2000 233.8900 230.9200 [25] 231.1100 231.6600 231.3700 232.3800 232.8700 233.4200 226.3600 224.5600 [33] 234.4000 235.9400 233.1000 232.9200 232.6300 234.9800 230.2400 231.8300 [41] 223.1300 227.0000 219.6700 217.9900 219.8600 216.4700 218.0500 223.3800 [49] 212.2100 216.4900 221.2900 225.4500 218.3700 216.4200 221.5400 216.4100 [57] 212.8700 217.0300 222.6500 230.4000 227.9100 228.7200 228.6600 228.2700 [65] 225.9000 216.8400 222.7300 225.5600 228.7600 226.3200 230.9600 226.5200 [73] 227.8400 220.7700 217.6600 224.6800 223.8900 226.8700 228.3700 225.2700 [81] 234.8300 238.0800 228.5500 228.7100 228.9600 227.7500 227.7200 229.2900 [89] 223.9300 226.7900 229.5200 226.3900 226.7900 230.8600 232.3300 234.2600 [97] 232.2000 228.8300 232.4100 228.2100 227.6600 229.3500 228.7500 218.7400 [105] 222.1600 226.2100 225.3600 223.2800 227.4400 226.1200 224.6000 225.2300 [113] 217.2100 216.9800 215.8300 219.9200 208.5300 206.4500 204.6700 202.9700 [121] 212.0900 210.0900 203.5500 204.3800 198.5900 199.4700 192.6200 193.4000 [129] 198.7600 191.7700 191.9300 203.5600 205.0000 205.9000 204.4000 199.6400 [137] 207.0900 212.9600 210.4200 210.8800 203.7000 204.1500 208.4100 211.1900 [145] 209.8000 212.5400 213.5900 217.2900 221.8000 213.1900 211.2000 211.3200 [153] 204.0900 202.0800 192.9400 198.5500 195.5900 199.0200 199.2600 196.7100 [161] 197.7600 190.0200 191.5800 187.2000 183.7900 185.2400 183.3300 197.2000 [169] 195.4000 189.0200 186.5200 193.5400 188.9000 192.8900 194.8700 192.9700 [177] 184.5200 182.4900 192.2500 192.3100 195.7300 195.1400 200.1500 195.9300 [185] 190.3800 190.9800 189.7200 192.6900 189.9300 188.3900 193.5200 192.6900 [193] 187.8800 182.1900 174.1600 173.6600 177.0200 178.3900 171.3500 165.7200 [201] 165.7200 165.7200 165.7200 167.0000 170.2300 158.6700 174.0200 168.3400 [209] 172.9700 175.0500 173.2000 167.5200 165.8800 161.8400 179.0800 175.5300 [217] 175.3800 172.4400 173.6200 171.8500 170.7100 173.0000 167.8500 163.6000 [225] 168.3300 172.2000 171.4900 178.3600 175.3900 169.0500 172.3600 177.3000 [233] 181.5600 177.0200 175.7900 175.7600 186.2200 182.1500 182.7300 180.7700 [241] 180.9000 181.1700 182.8600 183.9300 184.9900 187.5700 185.1300 184.8800 [249] 179.8900 172.4200 175.8200 171.7200 171.0500 176.0800 178.1400 183.3900 [257] 182.1200 177.1000 174.2300 176.2400 178.5000 168.1900 166.4200 174.2600 [265] 171.3400 172.6000 181.7300 181.0800 174.9100 171.9300 172.2500 167.2900 [273] 171.2500 169.2300 171.2400 160.7700 160.4100 160.8900 161.4300 154.2600 [281] 158.3100 159.7600 157.7100 160.6600 156.8900 157.1800 152.5400 150.5100 [289] 146.6300 155.6400 157.3500 158.4300 160.3700 163.1000 163.7700 157.3500 [297] 155.8100 153.3600 162.8900 168.2900 162.1800 161.0400 159.7400 170.5500 [305] 171.6200 168.0800 180.7900 181.6400 180.5500 180.9000 184.5500 189.1300 [313] 183.5600 181.0400 185.0200 182.6700 186.7300 181.8400 186.4000 187.6400 [321] 191.9000 188.1600 186.5500 189.4700 188.8000 181.6600 184.7400 195.0400 [329] 182.3800 187.9600 194.5800 191.6400 177.0500 186.4600 172.7600 178.9500 [337] 193.8400 176.3700 167.3500 168.3700 174.8800 175.7700 155.8100 175.0300 [345] 179.9600 171.5900 170.9300 165.2200 149.6000 149.6300 137.7100 134.9800 [353] 164.0200 154.6600 133.3600 143.2000 148.7900 162.5600 154.1900 145.0900 [361] 148.4400 146.2400 141.9100 165.9200 156.8500 165.0100 159.8900 167.0700 [369] 174.1500 151.0200 148.1600 153.9100 150.9400 151.9700 140.1800 143.3895 [377] 145.8176 145.2604 144.7636 144.8580 144.9590 144.9436 144.9231 144.9255 [385] 144.9297 144.9293 144.9285 144.9285 144.9287 144.9287 144.9287 144.9287 [393] 144.9287 144.9287 144.9287 > (perc.se <- (ub-forecast$pred)/1.96/forecast$pred) Time Series: Start = 376 End = 395 Frequency = 1 [1] 0.03903104 0.04726181 0.05138531 0.05594312 0.06057594 0.06474972 [7] 0.06860444 0.07227451 0.07577719 0.07912140 0.08232832 0.08541555 [13] 0.08839526 0.09127764 0.09407170 0.09678516 0.09942460 0.10199575 [19] 0.10450366 0.10695278 > postscript(file="/var/www/html/freestat/rcomp/tmp/10aok1229552441.ps",horizontal=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.se <- array(0, dim=fx) > perf.mse <- 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.mape[i] = perf.mape[i] + abs(perf.pe[i]) + perf.se[i] = (x[nx+i] - forecast$pred[i])^2 + perf.mse[i] = perf.mse[i] + perf.se[i] + 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 = perf.mape / fx > perf.mse = perf.mse / fx > perf.rmse = sqrt(perf.mse) > postscript(file="/var/www/html/freestat/rcomp/tmp/2mw6a1229552441.ps",horizontal=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:12] <- x[(nx+1):lx] Warning message: In NextMethod("[<-") : number of items to replace is not a multiple of replacement length > lines(dum, lty=1) > lines(ub,lty=3) > lines(lb,lty=3) > dev.off() null device 1 > > #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/3jce61229552441.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.mape[i],4)) + a<-table.element(a,round(perf.se[i],4)) + a<-table.element(a,round(perf.mse[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/freestat/rcomp/tmp/4k9sf1229552441.tab") > > system("convert tmp/10aok1229552441.ps tmp/10aok1229552441.png") > system("convert tmp/2mw6a1229552441.ps tmp/2mw6a1229552441.png") > > > proc.time() user system elapsed 0.916 0.460 1.011