R version 2.12.1 (2010-12-16) 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 <- array(list(33907 + ,71433 + ,152 + ,74272 + ,99 + ,765 + ,35981 + ,53655 + ,99 + ,78867 + ,128 + ,1371 + ,36588 + ,70556 + ,92 + ,80176 + ,57 + ,1880 + ,16967 + ,74702 + ,138 + ,36541 + ,95 + ,232 + ,25333 + ,61201 + ,106 + ,55107 + ,205 + ,230 + ,21027 + ,686 + ,95 + ,45527 + ,51 + ,828 + ,21114 + ,87586 + ,145 + ,46001 + ,59 + ,1833 + ,28777 + ,6615 + ,181 + ,62854 + ,194 + ,906 + ,35612 + ,89725 + ,190 + ,78112 + ,27 + ,1781 + ,24183 + ,40420 + ,150 + ,52653 + ,9 + ,1264 + ,22262 + ,49569 + ,186 + ,48467 + ,24 + ,1123 + ,20637 + ,13963 + ,174 + ,44873 + ,189 + ,1461 + ,29948 + ,62508 + ,151 + ,65605 + ,37 + ,820 + ,22093 + ,90901 + ,112 + ,48016 + ,81 + ,107 + ,36997 + ,89418 + ,143 + ,81110 + ,72 + ,1349 + ,31089 + ,83237 + ,120 + ,68019 + ,81 + ,870 + ,19477 + ,22183 + ,169 + ,42198 + ,90 + ,1471 + ,31301 + ,24346 + ,135 + ,68531 + ,216 + ,731 + ,18497 + ,74341 + ,161 + ,40071 + ,216 + ,1945 + ,30142 + ,24188 + ,98 + ,65849 + ,13 + ,521 + ,21326 + ,11781 + ,142 + ,46362 + ,153 + ,1920 + ,16779 + ,23072 + ,190 + ,36313 + ,185 + ,1924 + ,38068 + ,49119 + ,169 + ,83521 + ,131 + ,100 + ,29707 + ,67776 + ,130 + ,64932 + ,136 + ,34 + ,35016 + ,86910 + ,160 + ,76730 + ,182 + ,325 + ,26131 + ,69358 + ,176 + ,56982 + ,139 + ,1677 + ,29251 + ,16144 + ,111 + ,63793 + ,42 + ,1779 + ,22855 + ,77863 + ,165 + ,49740 + ,213 + ,477 + ,31806 + ,89070 + ,117 + ,69447 + ,184 + ,1007 + ,34124 + ,34790 + ,122 + ,74708 + ,44 + ,1527) + ,dim=c(6 + ,30) + ,dimnames=list(c('Omzet' + ,'Uitgaven' + ,'Prijs' + ,'Gemiddeld_Budget' + ,'Consumentevertrouwen' + ,'Uitgaven') + ,1:30)) > y <- array(NA,dim=c(6,30),dimnames=list(c('Omzet','Uitgaven','Prijs','Gemiddeld_Budget','Consumentevertrouwen','Uitgaven'),1:30)) > 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' > library(lattice) > library(lmtest) Loading required package: zoo > n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test > 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 Omzet Uitgaven Prijs Gemiddeld_Budget Consumentevertrouwen Uitgaven 1 33907 71433 152 74272 99 765 2 35981 53655 99 78867 128 1371 3 36588 70556 92 80176 57 1880 4 16967 74702 138 36541 95 232 5 25333 61201 106 55107 205 230 6 21027 686 95 45527 51 828 7 21114 87586 145 46001 59 1833 8 28777 6615 181 62854 194 906 9 35612 89725 190 78112 27 1781 10 24183 40420 150 52653 9 1264 11 22262 49569 186 48467 24 1123 12 20637 13963 174 44873 189 1461 13 29948 62508 151 65605 37 820 14 22093 90901 112 48016 81 107 15 36997 89418 143 81110 72 1349 16 31089 83237 120 68019 81 870 17 19477 22183 169 42198 90 1471 18 31301 24346 135 68531 216 731 19 18497 74341 161 40071 216 1945 20 30142 24188 98 65849 13 521 21 21326 11781 142 46362 153 1920 22 16779 23072 190 36313 185 1924 23 38068 49119 169 83521 131 100 24 29707 67776 130 64932 136 34 25 35016 86910 160 76730 182 325 26 26131 69358 176 56982 139 1677 27 29251 16144 111 63793 42 1779 28 22855 77863 165 49740 213 477 29 31806 89070 117 69447 184 1007 30 34124 34790 122 74708 44 1527 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) Uitgaven Prijs 5.649e+02 -1.152e-04 -6.698e-01 Gemiddeld_Budget Consumentevertrouwen 4.503e-01 5.220e-02 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -60.101 -15.438 2.706 19.936 49.517 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.649e+02 3.847e+01 14.682 8.5e-14 *** Uitgaven -1.152e-04 1.867e-04 -0.617 0.54283 Prijs -6.698e-01 1.832e-01 -3.657 0.00119 ** Gemiddeld_Budget 4.503e-01 3.868e-04 1164.068 < 2e-16 *** Consumentevertrouwen 5.220e-02 7.946e-02 0.657 0.51719 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 28.35 on 25 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 3.94e+05 on 4 and 25 DF, p-value: < 2.2e-16 > if (n > n25) { + kp3 <- k + 3 + nmkm3 <- n - k - 3 + gqarr <- array(NA, dim=c(nmkm3-kp3+1,3)) + numgqtests <- 0 + numsignificant1 <- 0 + numsignificant5 <- 0 + numsignificant10 <- 0 + for (mypoint in kp3:nmkm3) { + j <- 0 + numgqtests <- numgqtests + 1 + for (myalt in c('greater', 'two.sided', 'less')) { + j <- j + 1 + gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value + } + if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1 + if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1 + if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1 + } + gqarr + } [,1] [,2] [,3] [1,] 0.8943539 0.2112922 0.10564612 [2,] 0.8168103 0.3663793 0.18318966 [3,] 0.7181211 0.5637577 0.28187886 [4,] 0.7522842 0.4954317 0.24771584 [5,] 0.9086223 0.1827553 0.09137765 [6,] 0.8730556 0.2538889 0.12694445 [7,] 0.8275090 0.3449820 0.17249102 [8,] 0.8988706 0.2022588 0.10112940 [9,] 0.8666983 0.2666034 0.13330171 [10,] 0.8429146 0.3141707 0.15708535 [11,] 0.7509154 0.4981692 0.24908462 [12,] 0.7021376 0.5957248 0.29786238 [13,] 0.6184380 0.7631240 0.38156202 > postscript(file="/var/www/rcomp/tmp/17dow1323435270.ps",horizontal=F,onefile=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/rcomp/tmp/2x9gr1323435270.ps",horizontal=F,onefile=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/rcomp/tmp/3yziq1323435270.ps",horizontal=F,onefile=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/rcomp/tmp/4hrv11323435270.ps",horizontal=F,onefile=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/rcomp/tmp/54u811323435270.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > qqnorm(mysum$resid, main='Residual Normal Q-Q Plot') > qqline(mysum$resid) > grid() > dev.off() null device 1 > (myerror <- as.ts(mysum$resid)) Time Series: Start = 1 End = 30 Frequency = 1 1 2 3 4 5 6 7 3.762231 -30.337666 -11.789675 44.498369 21.861975 23.246690 -60.101430 8 9 10 11 12 13 14 22.060319 11.004229 14.158501 2.416020 -25.029691 -51.119563 -11.286314 15 16 17 18 19 20 21 12.198575 -17.768770 22.236283 -40.097688 -5.943819 -5.649726 -26.302948 22 23 24 25 26 27 28 -16.654740 7.264523 -7.707752 8.781999 26.857615 35.393460 1.537072 29 30 49.516756 2.995167 > postscript(file="/var/www/rcomp/tmp/6gx621323435270.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > dum <- cbind(lag(myerror,k=1),myerror) > dum Time Series: Start = 0 End = 30 Frequency = 1 lag(myerror, k = 1) myerror 0 3.762231 NA 1 -30.337666 3.762231 2 -11.789675 -30.337666 3 44.498369 -11.789675 4 21.861975 44.498369 5 23.246690 21.861975 6 -60.101430 23.246690 7 22.060319 -60.101430 8 11.004229 22.060319 9 14.158501 11.004229 10 2.416020 14.158501 11 -25.029691 2.416020 12 -51.119563 -25.029691 13 -11.286314 -51.119563 14 12.198575 -11.286314 15 -17.768770 12.198575 16 22.236283 -17.768770 17 -40.097688 22.236283 18 -5.943819 -40.097688 19 -5.649726 -5.943819 20 -26.302948 -5.649726 21 -16.654740 -26.302948 22 7.264523 -16.654740 23 -7.707752 7.264523 24 8.781999 -7.707752 25 26.857615 8.781999 26 35.393460 26.857615 27 1.537072 35.393460 28 49.516756 1.537072 29 2.995167 49.516756 30 NA 2.995167 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -30.337666 3.762231 [2,] -11.789675 -30.337666 [3,] 44.498369 -11.789675 [4,] 21.861975 44.498369 [5,] 23.246690 21.861975 [6,] -60.101430 23.246690 [7,] 22.060319 -60.101430 [8,] 11.004229 22.060319 [9,] 14.158501 11.004229 [10,] 2.416020 14.158501 [11,] -25.029691 2.416020 [12,] -51.119563 -25.029691 [13,] -11.286314 -51.119563 [14,] 12.198575 -11.286314 [15,] -17.768770 12.198575 [16,] 22.236283 -17.768770 [17,] -40.097688 22.236283 [18,] -5.943819 -40.097688 [19,] -5.649726 -5.943819 [20,] -26.302948 -5.649726 [21,] -16.654740 -26.302948 [22,] 7.264523 -16.654740 [23,] -7.707752 7.264523 [24,] 8.781999 -7.707752 [25,] 26.857615 8.781999 [26,] 35.393460 26.857615 [27,] 1.537072 35.393460 [28,] 49.516756 1.537072 [29,] 2.995167 49.516756 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -30.337666 3.762231 2 -11.789675 -30.337666 3 44.498369 -11.789675 4 21.861975 44.498369 5 23.246690 21.861975 6 -60.101430 23.246690 7 22.060319 -60.101430 8 11.004229 22.060319 9 14.158501 11.004229 10 2.416020 14.158501 11 -25.029691 2.416020 12 -51.119563 -25.029691 13 -11.286314 -51.119563 14 12.198575 -11.286314 15 -17.768770 12.198575 16 22.236283 -17.768770 17 -40.097688 22.236283 18 -5.943819 -40.097688 19 -5.649726 -5.943819 20 -26.302948 -5.649726 21 -16.654740 -26.302948 22 7.264523 -16.654740 23 -7.707752 7.264523 24 8.781999 -7.707752 25 26.857615 8.781999 26 35.393460 26.857615 27 1.537072 35.393460 28 49.516756 1.537072 29 2.995167 49.516756 > 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/rcomp/tmp/77e201323435270.ps",horizontal=F,onefile=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/rcomp/tmp/8urv11323435270.ps",horizontal=F,onefile=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/rcomp/tmp/94kx11323435270.ps",horizontal=F,onefile=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 > if (n > n25) { + postscript(file="/var/www/rcomp/tmp/10uess1323435270.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) + plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint') + grid() + 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, '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='') + } + } Error: subscript out of bounds Execution halted