R version 2.13.0 (2011-04-13) Copyright (C) 2011 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(229.00,48856.87,1000300.00,38894.00,53029.88,258291900.00,7544.00,40497.17,36129100.00,2592.00,36596.73,7882600.00,160.00,37061.55,385300.00,2073.00,39840.23,3931000.00,6483.00,50608.69,12657200.00,41.00,119660.32,48774000.00,829.00,76893.41,39529600.00,782.00,52288.86,9918700.00,1524.00,50764.75,9191000.00,7104.00,52145.04,37266300.00,1994.00,52509.60,13025100.00,2670.00,57360.74,11164400.00,973.00,61950.96,21304600.00,4125.00,36467.10,7133000.00,176.00,94484.11,36607100.00,66619.00,38666.46,43485600.00,136486.00,40168.91,356407900.00,18947.00,43649.10,81571400.00,43871.00,54290.36,203647600.00,73668.00,27303.67,71188900.00,44392.00,17743.22,10419300.00,19688.00,48363.52,65421400.00,118074.00,46018.40,79003600.00),dim=c(3,25),dimnames=list(c('#Ondernemingen','Personeelskost','Omzet'),1:25)) > y <- array(NA,dim=c(3,25),dimnames=list(c('#Ondernemingen','Personeelskost','Omzet'),1:25)) > 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) > 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 #Ondernemingen Personeelskost Omzet 1 229 48856.87 1000300 2 38894 53029.88 258291900 3 7544 40497.17 36129100 4 2592 36596.73 7882600 5 160 37061.55 385300 6 2073 39840.23 3931000 7 6483 50608.69 12657200 8 41 119660.32 48774000 9 829 76893.41 39529600 10 782 52288.86 9918700 11 1524 50764.75 9191000 12 7104 52145.04 37266300 13 1994 52509.60 13025100 14 2670 57360.74 11164400 15 973 61950.96 21304600 16 4125 36467.10 7133000 17 176 94484.11 36607100 18 66619 38666.46 43485600 19 136486 40168.91 356407900 20 18947 43649.10 81571400 21 43871 54290.36 203647600 22 73668 27303.67 71188900 23 44392 17743.22 10419300 24 19688 48363.52 65421400 25 118074 46018.40 79003600 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) Personeelskost Omzet 3.658e+04 -5.825e-01 2.930e-04 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -42481 -13305 -7892 7904 85148 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.658e+04 1.459e+04 2.508 0.020 * Personeelskost -5.825e-01 2.561e-01 -2.274 0.033 * Omzet 2.930e-04 6.089e-05 4.812 8.31e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 26030 on 22 degrees of freedom Multiple R-squared: 0.5666, Adjusted R-squared: 0.5272 F-statistic: 14.38 on 2 and 22 DF, p-value: 0.0001012 > 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/wessaorg/rcomp/tmp/1tg3g1321714524.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/wessaorg/rcomp/tmp/22cjz1321714524.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/wessaorg/rcomp/tmp/3cvd11321714524.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/wessaorg/rcomp/tmp/4r3z21321714524.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/wessaorg/rcomp/tmp/5q1xy1321714524.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 = 25 Frequency = 1 1 2 3 4 5 6 7 -8187.202 -42481.117 -16034.775 -14982.147 -14946.595 -12453.990 -4328.410 8 9 10 11 12 13 14 18868.449 -2545.905 -8248.322 -8180.867 -10023.281 -7817.958 -3771.034 15 16 17 18 19 20 21 -5765.505 -13305.013 7903.740 39818.319 18870.357 -16110.981 -20758.449 22 23 24 25 32131.255 15092.684 -7891.752 85148.498 > postscript(file="/var/wessaorg/rcomp/tmp/64pci1321714524.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 = 25 Frequency = 1 lag(myerror, k = 1) myerror 0 -8187.202 NA 1 -42481.117 -8187.202 2 -16034.775 -42481.117 3 -14982.147 -16034.775 4 -14946.595 -14982.147 5 -12453.990 -14946.595 6 -4328.410 -12453.990 7 18868.449 -4328.410 8 -2545.905 18868.449 9 -8248.322 -2545.905 10 -8180.867 -8248.322 11 -10023.281 -8180.867 12 -7817.958 -10023.281 13 -3771.034 -7817.958 14 -5765.505 -3771.034 15 -13305.013 -5765.505 16 7903.740 -13305.013 17 39818.319 7903.740 18 18870.357 39818.319 19 -16110.981 18870.357 20 -20758.449 -16110.981 21 32131.255 -20758.449 22 15092.684 32131.255 23 -7891.752 15092.684 24 85148.498 -7891.752 25 NA 85148.498 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -42481.117 -8187.202 [2,] -16034.775 -42481.117 [3,] -14982.147 -16034.775 [4,] -14946.595 -14982.147 [5,] -12453.990 -14946.595 [6,] -4328.410 -12453.990 [7,] 18868.449 -4328.410 [8,] -2545.905 18868.449 [9,] -8248.322 -2545.905 [10,] -8180.867 -8248.322 [11,] -10023.281 -8180.867 [12,] -7817.958 -10023.281 [13,] -3771.034 -7817.958 [14,] -5765.505 -3771.034 [15,] -13305.013 -5765.505 [16,] 7903.740 -13305.013 [17,] 39818.319 7903.740 [18,] 18870.357 39818.319 [19,] -16110.981 18870.357 [20,] -20758.449 -16110.981 [21,] 32131.255 -20758.449 [22,] 15092.684 32131.255 [23,] -7891.752 15092.684 [24,] 85148.498 -7891.752 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -42481.117 -8187.202 2 -16034.775 -42481.117 3 -14982.147 -16034.775 4 -14946.595 -14982.147 5 -12453.990 -14946.595 6 -4328.410 -12453.990 7 18868.449 -4328.410 8 -2545.905 18868.449 9 -8248.322 -2545.905 10 -8180.867 -8248.322 11 -10023.281 -8180.867 12 -7817.958 -10023.281 13 -3771.034 -7817.958 14 -5765.505 -3771.034 15 -13305.013 -5765.505 16 7903.740 -13305.013 17 39818.319 7903.740 18 18870.357 39818.319 19 -16110.981 18870.357 20 -20758.449 -16110.981 21 32131.255 -20758.449 22 15092.684 32131.255 23 -7891.752 15092.684 24 85148.498 -7891.752 > 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/wessaorg/rcomp/tmp/7n5mk1321714524.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/wessaorg/rcomp/tmp/8nwwg1321714524.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/wessaorg/rcomp/tmp/9czxf1321714524.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/wessaorg/rcomp/tmp/10cgfo1321714524.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11csq31321714524.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/wessaorg/rcomp/tmp/12wp5m1321714524.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/wessaorg/rcomp/tmp/13r3lw1321714524.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/wessaorg/rcomp/tmp/14syeg1321714524.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/wessaorg/rcomp/tmp/15401r1321714525.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/wessaorg/rcomp/tmp/16wx791321714525.tab") + } > > try(system("convert tmp/1tg3g1321714524.ps tmp/1tg3g1321714524.png",intern=TRUE)) character(0) > try(system("convert tmp/22cjz1321714524.ps tmp/22cjz1321714524.png",intern=TRUE)) character(0) > try(system("convert tmp/3cvd11321714524.ps tmp/3cvd11321714524.png",intern=TRUE)) character(0) > try(system("convert tmp/4r3z21321714524.ps tmp/4r3z21321714524.png",intern=TRUE)) character(0) > try(system("convert tmp/5q1xy1321714524.ps tmp/5q1xy1321714524.png",intern=TRUE)) character(0) > try(system("convert tmp/64pci1321714524.ps tmp/64pci1321714524.png",intern=TRUE)) character(0) > try(system("convert tmp/7n5mk1321714524.ps tmp/7n5mk1321714524.png",intern=TRUE)) character(0) > try(system("convert tmp/8nwwg1321714524.ps tmp/8nwwg1321714524.png",intern=TRUE)) character(0) > try(system("convert tmp/9czxf1321714524.ps tmp/9czxf1321714524.png",intern=TRUE)) character(0) > try(system("convert tmp/10cgfo1321714524.ps tmp/10cgfo1321714524.png",intern=TRUE)) convert: unable to open image `tmp/10cgfo1321714524.ps': No such file or directory @ blob.c/OpenBlob/2480. convert: missing an image filename `tmp/10cgfo1321714524.png' @ convert.c/ConvertImageCommand/2838. character(0) Warning message: running command 'convert tmp/10cgfo1321714524.ps tmp/10cgfo1321714524.png' had status 1 > > > proc.time() user system elapsed 2.488 0.456 2.960