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(96560 + ,76 + ,129 + ,17 + ,22996 + ,0 + ,0 + ,78 + ,62 + ,72 + ,112611 + ,41 + ,36 + ,20 + ,26706 + ,3 + ,7 + ,44 + ,56 + ,64 + ,98146 + ,40 + ,37 + ,17 + ,27114 + ,0 + ,1 + ,80 + ,49 + ,66 + ,121848 + ,39 + ,45 + ,17 + ,30594 + ,2 + ,0 + ,73 + ,63 + ,78 + ,31774 + ,23 + ,48 + ,17 + ,4143 + ,1 + ,2 + ,107 + ,67 + ,71 + ,65475 + ,18 + ,24 + ,13 + ,69008 + ,1 + ,0 + ,42 + ,59 + ,71 + ,108446 + ,60 + ,90 + ,17 + ,46300 + ,0 + ,2 + ,76 + ,40 + ,59 + ,76302 + ,31 + ,26 + ,20 + ,30976 + ,2 + ,2 + ,69 + ,34 + ,65 + ,30989 + ,14 + ,35 + ,17 + ,4154 + ,-2 + ,0 + ,62 + ,37 + ,48 + ,150580 + ,77 + ,124 + ,22 + ,45588 + ,-4 + ,0 + ,46 + ,61 + ,72 + ,59382 + ,49 + ,49 + ,12 + ,26263 + ,1 + ,0 + ,133 + ,60 + ,66 + ,341570 + ,168 + ,190 + ,21 + ,117105 + ,0 + ,0 + ,71 + ,57 + ,68 + ,133328 + ,55 + ,79 + ,20 + ,40909 + ,-3 + ,0 + ,46 + ,56 + ,75 + ,101523 + ,42 + ,76 + ,22 + ,61056 + ,0 + ,1 + ,131 + ,67 + ,73 + ,92499 + ,32 + ,57 + ,18 + ,21399 + ,-2 + ,2 + ,47 + ,38 + ,59 + ,118612 + ,46 + ,72 + ,12 + ,30080 + ,-2 + ,0 + ,15 + ,49 + ,72 + ,98104 + ,54 + ,132 + ,17 + ,25568 + ,-3 + ,0 + ,37 + ,32 + ,65 + ,84105 + ,20 + ,45 + ,17 + ,20055 + ,0 + ,1 + ,0 + ,63 + ,69 + ,237213 + ,84 + ,74 + ,38 + ,66198 + ,0 + ,3 + ,79 + ,67 + ,71 + ,133131 + ,55 + ,52 + ,30 + ,57793 + ,4 + ,3 + ,77 + ,43 + ,54 + ,344297 + ,75 + ,86 + ,30 + ,67654 + ,1 + ,5 + ,101 + ,84 + ,84 + ,174415 + ,100 + ,63 + ,31 + ,82753 + ,3 + ,0 + ,105 + ,49 + ,66 + ,294424 + ,77 + ,59 + ,33 + ,101494 + ,4 + ,2 + ,124 + ,58 + ,73 + ,362301 + ,119 + ,715 + ,34 + ,100708 + ,2 + ,0 + ,83 + ,63 + ,69 + ,167488 + ,45 + ,77 + ,28 + ,83737 + ,0 + ,0 + ,106 + ,29 + ,70 + ,152299 + ,53 + ,63 + ,33 + ,61370 + ,2 + ,2 + ,25 + ,58 + ,72 + ,243511 + ,71 + ,65 + ,42 + ,101338 + ,0 + ,2 + ,16 + ,62 + ,63 + ,132487 + ,41 + ,97 + ,36 + ,40735 + ,3 + ,1 + ,22 + ,54 + ,66 + ,172494 + ,52 + ,52 + ,43 + ,86687 + ,3 + ,0 + ,29 + ,53 + ,60 + ,224330 + ,83 + ,52 + ,39 + ,130115 + ,0 + ,0 + ,5 + ,66 + ,66 + ,181633 + ,70 + ,48 + ,30 + ,64466 + ,6 + ,0 + ,27 + ,53 + ,71 + ,210907 + ,56 + ,81 + ,30 + ,112285 + ,2 + ,0 + ,29 + ,26 + ,50 + ,236785 + ,119 + ,75 + ,31 + ,101481 + ,0 + ,5 + ,43 + ,43 + ,52 + ,244052 + ,68 + ,66 + ,44 + ,143558 + ,2 + ,0 + ,158 + ,53 + ,70 + ,143756 + ,46 + ,57 + ,34 + ,69094 + ,4 + ,4 + ,102 + ,54 + ,60 + ,182079 + ,63 + ,63 + ,33 + ,102860 + ,2 + ,0 + ,123 + ,47 + ,76 + ,100750 + ,72 + ,67 + ,30 + ,140867 + ,3 + ,0 + ,105 + ,43 + ,60 + ,152474 + ,65 + ,45 + ,32 + ,65567 + ,0 + ,1 + ,33 + ,57 + ,70 + ,97839 + ,38 + ,42 + ,24 + ,94785 + ,-1 + ,2 + ,96 + ,41 + ,75 + ,149061 + ,44 + ,66 + ,26 + ,116174 + ,0 + ,6 + ,56 + ,37 + ,54 + ,324799 + ,154 + ,108 + ,47 + ,97668 + ,0 + ,5 + ,59 + ,52 + ,65 + ,230964 + ,53 + ,43 + ,30 + ,133824 + ,-1 + ,0 + ,91 + ,52 + ,73 + ,174724 + ,92 + ,135 + ,34 + ,69112 + ,0 + ,1 + ,76 + ,67 + ,42 + ,223632 + ,73 + ,52 + ,33 + ,72654 + ,-1 + ,1 + ,94 + ,70 + ,65 + ,106408 + ,30 + ,32 + ,14 + ,31081 + ,1 + ,2 + ,41 + ,68 + ,75 + ,265769 + ,146 + ,37 + ,32 + ,83122 + ,-2 + ,5 + ,67 + ,43 + ,66 + ,149112 + ,56 + ,65 + ,35 + ,60578 + ,-4 + ,2 + ,100 + ,56 + ,70 + ,152871 + ,58 + ,74 + ,28 + ,79892 + ,2 + ,5 + ,67 + ,74 + ,81 + ,183167 + ,66 + ,66 + ,39 + ,82875 + ,-4 + ,1 + ,135 + ,58 + ,71 + ,218946 + ,41 + ,112 + ,29 + ,80670 + ,2 + ,4 + ,58 + ,63 + ,68 + ,196553 + ,57 + ,50 + ,29 + ,95260 + ,-3 + ,0 + ,56 + ,64 + ,67 + ,143246 + ,103 + ,42 + ,27 + ,106671 + ,2 + ,0 + ,59 + ,53 + ,76 + ,193339 + ,78 + ,47 + ,35 + ,84651 + ,2 + ,1 + ,116 + ,51 + ,71 + ,130585 + ,46 + ,57 + ,29 + ,95364 + ,-4 + ,2 + ,98 + ,54 + ,70 + ,148446 + ,91 + ,63 + ,37 + ,126846 + ,3 + ,8 + ,32 + ,48 + ,65 + ,243060 + ,63 + ,110 + ,29 + ,111813 + ,-1 + ,4 + ,63 + ,50 + ,68 + ,317394 + ,86 + ,53 + ,31 + ,91413 + ,-3 + ,0 + ,113 + ,45 + ,70 + ,244749 + ,95 + ,144 + ,33 + ,76643 + ,0 + ,1 + ,111 + ,61 + ,64 + ,128423 + ,64 + ,89 + ,38 + ,92696 + ,2 + ,10 + ,120 + ,56 + ,70 + ,229242 + ,247 + ,128 + ,31 + ,91721 + ,2 + ,0 + ,25 + ,46 + ,66 + ,324598 + ,110 + ,128 + ,37 + ,135777 + ,2 + ,1 + ,109 + ,51 + ,59 + ,195838 + ,67 + ,50 + ,31 + ,102372 + ,-2 + ,0 + ,37 + ,37 + ,78 + ,254488 + ,83 + ,50 + ,39 + ,103772 + ,0 + ,2 + ,54 + ,42 + ,67 + ,271856 + ,103 + ,91 + ,37 + ,54990 + ,-3 + ,2 + ,55 + ,69 + ,67 + ,95227 + ,34 + ,70 + ,32 + ,34777 + ,3 + ,0 + ,17 + ,56 + ,61) + ,dim=c(10 + ,65) + ,dimnames=list(c('tijd' + ,'login' + ,'vieuws' + ,'revieuws' + ,'size' + ,'test' + ,'shared' + ,'blogged' + ,'intrisieke' + ,'extrisieke') + ,1:65)) > y <- array(NA,dim=c(10,65),dimnames=list(c('tijd','login','vieuws','revieuws','size','test','shared','blogged','intrisieke','extrisieke'),1:65)) > 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 = '4' > 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 revieuws tijd login vieuws size test shared blogged intrisieke 1 17 96560 76 129 22996 0 0 78 62 2 20 112611 41 36 26706 3 7 44 56 3 17 98146 40 37 27114 0 1 80 49 4 17 121848 39 45 30594 2 0 73 63 5 17 31774 23 48 4143 1 2 107 67 6 13 65475 18 24 69008 1 0 42 59 7 17 108446 60 90 46300 0 2 76 40 8 20 76302 31 26 30976 2 2 69 34 9 17 30989 14 35 4154 -2 0 62 37 10 22 150580 77 124 45588 -4 0 46 61 11 12 59382 49 49 26263 1 0 133 60 12 21 341570 168 190 117105 0 0 71 57 13 20 133328 55 79 40909 -3 0 46 56 14 22 101523 42 76 61056 0 1 131 67 15 18 92499 32 57 21399 -2 2 47 38 16 12 118612 46 72 30080 -2 0 15 49 17 17 98104 54 132 25568 -3 0 37 32 18 17 84105 20 45 20055 0 1 0 63 19 38 237213 84 74 66198 0 3 79 67 20 30 133131 55 52 57793 4 3 77 43 21 30 344297 75 86 67654 1 5 101 84 22 31 174415 100 63 82753 3 0 105 49 23 33 294424 77 59 101494 4 2 124 58 24 34 362301 119 715 100708 2 0 83 63 25 28 167488 45 77 83737 0 0 106 29 26 33 152299 53 63 61370 2 2 25 58 27 42 243511 71 65 101338 0 2 16 62 28 36 132487 41 97 40735 3 1 22 54 29 43 172494 52 52 86687 3 0 29 53 30 39 224330 83 52 130115 0 0 5 66 31 30 181633 70 48 64466 6 0 27 53 32 30 210907 56 81 112285 2 0 29 26 33 31 236785 119 75 101481 0 5 43 43 34 44 244052 68 66 143558 2 0 158 53 35 34 143756 46 57 69094 4 4 102 54 36 33 182079 63 63 102860 2 0 123 47 37 30 100750 72 67 140867 3 0 105 43 38 32 152474 65 45 65567 0 1 33 57 39 24 97839 38 42 94785 -1 2 96 41 40 26 149061 44 66 116174 0 6 56 37 41 47 324799 154 108 97668 0 5 59 52 42 30 230964 53 43 133824 -1 0 91 52 43 34 174724 92 135 69112 0 1 76 67 44 33 223632 73 52 72654 -1 1 94 70 45 14 106408 30 32 31081 1 2 41 68 46 32 265769 146 37 83122 -2 5 67 43 47 35 149112 56 65 60578 -4 2 100 56 48 28 152871 58 74 79892 2 5 67 74 49 39 183167 66 66 82875 -4 1 135 58 50 29 218946 41 112 80670 2 4 58 63 51 29 196553 57 50 95260 -3 0 56 64 52 27 143246 103 42 106671 2 0 59 53 53 35 193339 78 47 84651 2 1 116 51 54 29 130585 46 57 95364 -4 2 98 54 55 37 148446 91 63 126846 3 8 32 48 56 29 243060 63 110 111813 -1 4 63 50 57 31 317394 86 53 91413 -3 0 113 45 58 33 244749 95 144 76643 0 1 111 61 59 38 128423 64 89 92696 2 10 120 56 60 31 229242 247 128 91721 2 0 25 46 61 37 324598 110 128 135777 2 1 109 51 62 31 195838 67 50 102372 -2 0 37 37 63 39 254488 83 50 103772 0 2 54 42 64 37 271856 103 91 54990 -3 2 55 69 65 32 95227 34 70 34777 3 0 17 56 extrisieke 1 72 2 64 3 66 4 78 5 71 6 71 7 59 8 65 9 48 10 72 11 66 12 68 13 75 14 73 15 59 16 72 17 65 18 69 19 71 20 54 21 84 22 66 23 73 24 69 25 70 26 72 27 63 28 66 29 60 30 66 31 71 32 50 33 52 34 70 35 60 36 76 37 60 38 70 39 75 40 54 41 65 42 73 43 42 44 65 45 75 46 66 47 70 48 81 49 71 50 68 51 67 52 76 53 71 54 70 55 65 56 68 57 70 58 64 59 70 60 66 61 59 62 78 63 67 64 67 65 61 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) tijd login vieuws size test 2.182e+01 4.462e-05 -5.236e-03 -1.225e-02 9.976e-05 3.258e-01 shared blogged intrisieke extrisieke 3.658e-01 6.224e-04 1.092e-01 -2.043e-01 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -16.909 -4.186 0.266 3.565 11.840 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.182e+01 7.187e+00 3.036 0.00366 ** tijd 4.462e-05 1.641e-05 2.719 0.00873 ** login -5.236e-03 2.704e-02 -0.194 0.84719 vieuws -1.225e-02 1.004e-02 -1.220 0.22762 size 9.976e-05 3.022e-05 3.301 0.00169 ** test 3.258e-01 3.473e-01 0.938 0.35225 shared 3.658e-01 3.543e-01 1.032 0.30642 blogged 6.224e-04 2.150e-02 0.029 0.97701 intrisieke 1.092e-01 7.865e-02 1.388 0.17065 extrisieke -2.043e-01 1.117e-01 -1.829 0.07280 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 6.057 on 55 degrees of freedom Multiple R-squared: 0.5676, Adjusted R-squared: 0.4968 F-statistic: 8.021 on 9 and 55 DF, p-value: 1.816e-07 > 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.04807880 0.09615759 0.95192120 [2,] 0.03029046 0.06058093 0.96970954 [3,] 0.03149581 0.06299162 0.96850419 [4,] 0.08279664 0.16559328 0.91720336 [5,] 0.05348136 0.10696271 0.94651864 [6,] 0.03610286 0.07220572 0.96389714 [7,] 0.20528535 0.41057071 0.79471465 [8,] 0.33939918 0.67879836 0.66060082 [9,] 0.52334929 0.95330142 0.47665071 [10,] 0.56613405 0.86773189 0.43386595 [11,] 0.50151551 0.99696897 0.49848449 [12,] 0.46299812 0.92599623 0.53700188 [13,] 0.38651977 0.77303955 0.61348023 [14,] 0.51947649 0.96104702 0.48052351 [15,] 0.59141871 0.81716257 0.40858129 [16,] 0.73400614 0.53198771 0.26599386 [17,] 0.84345641 0.31308718 0.15654359 [18,] 0.84051920 0.31896160 0.15948080 [19,] 0.78688599 0.42622802 0.21311401 [20,] 0.79163293 0.41673414 0.20836707 [21,] 0.74956693 0.50086614 0.25043307 [22,] 0.78630970 0.42738059 0.21369030 [23,] 0.72537077 0.54925846 0.27462923 [24,] 0.66983125 0.66033751 0.33016875 [25,] 0.59855223 0.80289553 0.40144777 [26,] 0.58694826 0.82610347 0.41305174 [27,] 0.52380681 0.95238637 0.47619319 [28,] 0.62821459 0.74357081 0.37178541 [29,] 0.82513299 0.34973402 0.17486701 [30,] 0.77030401 0.45939199 0.22969599 [31,] 0.72794740 0.54410520 0.27205260 [32,] 0.64173687 0.71652627 0.35826313 [33,] 0.90540179 0.18919641 0.09459821 [34,] 0.94547508 0.10904984 0.05452492 [35,] 0.92938198 0.14123603 0.07061802 [36,] 0.87536758 0.24926483 0.12463242 [37,] 0.92879239 0.14241523 0.07120761 [38,] 0.90749607 0.18500786 0.09250393 [39,] 0.82126308 0.35747384 0.17873692 [40,] 0.76581884 0.46836233 0.23418116 > postscript(file="/var/wessaorg/rcomp/tmp/11onx1323622899.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/2z2tp1323622899.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/3n2451323622899.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/4yc4z1323622899.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/5umh81323622899.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 = 65 Frequency = 1 1 2 3 4 5 6 -1.5513905 -5.4560134 -3.5219096 -4.1924769 0.1245580 -10.5245475 7 8 9 10 11 12 -5.9519275 0.3086220 0.2659772 0.1608185 -7.7072982 -16.9092836 13 14 15 16 17 18 -0.4361750 -2.1372656 -1.4192943 -6.9872385 0.8931564 -3.0638747 19 20 21 22 23 24 6.3830464 1.2899020 -6.7117651 1.5302554 -4.4852195 1.8672371 25 26 27 28 29 30 2.6023233 6.2911587 5.7349262 11.8402170 11.2223103 1.5381482 31 32 33 34 35 36 1.3474465 -4.4380069 -5.8921010 5.8998918 3.3459802 3.5646007 37 38 39 40 41 42 -2.6482087 5.4199358 -0.5458408 -8.2567082 8.8117631 -5.1643400 43 44 45 46 47 48 0.4789592 0.5125730 -8.3031223 -1.1815383 10.2704087 -1.4528963 49 50 51 52 53 54 10.9211479 -4.1854753 -2.5378514 -1.7460020 4.9415631 1.6962353 55 56 57 58 59 60 1.2676223 -5.8837745 -2.7054738 0.8567922 6.4323036 0.4584340 61 62 63 64 65 -5.3023450 3.7181214 4.8585722 5.5856514 8.8587054 > postscript(file="/var/wessaorg/rcomp/tmp/6694n1323622899.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 = 65 Frequency = 1 lag(myerror, k = 1) myerror 0 -1.5513905 NA 1 -5.4560134 -1.5513905 2 -3.5219096 -5.4560134 3 -4.1924769 -3.5219096 4 0.1245580 -4.1924769 5 -10.5245475 0.1245580 6 -5.9519275 -10.5245475 7 0.3086220 -5.9519275 8 0.2659772 0.3086220 9 0.1608185 0.2659772 10 -7.7072982 0.1608185 11 -16.9092836 -7.7072982 12 -0.4361750 -16.9092836 13 -2.1372656 -0.4361750 14 -1.4192943 -2.1372656 15 -6.9872385 -1.4192943 16 0.8931564 -6.9872385 17 -3.0638747 0.8931564 18 6.3830464 -3.0638747 19 1.2899020 6.3830464 20 -6.7117651 1.2899020 21 1.5302554 -6.7117651 22 -4.4852195 1.5302554 23 1.8672371 -4.4852195 24 2.6023233 1.8672371 25 6.2911587 2.6023233 26 5.7349262 6.2911587 27 11.8402170 5.7349262 28 11.2223103 11.8402170 29 1.5381482 11.2223103 30 1.3474465 1.5381482 31 -4.4380069 1.3474465 32 -5.8921010 -4.4380069 33 5.8998918 -5.8921010 34 3.3459802 5.8998918 35 3.5646007 3.3459802 36 -2.6482087 3.5646007 37 5.4199358 -2.6482087 38 -0.5458408 5.4199358 39 -8.2567082 -0.5458408 40 8.8117631 -8.2567082 41 -5.1643400 8.8117631 42 0.4789592 -5.1643400 43 0.5125730 0.4789592 44 -8.3031223 0.5125730 45 -1.1815383 -8.3031223 46 10.2704087 -1.1815383 47 -1.4528963 10.2704087 48 10.9211479 -1.4528963 49 -4.1854753 10.9211479 50 -2.5378514 -4.1854753 51 -1.7460020 -2.5378514 52 4.9415631 -1.7460020 53 1.6962353 4.9415631 54 1.2676223 1.6962353 55 -5.8837745 1.2676223 56 -2.7054738 -5.8837745 57 0.8567922 -2.7054738 58 6.4323036 0.8567922 59 0.4584340 6.4323036 60 -5.3023450 0.4584340 61 3.7181214 -5.3023450 62 4.8585722 3.7181214 63 5.5856514 4.8585722 64 8.8587054 5.5856514 65 NA 8.8587054 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -5.4560134 -1.5513905 [2,] -3.5219096 -5.4560134 [3,] -4.1924769 -3.5219096 [4,] 0.1245580 -4.1924769 [5,] -10.5245475 0.1245580 [6,] -5.9519275 -10.5245475 [7,] 0.3086220 -5.9519275 [8,] 0.2659772 0.3086220 [9,] 0.1608185 0.2659772 [10,] -7.7072982 0.1608185 [11,] -16.9092836 -7.7072982 [12,] -0.4361750 -16.9092836 [13,] -2.1372656 -0.4361750 [14,] -1.4192943 -2.1372656 [15,] -6.9872385 -1.4192943 [16,] 0.8931564 -6.9872385 [17,] -3.0638747 0.8931564 [18,] 6.3830464 -3.0638747 [19,] 1.2899020 6.3830464 [20,] -6.7117651 1.2899020 [21,] 1.5302554 -6.7117651 [22,] -4.4852195 1.5302554 [23,] 1.8672371 -4.4852195 [24,] 2.6023233 1.8672371 [25,] 6.2911587 2.6023233 [26,] 5.7349262 6.2911587 [27,] 11.8402170 5.7349262 [28,] 11.2223103 11.8402170 [29,] 1.5381482 11.2223103 [30,] 1.3474465 1.5381482 [31,] -4.4380069 1.3474465 [32,] -5.8921010 -4.4380069 [33,] 5.8998918 -5.8921010 [34,] 3.3459802 5.8998918 [35,] 3.5646007 3.3459802 [36,] -2.6482087 3.5646007 [37,] 5.4199358 -2.6482087 [38,] -0.5458408 5.4199358 [39,] -8.2567082 -0.5458408 [40,] 8.8117631 -8.2567082 [41,] -5.1643400 8.8117631 [42,] 0.4789592 -5.1643400 [43,] 0.5125730 0.4789592 [44,] -8.3031223 0.5125730 [45,] -1.1815383 -8.3031223 [46,] 10.2704087 -1.1815383 [47,] -1.4528963 10.2704087 [48,] 10.9211479 -1.4528963 [49,] -4.1854753 10.9211479 [50,] -2.5378514 -4.1854753 [51,] -1.7460020 -2.5378514 [52,] 4.9415631 -1.7460020 [53,] 1.6962353 4.9415631 [54,] 1.2676223 1.6962353 [55,] -5.8837745 1.2676223 [56,] -2.7054738 -5.8837745 [57,] 0.8567922 -2.7054738 [58,] 6.4323036 0.8567922 [59,] 0.4584340 6.4323036 [60,] -5.3023450 0.4584340 [61,] 3.7181214 -5.3023450 [62,] 4.8585722 3.7181214 [63,] 5.5856514 4.8585722 [64,] 8.8587054 5.5856514 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -5.4560134 -1.5513905 2 -3.5219096 -5.4560134 3 -4.1924769 -3.5219096 4 0.1245580 -4.1924769 5 -10.5245475 0.1245580 6 -5.9519275 -10.5245475 7 0.3086220 -5.9519275 8 0.2659772 0.3086220 9 0.1608185 0.2659772 10 -7.7072982 0.1608185 11 -16.9092836 -7.7072982 12 -0.4361750 -16.9092836 13 -2.1372656 -0.4361750 14 -1.4192943 -2.1372656 15 -6.9872385 -1.4192943 16 0.8931564 -6.9872385 17 -3.0638747 0.8931564 18 6.3830464 -3.0638747 19 1.2899020 6.3830464 20 -6.7117651 1.2899020 21 1.5302554 -6.7117651 22 -4.4852195 1.5302554 23 1.8672371 -4.4852195 24 2.6023233 1.8672371 25 6.2911587 2.6023233 26 5.7349262 6.2911587 27 11.8402170 5.7349262 28 11.2223103 11.8402170 29 1.5381482 11.2223103 30 1.3474465 1.5381482 31 -4.4380069 1.3474465 32 -5.8921010 -4.4380069 33 5.8998918 -5.8921010 34 3.3459802 5.8998918 35 3.5646007 3.3459802 36 -2.6482087 3.5646007 37 5.4199358 -2.6482087 38 -0.5458408 5.4199358 39 -8.2567082 -0.5458408 40 8.8117631 -8.2567082 41 -5.1643400 8.8117631 42 0.4789592 -5.1643400 43 0.5125730 0.4789592 44 -8.3031223 0.5125730 45 -1.1815383 -8.3031223 46 10.2704087 -1.1815383 47 -1.4528963 10.2704087 48 10.9211479 -1.4528963 49 -4.1854753 10.9211479 50 -2.5378514 -4.1854753 51 -1.7460020 -2.5378514 52 4.9415631 -1.7460020 53 1.6962353 4.9415631 54 1.2676223 1.6962353 55 -5.8837745 1.2676223 56 -2.7054738 -5.8837745 57 0.8567922 -2.7054738 58 6.4323036 0.8567922 59 0.4584340 6.4323036 60 -5.3023450 0.4584340 61 3.7181214 -5.3023450 62 4.8585722 3.7181214 63 5.5856514 4.8585722 64 8.8587054 5.5856514 > 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/7sm4x1323622899.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/83l841323622899.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/9av191323622899.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/10yvzl1323622899.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/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/11x4rw1323622899.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/12hwwc1323622899.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/13idkp1323622899.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/14mpdf1323622899.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/15bhcf1323622899.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/162ljc1323622899.tab") + } > > try(system("convert tmp/11onx1323622899.ps tmp/11onx1323622899.png",intern=TRUE)) character(0) > try(system("convert tmp/2z2tp1323622899.ps tmp/2z2tp1323622899.png",intern=TRUE)) character(0) > try(system("convert tmp/3n2451323622899.ps tmp/3n2451323622899.png",intern=TRUE)) character(0) > try(system("convert tmp/4yc4z1323622899.ps tmp/4yc4z1323622899.png",intern=TRUE)) character(0) > try(system("convert tmp/5umh81323622899.ps tmp/5umh81323622899.png",intern=TRUE)) character(0) > try(system("convert tmp/6694n1323622899.ps tmp/6694n1323622899.png",intern=TRUE)) character(0) > try(system("convert tmp/7sm4x1323622899.ps tmp/7sm4x1323622899.png",intern=TRUE)) character(0) > try(system("convert tmp/83l841323622899.ps tmp/83l841323622899.png",intern=TRUE)) character(0) > try(system("convert tmp/9av191323622899.ps tmp/9av191323622899.png",intern=TRUE)) character(0) > try(system("convert tmp/10yvzl1323622899.ps tmp/10yvzl1323622899.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.509 0.534 4.091