R version 2.15.2 (2012-10-26) -- "Trick or Treat" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i686-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(121,66,72,29,141,87,95,83,68,105,102,66,124,59,76,38,136,90,95,85,64,118,104,73,145,78,90,33,140,113,114,102,69,129,124,114,135,70,84,22,109,105,107,86,67,124,118,107,128,65,75,20,109,100,100,84,71,128,109,102,142,88,90,31,128,116,112,93,58,129,129,125,130,75,77,27,162,89,101,64,57,128,105,80,131,62,60,28,147,87,100,81,69,125,100,95,141,85,92,28,148,111,111,100,76,125,125,120,140,82,88,25,103,110,107,96,74,130,116,117,142,83,83,21,102,104,105,93,77,125,112,99,140,78,69,24,100,85,104,102,81,122,97,64,132,81,73,28,117,96,106,78,77,129,107,82,132,75,78,33,139,99,105,92,64,124,114,97,151,91,92,31,122,117,114,99,67,144,130,121),dim=c(12,15),dimnames=list(c('Voedingsmiddelen','Tabaksproducten','Textiel','Kleding','Leer','Hout','Papier','Uitgeverijen','Cokes','Chemische','Rubber','Nietmetaalhoudende'),1:15)) > y <- array(NA,dim=c(12,15),dimnames=list(c('Voedingsmiddelen','Tabaksproducten','Textiel','Kleding','Leer','Hout','Papier','Uitgeverijen','Cokes','Chemische','Rubber','Nietmetaalhoudende'),1:15)) > 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 = '6' > par3 <- 'No Linear Trend' > par2 <- 'Do not include Seasonal Dummies' > par1 <- '6' > #'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) > library(lmtest) Loading required package: zoo Attaching package: 'zoo' The following object(s) are masked from 'package:base': as.Date, as.Date.numeric > 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 Hout Voedingsmiddelen Tabaksproducten Textiel Kleding Leer Papier 1 87 121 66 72 29 141 95 2 90 124 59 76 38 136 95 3 113 145 78 90 33 140 114 4 105 135 70 84 22 109 107 5 100 128 65 75 20 109 100 6 116 142 88 90 31 128 112 7 89 130 75 77 27 162 101 8 87 131 62 60 28 147 100 9 111 141 85 92 28 148 111 10 110 140 82 88 25 103 107 11 104 142 83 83 21 102 105 12 85 140 78 69 24 100 104 13 96 132 81 73 28 117 106 14 99 132 75 78 33 139 105 15 117 151 91 92 31 122 114 Uitgeverijen Cokes Chemische Rubber Nietmetaalhoudende 1 83 68 105 102 66 2 85 64 118 104 73 3 102 69 129 124 114 4 86 67 124 118 107 5 84 71 128 109 102 6 93 58 129 129 125 7 64 57 128 105 80 8 81 69 125 100 95 9 100 76 125 125 120 10 96 74 130 116 117 11 93 77 125 112 99 12 102 81 122 97 64 13 78 77 129 107 82 14 92 64 124 114 97 15 99 67 144 130 121 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) Voedingsmiddelen Tabaksproducten Textiel 7.805525 0.041937 0.005763 0.237616 Kleding Leer Papier Uitgeverijen 0.214847 -0.101546 0.010826 -0.096319 Cokes Chemische Rubber Nietmetaalhoudende 0.132943 -0.035538 0.508570 0.204355 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: 1 2 3 4 5 6 7 8 1.01459 -0.72391 2.04768 -1.74098 1.24207 0.15340 -0.04258 -0.28430 9 10 11 12 13 14 15 -1.38115 -0.00469 0.72043 -0.66760 0.25018 -0.45900 -0.12413 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.805525 23.804613 0.328 0.765 Voedingsmiddelen 0.041937 0.310967 0.135 0.901 Tabaksproducten 0.005763 0.140743 0.041 0.970 Textiel 0.237616 0.191176 1.243 0.302 Kleding 0.214847 0.225723 0.952 0.411 Leer -0.101546 0.050872 -1.996 0.140 Papier 0.010826 0.379701 0.029 0.979 Uitgeverijen -0.096319 0.165927 -0.580 0.602 Cokes 0.132943 0.162236 0.819 0.473 Chemische -0.035538 0.174826 -0.203 0.852 Rubber 0.508570 0.333614 1.524 0.225 Nietmetaalhoudende 0.204355 0.095757 2.134 0.123 Residual standard error: 2.128 on 3 degrees of freedom Multiple R-squared: 0.9923, Adjusted R-squared: 0.9642 F-statistic: 35.29 on 11 and 3 DF, p-value: 0.006808 > 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 + } > postscript(file="/var/fisher/rcomp/tmp/15zzn1353056517.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/fisher/rcomp/tmp/2rw3p1353056517.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/fisher/rcomp/tmp/37vog1353056517.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/fisher/rcomp/tmp/44sjs1353056517.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/fisher/rcomp/tmp/5mqt21353056517.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 = 15 Frequency = 1 1 2 3 4 5 6 1.014586596 -0.723914329 2.047679831 -1.740978967 1.242069320 0.153401054 7 8 9 10 11 12 -0.042576293 -0.284304005 -1.381154310 -0.004690297 0.720434109 -0.667604622 13 14 15 0.250183865 -0.458997962 -0.124133989 > postscript(file="/var/fisher/rcomp/tmp/6v0te1353056517.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 = 15 Frequency = 1 lag(myerror, k = 1) myerror 0 1.014586596 NA 1 -0.723914329 1.014586596 2 2.047679831 -0.723914329 3 -1.740978967 2.047679831 4 1.242069320 -1.740978967 5 0.153401054 1.242069320 6 -0.042576293 0.153401054 7 -0.284304005 -0.042576293 8 -1.381154310 -0.284304005 9 -0.004690297 -1.381154310 10 0.720434109 -0.004690297 11 -0.667604622 0.720434109 12 0.250183865 -0.667604622 13 -0.458997962 0.250183865 14 -0.124133989 -0.458997962 15 NA -0.124133989 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -0.723914329 1.014586596 [2,] 2.047679831 -0.723914329 [3,] -1.740978967 2.047679831 [4,] 1.242069320 -1.740978967 [5,] 0.153401054 1.242069320 [6,] -0.042576293 0.153401054 [7,] -0.284304005 -0.042576293 [8,] -1.381154310 -0.284304005 [9,] -0.004690297 -1.381154310 [10,] 0.720434109 -0.004690297 [11,] -0.667604622 0.720434109 [12,] 0.250183865 -0.667604622 [13,] -0.458997962 0.250183865 [14,] -0.124133989 -0.458997962 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -0.723914329 1.014586596 2 2.047679831 -0.723914329 3 -1.740978967 2.047679831 4 1.242069320 -1.740978967 5 0.153401054 1.242069320 6 -0.042576293 0.153401054 7 -0.284304005 -0.042576293 8 -1.381154310 -0.284304005 9 -0.004690297 -1.381154310 10 0.720434109 -0.004690297 11 -0.667604622 0.720434109 12 0.250183865 -0.667604622 13 -0.458997962 0.250183865 14 -0.124133989 -0.458997962 > 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/fisher/rcomp/tmp/78jqm1353056517.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/fisher/rcomp/tmp/8bn1a1353056517.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/fisher/rcomp/tmp/9tgkc1353056517.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') Warning messages: 1: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced 2: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced > par(opar) > dev.off() null device 1 > if (n > n25) { + postscript(file="/var/fisher/rcomp/tmp/10utfq1353056517.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() + } > > #Note: the /var/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/fisher/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/fisher/rcomp/tmp/11wssi1353056517.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/fisher/rcomp/tmp/12js6v1353056517.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/fisher/rcomp/tmp/1336mp1353056517.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/fisher/rcomp/tmp/144vg01353056517.tab") > if (n > n25) { + a<-table.start() + a<-table.row.start(a) + a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'p-values',header=TRUE) + a<-table.element(a,'Alternative Hypothesis',3,header=TRUE) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'breakpoint index',header=TRUE) + a<-table.element(a,'greater',header=TRUE) + a<-table.element(a,'2-sided',header=TRUE) + a<-table.element(a,'less',header=TRUE) + a<-table.row.end(a) + for (mypoint in kp3:nmkm3) { + a<-table.row.start(a) + a<-table.element(a,mypoint,header=TRUE) + a<-table.element(a,gqarr[mypoint-kp3+1,1]) + a<-table.element(a,gqarr[mypoint-kp3+1,2]) + a<-table.element(a,gqarr[mypoint-kp3+1,3]) + a<-table.row.end(a) + } + a<-table.end(a) + table.save(a,file="/var/fisher/rcomp/tmp/15tcn91353056517.tab") + a<-table.start() + a<-table.row.start(a) + a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'Description',header=TRUE) + a<-table.element(a,'# significant tests',header=TRUE) + a<-table.element(a,'% significant tests',header=TRUE) + a<-table.element(a,'OK/NOK',header=TRUE) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'1% type I error level',header=TRUE) + a<-table.element(a,numsignificant1) + a<-table.element(a,numsignificant1/numgqtests) + if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK' + a<-table.element(a,dum) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'5% type I error level',header=TRUE) + a<-table.element(a,numsignificant5) + a<-table.element(a,numsignificant5/numgqtests) + if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK' + a<-table.element(a,dum) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'10% type I error level',header=TRUE) + a<-table.element(a,numsignificant10) + a<-table.element(a,numsignificant10/numgqtests) + if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK' + a<-table.element(a,dum) + a<-table.row.end(a) + a<-table.end(a) + table.save(a,file="/var/fisher/rcomp/tmp/16vbip1353056517.tab") + } > > try(system("convert tmp/15zzn1353056517.ps tmp/15zzn1353056517.png",intern=TRUE)) character(0) > try(system("convert tmp/2rw3p1353056517.ps tmp/2rw3p1353056517.png",intern=TRUE)) character(0) > try(system("convert tmp/37vog1353056517.ps tmp/37vog1353056517.png",intern=TRUE)) character(0) > try(system("convert tmp/44sjs1353056517.ps tmp/44sjs1353056517.png",intern=TRUE)) character(0) > try(system("convert tmp/5mqt21353056517.ps tmp/5mqt21353056517.png",intern=TRUE)) character(0) > try(system("convert tmp/6v0te1353056517.ps tmp/6v0te1353056517.png",intern=TRUE)) character(0) > try(system("convert tmp/78jqm1353056517.ps tmp/78jqm1353056517.png",intern=TRUE)) character(0) > try(system("convert tmp/8bn1a1353056517.ps tmp/8bn1a1353056517.png",intern=TRUE)) character(0) > try(system("convert tmp/9tgkc1353056517.ps tmp/9tgkc1353056517.png",intern=TRUE)) character(0) > try(system("convert tmp/10utfq1353056517.ps tmp/10utfq1353056517.png",intern=TRUE)) convert: unable to open image `tmp/10utfq1353056517.ps': @ error/blob.c/OpenBlob/2587. convert: missing an image filename `tmp/10utfq1353056517.png' @ error/convert.c/ConvertImageCommand/3011. character(0) attr(,"status") [1] 1 Warning message: running command 'convert tmp/10utfq1353056517.ps tmp/10utfq1353056517.png' had status 1 > > > proc.time() user system elapsed 5.431 1.195 6.627