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 <- array(list(20.7246301,0,21.44580352,0,22.09413114,0,21.53321848,0,23.3470789,0,23.5656163,0,26.42117166,0,25.21193138,0,26.43574082,0,29.33500366,0,29.40056488,0,33.05013946,0,28.38072368,0,26.0059506,0,29.31314992,0,30.36212944,0,35.74543406,0,36.15337054,0,34.20838768,0,37.90895432,0,38.70297354,0,42.11944156,0,42.16314904,0,39.79566054,0,37.36261082,0,38.3533137,0,42.60022384,0,41.24529196,0,42.15586446,0,46.94183352,0,47.42990038,0,47.0583868,0,50.18347162,0,50.12519498,0,43.22669772,0,40.04333626,0,40.37114236,0,42.2141411,0,36.99838182,0,39.74466848,0,42.68035422,0,46.2935059,0,46.97097184,0,48.72655562,0,52.36884562,1,50.05234918,1,54.03701444,1,57.78128856,1,64.71620872,1,63.4122689,1,64.3592643,1,66.02743312,1,72.13919574,1,76.60464328,1,86.97060062,1,93.48301514,1,95.58825876,1,81.88596378,1,70.5511573,1,50.38015528,1,36.24807008,0),dim=c(2,61),dimnames=list(c('Olie','Dumivariabele'),1:61)) > y <- array(NA,dim=c(2,61),dimnames=list(c('Olie','Dumivariabele'),1:61)) > for (i in 1:dim(x)[1]) + { + for (j in 1:dim(x)[2]) + { + y[i,j] <- as.numeric(x[i,j]) + } + } > par3 = 'No Linear Trend' > par2 = 'Do not include Seasonal Dummies' > par1 = '1' > #'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!) > library(lattice) > par1 <- as.numeric(par1) > x <- t(y) > k <- length(x[1,]) > n <- length(x[,1]) > x1 <- cbind(x[,par1], x[,1:k!=par1]) > mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1]) > colnames(x1) <- mycolnames #colnames(x)[par1] > x <- x1 > if (par3 == 'First Differences'){ + x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep=''))) + for (i in 1:n-1) { + for (j in 1:k) { + x2[i,j] <- x[i+1,j] - x[i,j] + } + } + x <- x2 + } > if (par2 == 'Include Monthly Dummies'){ + x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep =''))) + for (i in 1:11){ + x2[seq(i,n,12),i] <- 1 + } + x <- cbind(x, x2) + } > if (par2 == 'Include Quarterly Dummies'){ + x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep =''))) + for (i in 1:3){ + x2[seq(i,n,4),i] <- 1 + } + x <- cbind(x, x2) + } > k <- length(x[1,]) > if (par3 == 'Linear Trend'){ + x <- cbind(x, c(1:n)) + colnames(x)[k+1] <- 't' + } > x Olie Dumivariabele 1 20.72463 0 2 21.44580 0 3 22.09413 0 4 21.53322 0 5 23.34708 0 6 23.56562 0 7 26.42117 0 8 25.21193 0 9 26.43574 0 10 29.33500 0 11 29.40056 0 12 33.05014 0 13 28.38072 0 14 26.00595 0 15 29.31315 0 16 30.36213 0 17 35.74543 0 18 36.15337 0 19 34.20839 0 20 37.90895 0 21 38.70297 0 22 42.11944 0 23 42.16315 0 24 39.79566 0 25 37.36261 0 26 38.35331 0 27 42.60022 0 28 41.24529 0 29 42.15586 0 30 46.94183 0 31 47.42990 0 32 47.05839 0 33 50.18347 0 34 50.12519 0 35 43.22670 0 36 40.04334 0 37 40.37114 0 38 42.21414 0 39 36.99838 0 40 39.74467 0 41 42.68035 0 42 46.29351 0 43 46.97097 0 44 48.72656 0 45 52.36885 1 46 50.05235 1 47 54.03701 1 48 57.78129 1 49 64.71621 1 50 63.41227 1 51 64.35926 1 52 66.02743 1 53 72.13920 1 54 76.60464 1 55 86.97060 1 56 93.48302 1 57 95.58826 1 58 81.88596 1 59 70.55116 1 60 50.38016 1 61 36.24807 0 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) Dumivariabele 36.23 32.54 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -18.720 -7.850 1.132 6.369 26.816 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 36.231 1.581 22.92 < 2e-16 *** Dumivariabele 32.541 3.087 10.54 3.43e-15 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 10.61 on 59 degrees of freedom Multiple R-squared: 0.6532, Adjusted R-squared: 0.6473 F-statistic: 111.1 on 1 and 59 DF, p-value: 3.430e-15 > postscript(file="/var/www/html/rcomp/tmp/104561229874652.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index') > points(x[,1]-mysum$resid) > grid() > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/2mrqa1229874652.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index') > grid() > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/3x4gq1229874652.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals') > grid() > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/4k3ep1229874652.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals') > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/5rgd81229874652.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > qqnorm(mysum$resid, main='Residual Normal Q-Q Plot') > grid() > dev.off() null device 1 > (myerror <- as.ts(mysum$resid)) Time Series: Start = 1 End = 61 Frequency = 1 1 2 3 4 5 6 -15.50644263 -14.78526921 -14.13694159 -14.69785425 -12.88399383 -12.66545643 7 8 9 10 11 12 -9.80990107 -11.01914135 -9.79533191 -6.89606907 -6.83050785 -3.18093327 13 14 15 16 17 18 -7.85034905 -10.22512213 -6.91792281 -5.86894329 -0.48563867 -0.07770219 19 20 21 22 23 24 -2.02268505 1.67788159 2.47190081 5.88836883 5.93207631 3.56458781 25 26 27 28 29 30 1.13153809 2.12224097 6.36915111 5.01421923 5.92479173 10.71076079 31 32 33 34 35 36 11.19882765 10.82731407 13.95239889 13.89412225 6.99562499 3.81226353 37 38 39 40 41 42 4.14006963 5.98306837 0.76730909 3.51359575 6.44928149 10.06243317 43 44 45 46 47 48 10.73989911 12.49548289 -16.40350830 -18.72000474 -14.73533948 -10.99106536 49 50 51 52 53 54 -4.05614520 -5.36008502 -4.41308962 -2.74492080 3.36684182 7.83228936 55 56 57 58 59 60 18.19824670 24.71066122 26.81590484 13.11360986 1.77880338 -18.39219864 61 0.01699735 > postscript(file="/var/www/html/rcomp/tmp/6lwzy1229874652.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > dum <- cbind(lag(myerror,k=1),myerror) > dum Time Series: Start = 0 End = 61 Frequency = 1 lag(myerror, k = 1) myerror 0 -15.50644263 NA 1 -14.78526921 -15.50644263 2 -14.13694159 -14.78526921 3 -14.69785425 -14.13694159 4 -12.88399383 -14.69785425 5 -12.66545643 -12.88399383 6 -9.80990107 -12.66545643 7 -11.01914135 -9.80990107 8 -9.79533191 -11.01914135 9 -6.89606907 -9.79533191 10 -6.83050785 -6.89606907 11 -3.18093327 -6.83050785 12 -7.85034905 -3.18093327 13 -10.22512213 -7.85034905 14 -6.91792281 -10.22512213 15 -5.86894329 -6.91792281 16 -0.48563867 -5.86894329 17 -0.07770219 -0.48563867 18 -2.02268505 -0.07770219 19 1.67788159 -2.02268505 20 2.47190081 1.67788159 21 5.88836883 2.47190081 22 5.93207631 5.88836883 23 3.56458781 5.93207631 24 1.13153809 3.56458781 25 2.12224097 1.13153809 26 6.36915111 2.12224097 27 5.01421923 6.36915111 28 5.92479173 5.01421923 29 10.71076079 5.92479173 30 11.19882765 10.71076079 31 10.82731407 11.19882765 32 13.95239889 10.82731407 33 13.89412225 13.95239889 34 6.99562499 13.89412225 35 3.81226353 6.99562499 36 4.14006963 3.81226353 37 5.98306837 4.14006963 38 0.76730909 5.98306837 39 3.51359575 0.76730909 40 6.44928149 3.51359575 41 10.06243317 6.44928149 42 10.73989911 10.06243317 43 12.49548289 10.73989911 44 -16.40350830 12.49548289 45 -18.72000474 -16.40350830 46 -14.73533948 -18.72000474 47 -10.99106536 -14.73533948 48 -4.05614520 -10.99106536 49 -5.36008502 -4.05614520 50 -4.41308962 -5.36008502 51 -2.74492080 -4.41308962 52 3.36684182 -2.74492080 53 7.83228936 3.36684182 54 18.19824670 7.83228936 55 24.71066122 18.19824670 56 26.81590484 24.71066122 57 13.11360986 26.81590484 58 1.77880338 13.11360986 59 -18.39219864 1.77880338 60 0.01699735 -18.39219864 61 NA 0.01699735 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -14.78526921 -15.50644263 [2,] -14.13694159 -14.78526921 [3,] -14.69785425 -14.13694159 [4,] -12.88399383 -14.69785425 [5,] -12.66545643 -12.88399383 [6,] -9.80990107 -12.66545643 [7,] -11.01914135 -9.80990107 [8,] -9.79533191 -11.01914135 [9,] -6.89606907 -9.79533191 [10,] -6.83050785 -6.89606907 [11,] -3.18093327 -6.83050785 [12,] -7.85034905 -3.18093327 [13,] -10.22512213 -7.85034905 [14,] -6.91792281 -10.22512213 [15,] -5.86894329 -6.91792281 [16,] -0.48563867 -5.86894329 [17,] -0.07770219 -0.48563867 [18,] -2.02268505 -0.07770219 [19,] 1.67788159 -2.02268505 [20,] 2.47190081 1.67788159 [21,] 5.88836883 2.47190081 [22,] 5.93207631 5.88836883 [23,] 3.56458781 5.93207631 [24,] 1.13153809 3.56458781 [25,] 2.12224097 1.13153809 [26,] 6.36915111 2.12224097 [27,] 5.01421923 6.36915111 [28,] 5.92479173 5.01421923 [29,] 10.71076079 5.92479173 [30,] 11.19882765 10.71076079 [31,] 10.82731407 11.19882765 [32,] 13.95239889 10.82731407 [33,] 13.89412225 13.95239889 [34,] 6.99562499 13.89412225 [35,] 3.81226353 6.99562499 [36,] 4.14006963 3.81226353 [37,] 5.98306837 4.14006963 [38,] 0.76730909 5.98306837 [39,] 3.51359575 0.76730909 [40,] 6.44928149 3.51359575 [41,] 10.06243317 6.44928149 [42,] 10.73989911 10.06243317 [43,] 12.49548289 10.73989911 [44,] -16.40350830 12.49548289 [45,] -18.72000474 -16.40350830 [46,] -14.73533948 -18.72000474 [47,] -10.99106536 -14.73533948 [48,] -4.05614520 -10.99106536 [49,] -5.36008502 -4.05614520 [50,] -4.41308962 -5.36008502 [51,] -2.74492080 -4.41308962 [52,] 3.36684182 -2.74492080 [53,] 7.83228936 3.36684182 [54,] 18.19824670 7.83228936 [55,] 24.71066122 18.19824670 [56,] 26.81590484 24.71066122 [57,] 13.11360986 26.81590484 [58,] 1.77880338 13.11360986 [59,] -18.39219864 1.77880338 [60,] 0.01699735 -18.39219864 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -14.78526921 -15.50644263 2 -14.13694159 -14.78526921 3 -14.69785425 -14.13694159 4 -12.88399383 -14.69785425 5 -12.66545643 -12.88399383 6 -9.80990107 -12.66545643 7 -11.01914135 -9.80990107 8 -9.79533191 -11.01914135 9 -6.89606907 -9.79533191 10 -6.83050785 -6.89606907 11 -3.18093327 -6.83050785 12 -7.85034905 -3.18093327 13 -10.22512213 -7.85034905 14 -6.91792281 -10.22512213 15 -5.86894329 -6.91792281 16 -0.48563867 -5.86894329 17 -0.07770219 -0.48563867 18 -2.02268505 -0.07770219 19 1.67788159 -2.02268505 20 2.47190081 1.67788159 21 5.88836883 2.47190081 22 5.93207631 5.88836883 23 3.56458781 5.93207631 24 1.13153809 3.56458781 25 2.12224097 1.13153809 26 6.36915111 2.12224097 27 5.01421923 6.36915111 28 5.92479173 5.01421923 29 10.71076079 5.92479173 30 11.19882765 10.71076079 31 10.82731407 11.19882765 32 13.95239889 10.82731407 33 13.89412225 13.95239889 34 6.99562499 13.89412225 35 3.81226353 6.99562499 36 4.14006963 3.81226353 37 5.98306837 4.14006963 38 0.76730909 5.98306837 39 3.51359575 0.76730909 40 6.44928149 3.51359575 41 10.06243317 6.44928149 42 10.73989911 10.06243317 43 12.49548289 10.73989911 44 -16.40350830 12.49548289 45 -18.72000474 -16.40350830 46 -14.73533948 -18.72000474 47 -10.99106536 -14.73533948 48 -4.05614520 -10.99106536 49 -5.36008502 -4.05614520 50 -4.41308962 -5.36008502 51 -2.74492080 -4.41308962 52 3.36684182 -2.74492080 53 7.83228936 3.36684182 54 18.19824670 7.83228936 55 24.71066122 18.19824670 56 26.81590484 24.71066122 57 13.11360986 26.81590484 58 1.77880338 13.11360986 59 -18.39219864 1.77880338 60 0.01699735 -18.39219864 > plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals') > lines(lowess(z)) > abline(lm(z)) > grid() > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/77eu61229874652.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function') > grid() > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/8yjcn1229874652.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function') > grid() > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/9xg6h1229874652.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0)) > plot(mylm, las = 1, sub='Residual Diagnostics') > par(opar) > 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, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE) > a<-table.row.end(a) > myeq <- colnames(x)[1] > myeq <- paste(myeq, '[t] = ', sep='') > for (i in 1:k){ + if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '') + myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ') + if (rownames(mysum$coefficients)[i] != '(Intercept)') { + myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='') + if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='') + } + } > myeq <- paste(myeq, ' + e[t]') > a<-table.row.start(a) > a<-table.element(a, myeq) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/10bsel1229874652.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Variable',header=TRUE) > a<-table.element(a,'Parameter',header=TRUE) > a<-table.element(a,'S.D.',header=TRUE) > a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE) > a<-table.element(a,'2-tail p-value',header=TRUE) > a<-table.element(a,'1-tail p-value',header=TRUE) > a<-table.row.end(a) > for (i in 1:k){ + a<-table.row.start(a) + a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE) + a<-table.element(a,mysum$coefficients[i,1]) + a<-table.element(a, round(mysum$coefficients[i,2],6)) + a<-table.element(a, round(mysum$coefficients[i,3],4)) + a<-table.element(a, round(mysum$coefficients[i,4],6)) + a<-table.element(a, round(mysum$coefficients[i,4]/2,6)) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/1146h71229874652.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'Multiple R',1,TRUE) > a<-table.element(a, sqrt(mysum$r.squared)) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'R-squared',1,TRUE) > a<-table.element(a, mysum$r.squared) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'Adjusted R-squared',1,TRUE) > a<-table.element(a, mysum$adj.r.squared) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'F-TEST (value)',1,TRUE) > a<-table.element(a, mysum$fstatistic[1]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE) > a<-table.element(a, mysum$fstatistic[2]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE) > a<-table.element(a, mysum$fstatistic[3]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'p-value',1,TRUE) > a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3])) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'Residual Standard Deviation',1,TRUE) > a<-table.element(a, mysum$sigma) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'Sum Squared Residuals',1,TRUE) > a<-table.element(a, sum(myerror*myerror)) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/12dvuj1229874652.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a, 'Time or Index', 1, TRUE) > a<-table.element(a, 'Actuals', 1, TRUE) > a<-table.element(a, 'Interpolation
Forecast', 1, TRUE) > a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE) > a<-table.row.end(a) > for (i in 1:n) { + a<-table.row.start(a) + a<-table.element(a,i, 1, TRUE) + a<-table.element(a,x[i]) + a<-table.element(a,x[i]-mysum$resid[i]) + a<-table.element(a,mysum$resid[i]) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/1393lj1229874652.tab") > > system("convert tmp/104561229874652.ps tmp/104561229874652.png") > system("convert tmp/2mrqa1229874652.ps tmp/2mrqa1229874652.png") > system("convert tmp/3x4gq1229874652.ps tmp/3x4gq1229874652.png") > system("convert tmp/4k3ep1229874652.ps tmp/4k3ep1229874652.png") > system("convert tmp/5rgd81229874652.ps tmp/5rgd81229874652.png") > system("convert tmp/6lwzy1229874652.ps tmp/6lwzy1229874652.png") > system("convert tmp/77eu61229874652.ps tmp/77eu61229874652.png") > system("convert tmp/8yjcn1229874652.ps tmp/8yjcn1229874652.png") > system("convert tmp/9xg6h1229874652.ps tmp/9xg6h1229874652.png") > > > proc.time() user system elapsed 1.941 1.433 3.168