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(1536
+ ,78
+ ,20
+ ,17
+ ,66
+ ,30
+ ,1134
+ ,46
+ ,38
+ ,17
+ ,68
+ ,42
+ ,192
+ ,18
+ ,0
+ ,0
+ ,0
+ ,0
+ ,2032
+ ,84
+ ,49
+ ,22
+ ,68
+ ,54
+ ,3230
+ ,124
+ ,74
+ ,30
+ ,120
+ ,86
+ ,5723
+ ,214
+ ,104
+ ,29
+ ,112
+ ,157
+ ,1321
+ ,49
+ ,37
+ ,19
+ ,72
+ ,36
+ ,1077
+ ,46
+ ,49
+ ,25
+ ,96
+ ,48
+ ,1462
+ ,37
+ ,42
+ ,30
+ ,109
+ ,45
+ ,2568
+ ,86
+ ,62
+ ,26
+ ,104
+ ,77
+ ,1810
+ ,69
+ ,50
+ ,20
+ ,54
+ ,49
+ ,1788
+ ,58
+ ,65
+ ,25
+ ,98
+ ,77
+ ,1334
+ ,85
+ ,28
+ ,15
+ ,49
+ ,28
+ ,2415
+ ,84
+ ,48
+ ,22
+ ,88
+ ,84
+ ,1155
+ ,43
+ ,42
+ ,12
+ ,45
+ ,31
+ ,1374
+ ,67
+ ,47
+ ,19
+ ,74
+ ,28
+ ,1503
+ ,49
+ ,71
+ ,28
+ ,112
+ ,99
+ ,999
+ ,47
+ ,0
+ ,12
+ ,45
+ ,2
+ ,2189
+ ,76
+ ,50
+ ,28
+ ,110
+ ,41
+ ,633
+ ,20
+ ,12
+ ,13
+ ,39
+ ,25
+ ,837
+ ,48
+ ,16
+ ,14
+ ,55
+ ,16
+ ,2167
+ ,81
+ ,76
+ ,27
+ ,102
+ ,96
+ ,1451
+ ,57
+ ,29
+ ,25
+ ,96
+ ,23
+ ,1790
+ ,45
+ ,38
+ ,30
+ ,86
+ ,33
+ ,1645
+ ,72
+ ,50
+ ,18
+ ,67
+ ,46
+ ,1179
+ ,22
+ ,33
+ ,17
+ ,64
+ ,59
+ ,1688
+ ,138
+ ,45
+ ,22
+ ,82
+ ,72
+ ,1100
+ ,74
+ ,59
+ ,28
+ ,100
+ ,72
+ ,2258
+ ,101
+ ,49
+ ,25
+ ,95
+ ,62
+ ,1767
+ ,35
+ ,40
+ ,16
+ ,63
+ ,55
+ ,1300
+ ,39
+ ,40
+ ,23
+ ,87
+ ,27
+ ,1432
+ ,38
+ ,51
+ ,20
+ ,65
+ ,41
+ ,1780
+ ,87
+ ,41
+ ,11
+ ,43
+ ,51
+ ,2475
+ ,102
+ ,73
+ ,20
+ ,80
+ ,26
+ ,1930
+ ,42
+ ,43
+ ,21
+ ,84
+ ,65
+ ,1
+ ,1
+ ,0
+ ,0
+ ,0
+ ,0
+ ,1782
+ ,54
+ ,46
+ ,27
+ ,105
+ ,28
+ ,1505
+ ,46
+ ,44
+ ,14
+ ,51
+ ,44
+ ,1820
+ ,41
+ ,31
+ ,29
+ ,98
+ ,36
+ ,1648
+ ,49
+ ,71
+ ,31
+ ,124
+ ,100
+ ,1668
+ ,56
+ ,61
+ ,19
+ ,75
+ ,104
+ ,1366
+ ,47
+ ,28
+ ,30
+ ,120
+ ,35
+ ,864
+ ,25
+ ,21
+ ,23
+ ,84
+ ,69
+ ,1602
+ ,62
+ ,42
+ ,20
+ ,78
+ ,73
+ ,1023
+ ,41
+ ,44
+ ,22
+ ,87
+ ,106
+ ,962
+ ,72
+ ,34
+ ,19
+ ,70
+ ,53
+ ,629
+ ,26
+ ,15
+ ,32
+ ,97
+ ,43
+ ,1568
+ ,77
+ ,46
+ ,18
+ ,72
+ ,49
+ ,1715
+ ,75
+ ,43
+ ,26
+ ,104
+ ,38
+ ,2093
+ ,51
+ ,47
+ ,25
+ ,93
+ ,51
+ ,658
+ ,28
+ ,12
+ ,22
+ ,82
+ ,14
+ ,1198
+ ,53
+ ,42
+ ,19
+ ,73
+ ,40
+ ,2059
+ ,64
+ ,56
+ ,24
+ ,87
+ ,79
+ ,1574
+ ,65
+ ,41
+ ,26
+ ,95
+ ,52
+ ,1447
+ ,48
+ ,48
+ ,27
+ ,105
+ ,44
+ ,1342
+ ,44
+ ,30
+ ,10
+ ,37
+ ,34
+ ,1526
+ ,54
+ ,44
+ ,26
+ ,96
+ ,47
+ ,669
+ ,16
+ ,25
+ ,21
+ ,80
+ ,32
+ ,859
+ ,55
+ ,42
+ ,21
+ ,83
+ ,31
+ ,2315
+ ,71
+ ,28
+ ,34
+ ,124
+ ,40
+ ,1326
+ ,47
+ ,33
+ ,29
+ ,116
+ ,42
+ ,1567
+ ,62
+ ,32
+ ,18
+ ,72
+ ,34
+ ,1080
+ ,44
+ ,28
+ ,16
+ ,55
+ ,40
+ ,896
+ ,28
+ ,31
+ ,23
+ ,86
+ ,35
+ ,855
+ ,25
+ ,13
+ ,22
+ ,85
+ ,11
+ ,1229
+ ,37
+ ,38
+ ,29
+ ,107
+ ,43
+ ,1939
+ ,60
+ ,39
+ ,31
+ ,124
+ ,53
+ ,2293
+ ,57
+ ,68
+ ,21
+ ,78
+ ,82
+ ,818
+ ,30
+ ,32
+ ,21
+ ,83
+ ,41)
+ ,dim=c(6
+ ,69)
+ ,dimnames=list(c('Pageviews'
+ ,'#Logins'
+ ,'Blogged_Computations'
+ ,'Reviewed_Compendiums'
+ ,'Feedback_in_PR'
+ ,'Included_Hyperlinks')
+ ,1:69))
> y <- array(NA,dim=c(6,69),dimnames=list(c('Pageviews','#Logins','Blogged_Computations','Reviewed_Compendiums','Feedback_in_PR','Included_Hyperlinks'),1:69))
> 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
Pageviews #Logins Blogged_Computations Reviewed_Compendiums Feedback_in_PR
1 1536 78 20 17 66
2 1134 46 38 17 68
3 192 18 0 0 0
4 2032 84 49 22 68
5 3230 124 74 30 120
6 5723 214 104 29 112
7 1321 49 37 19 72
8 1077 46 49 25 96
9 1462 37 42 30 109
10 2568 86 62 26 104
11 1810 69 50 20 54
12 1788 58 65 25 98
13 1334 85 28 15 49
14 2415 84 48 22 88
15 1155 43 42 12 45
16 1374 67 47 19 74
17 1503 49 71 28 112
18 999 47 0 12 45
19 2189 76 50 28 110
20 633 20 12 13 39
21 837 48 16 14 55
22 2167 81 76 27 102
23 1451 57 29 25 96
24 1790 45 38 30 86
25 1645 72 50 18 67
26 1179 22 33 17 64
27 1688 138 45 22 82
28 1100 74 59 28 100
29 2258 101 49 25 95
30 1767 35 40 16 63
31 1300 39 40 23 87
32 1432 38 51 20 65
33 1780 87 41 11 43
34 2475 102 73 20 80
35 1930 42 43 21 84
36 1 1 0 0 0
37 1782 54 46 27 105
38 1505 46 44 14 51
39 1820 41 31 29 98
40 1648 49 71 31 124
41 1668 56 61 19 75
42 1366 47 28 30 120
43 864 25 21 23 84
44 1602 62 42 20 78
45 1023 41 44 22 87
46 962 72 34 19 70
47 629 26 15 32 97
48 1568 77 46 18 72
49 1715 75 43 26 104
50 2093 51 47 25 93
51 658 28 12 22 82
52 1198 53 42 19 73
53 2059 64 56 24 87
54 1574 65 41 26 95
55 1447 48 48 27 105
56 1342 44 30 10 37
57 1526 54 44 26 96
58 669 16 25 21 80
59 859 55 42 21 83
60 2315 71 28 34 124
61 1326 47 33 29 116
62 1567 62 32 18 72
63 1080 44 28 16 55
64 896 28 31 23 86
65 855 25 13 22 85
66 1229 37 38 29 107
67 1939 60 39 31 124
68 2293 57 68 21 78
69 818 30 32 21 83
Included_Hyperlinks
1 30
2 42
3 0
4 54
5 86
6 157
7 36
8 48
9 45
10 77
11 49
12 77
13 28
14 84
15 31
16 28
17 99
18 2
19 41
20 25
21 16
22 96
23 23
24 33
25 46
26 59
27 72
28 72
29 62
30 55
31 27
32 41
33 51
34 26
35 65
36 0
37 28
38 44
39 36
40 100
41 104
42 35
43 69
44 73
45 106
46 53
47 43
48 49
49 38
50 51
51 14
52 40
53 79
54 52
55 44
56 34
57 47
58 32
59 31
60 40
61 42
62 34
63 40
64 35
65 11
66 43
67 53
68 82
69 41
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `#Logins` Blogged_Computations
-146.191 14.594 11.220
Reviewed_Compendiums Feedback_in_PR Included_Hyperlinks
21.763 -2.146 1.758
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1114.02 -125.23 -24.04 189.76 912.48
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -146.191 152.901 -0.956 0.3427
`#Logins` 14.594 1.846 7.906 5.21e-11 ***
Blogged_Computations 11.220 4.319 2.598 0.0117 *
Reviewed_Compendiums 21.763 24.824 0.877 0.3840
Feedback_in_PR -2.146 6.609 -0.325 0.7464
Included_Hyperlinks 1.758 2.582 0.681 0.4986
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 356.7 on 63 degrees of freedom
Multiple R-squared: 0.7991, Adjusted R-squared: 0.7832
F-statistic: 50.13 on 5 and 63 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.29427967 0.58855934 0.70572033
[2,] 0.18333264 0.36666528 0.81666736
[3,] 0.09676495 0.19352989 0.90323505
[4,] 0.06672194 0.13344389 0.93327806
[5,] 0.07686355 0.15372711 0.92313645
[6,] 0.07517626 0.15035251 0.92482374
[7,] 0.05838383 0.11676766 0.94161617
[8,] 0.03380764 0.06761529 0.96619236
[9,] 0.12786544 0.25573089 0.87213456
[10,] 0.11339151 0.22678303 0.88660849
[11,] 0.11397799 0.22795598 0.88602201
[12,] 0.08194520 0.16389041 0.91805480
[13,] 0.05922479 0.11844958 0.94077521
[14,] 0.05050429 0.10100857 0.94949571
[15,] 0.03177582 0.06355164 0.96822418
[16,] 0.03776472 0.07552944 0.96223528
[17,] 0.02354591 0.04709183 0.97645409
[18,] 0.02307591 0.04615181 0.97692409
[19,] 0.71663572 0.56672856 0.28336428
[20,] 0.96032071 0.07935857 0.03967929
[21,] 0.94126690 0.11746620 0.05873310
[22,] 0.97277518 0.05444963 0.02722482
[23,] 0.96056136 0.07887728 0.03943864
[24,] 0.94328154 0.11343693 0.05671846
[25,] 0.91921307 0.16157387 0.08078693
[26,] 0.89345808 0.21308383 0.10654192
[27,] 0.94357076 0.11285848 0.05642924
[28,] 0.92475793 0.15048413 0.07524207
[29,] 0.90321066 0.19357868 0.09678934
[30,] 0.87963741 0.24072517 0.12036259
[31,] 0.90999013 0.18001975 0.09000987
[32,] 0.90995483 0.18009034 0.09004517
[33,] 0.87793012 0.24413975 0.12206988
[34,] 0.83280910 0.33438180 0.16719090
[35,] 0.78642546 0.42714907 0.21357454
[36,] 0.72792160 0.54415680 0.27207840
[37,] 0.74486961 0.51026078 0.25513039
[38,] 0.91827494 0.16345013 0.08172506
[39,] 0.95769047 0.08461907 0.04230953
[40,] 0.94860219 0.10279563 0.05139781
[41,] 0.91992491 0.16015018 0.08007509
[42,] 0.97131458 0.05737084 0.02868542
[43,] 0.95138875 0.09722250 0.04861125
[44,] 0.92509455 0.14981089 0.07490545
[45,] 0.89702017 0.20595966 0.10297983
[46,] 0.92553263 0.14893474 0.07446737
[47,] 0.89472167 0.21055665 0.10527833
[48,] 0.88352349 0.23295302 0.11647651
[49,] 0.80741637 0.38516725 0.19258363
[50,] 0.69706258 0.60587484 0.30293742
[51,] 0.84852984 0.30294032 0.15147016
[52,] 0.77452142 0.45095716 0.22547858
> postscript(file="/var/wessaorg/rcomp/tmp/1o4d61323867229.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/29dl21323867229.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/32eio1323867229.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/4a2dg1323867229.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/55l4f1323867229.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 = 69
Frequency = 1
1 2 3 4 5 6
38.419519 -115.338913 75.500320 -25.238633 189.758242 912.483703
7 8 9 10 11 12
14.706392 -420.311331 98.938255 285.499819 -17.284422 -110.645866
13 14 15 16 17 18
-344.947078 359.174674 -16.654925 -288.831315 -405.532471 291.194481
19 20 21 22 23 24
219.715951 109.517704 -111.593816 -259.067345 61.504762 326.794506
25 26 27 28 29 30
-149.364311 297.548325 -1114.022238 -1017.034182 -68.734160 643.933276
31 32 33 34 35 36
66.950345 83.580443 -40.254699 4.286561 589.797152 132.597388
37 38 39 40 41 42
212.546029 213.623987 535.957591 -301.821950 -122.828573 55.270022
43 44 45 46 47 48
-31.818139 -24.039647 -401.223275 -680.469429 -336.349389 -248.994772
49 50 51 52 53 54
-125.228034 533.446739 -66.468149 -227.655286 168.409172 -141.775092
55 56 57 58 59 60
-85.455802 311.474527 -51.966821 -40.374171 -602.084481 566.760764
61 62 63 64 65 66
-39.958494 152.364501 -30.573604 -91.743382 174.806092 -68.195316
67 68 69
270.308567 410.621730 -183.612029
> postscript(file="/var/wessaorg/rcomp/tmp/6yr8u1323867229.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 38.419519 NA
1 -115.338913 38.419519
2 75.500320 -115.338913
3 -25.238633 75.500320
4 189.758242 -25.238633
5 912.483703 189.758242
6 14.706392 912.483703
7 -420.311331 14.706392
8 98.938255 -420.311331
9 285.499819 98.938255
10 -17.284422 285.499819
11 -110.645866 -17.284422
12 -344.947078 -110.645866
13 359.174674 -344.947078
14 -16.654925 359.174674
15 -288.831315 -16.654925
16 -405.532471 -288.831315
17 291.194481 -405.532471
18 219.715951 291.194481
19 109.517704 219.715951
20 -111.593816 109.517704
21 -259.067345 -111.593816
22 61.504762 -259.067345
23 326.794506 61.504762
24 -149.364311 326.794506
25 297.548325 -149.364311
26 -1114.022238 297.548325
27 -1017.034182 -1114.022238
28 -68.734160 -1017.034182
29 643.933276 -68.734160
30 66.950345 643.933276
31 83.580443 66.950345
32 -40.254699 83.580443
33 4.286561 -40.254699
34 589.797152 4.286561
35 132.597388 589.797152
36 212.546029 132.597388
37 213.623987 212.546029
38 535.957591 213.623987
39 -301.821950 535.957591
40 -122.828573 -301.821950
41 55.270022 -122.828573
42 -31.818139 55.270022
43 -24.039647 -31.818139
44 -401.223275 -24.039647
45 -680.469429 -401.223275
46 -336.349389 -680.469429
47 -248.994772 -336.349389
48 -125.228034 -248.994772
49 533.446739 -125.228034
50 -66.468149 533.446739
51 -227.655286 -66.468149
52 168.409172 -227.655286
53 -141.775092 168.409172
54 -85.455802 -141.775092
55 311.474527 -85.455802
56 -51.966821 311.474527
57 -40.374171 -51.966821
58 -602.084481 -40.374171
59 566.760764 -602.084481
60 -39.958494 566.760764
61 152.364501 -39.958494
62 -30.573604 152.364501
63 -91.743382 -30.573604
64 174.806092 -91.743382
65 -68.195316 174.806092
66 270.308567 -68.195316
67 410.621730 270.308567
68 -183.612029 410.621730
69 NA -183.612029
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -115.338913 38.419519
[2,] 75.500320 -115.338913
[3,] -25.238633 75.500320
[4,] 189.758242 -25.238633
[5,] 912.483703 189.758242
[6,] 14.706392 912.483703
[7,] -420.311331 14.706392
[8,] 98.938255 -420.311331
[9,] 285.499819 98.938255
[10,] -17.284422 285.499819
[11,] -110.645866 -17.284422
[12,] -344.947078 -110.645866
[13,] 359.174674 -344.947078
[14,] -16.654925 359.174674
[15,] -288.831315 -16.654925
[16,] -405.532471 -288.831315
[17,] 291.194481 -405.532471
[18,] 219.715951 291.194481
[19,] 109.517704 219.715951
[20,] -111.593816 109.517704
[21,] -259.067345 -111.593816
[22,] 61.504762 -259.067345
[23,] 326.794506 61.504762
[24,] -149.364311 326.794506
[25,] 297.548325 -149.364311
[26,] -1114.022238 297.548325
[27,] -1017.034182 -1114.022238
[28,] -68.734160 -1017.034182
[29,] 643.933276 -68.734160
[30,] 66.950345 643.933276
[31,] 83.580443 66.950345
[32,] -40.254699 83.580443
[33,] 4.286561 -40.254699
[34,] 589.797152 4.286561
[35,] 132.597388 589.797152
[36,] 212.546029 132.597388
[37,] 213.623987 212.546029
[38,] 535.957591 213.623987
[39,] -301.821950 535.957591
[40,] -122.828573 -301.821950
[41,] 55.270022 -122.828573
[42,] -31.818139 55.270022
[43,] -24.039647 -31.818139
[44,] -401.223275 -24.039647
[45,] -680.469429 -401.223275
[46,] -336.349389 -680.469429
[47,] -248.994772 -336.349389
[48,] -125.228034 -248.994772
[49,] 533.446739 -125.228034
[50,] -66.468149 533.446739
[51,] -227.655286 -66.468149
[52,] 168.409172 -227.655286
[53,] -141.775092 168.409172
[54,] -85.455802 -141.775092
[55,] 311.474527 -85.455802
[56,] -51.966821 311.474527
[57,] -40.374171 -51.966821
[58,] -602.084481 -40.374171
[59,] 566.760764 -602.084481
[60,] -39.958494 566.760764
[61,] 152.364501 -39.958494
[62,] -30.573604 152.364501
[63,] -91.743382 -30.573604
[64,] 174.806092 -91.743382
[65,] -68.195316 174.806092
[66,] 270.308567 -68.195316
[67,] 410.621730 270.308567
[68,] -183.612029 410.621730
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -115.338913 38.419519
2 75.500320 -115.338913
3 -25.238633 75.500320
4 189.758242 -25.238633
5 912.483703 189.758242
6 14.706392 912.483703
7 -420.311331 14.706392
8 98.938255 -420.311331
9 285.499819 98.938255
10 -17.284422 285.499819
11 -110.645866 -17.284422
12 -344.947078 -110.645866
13 359.174674 -344.947078
14 -16.654925 359.174674
15 -288.831315 -16.654925
16 -405.532471 -288.831315
17 291.194481 -405.532471
18 219.715951 291.194481
19 109.517704 219.715951
20 -111.593816 109.517704
21 -259.067345 -111.593816
22 61.504762 -259.067345
23 326.794506 61.504762
24 -149.364311 326.794506
25 297.548325 -149.364311
26 -1114.022238 297.548325
27 -1017.034182 -1114.022238
28 -68.734160 -1017.034182
29 643.933276 -68.734160
30 66.950345 643.933276
31 83.580443 66.950345
32 -40.254699 83.580443
33 4.286561 -40.254699
34 589.797152 4.286561
35 132.597388 589.797152
36 212.546029 132.597388
37 213.623987 212.546029
38 535.957591 213.623987
39 -301.821950 535.957591
40 -122.828573 -301.821950
41 55.270022 -122.828573
42 -31.818139 55.270022
43 -24.039647 -31.818139
44 -401.223275 -24.039647
45 -680.469429 -401.223275
46 -336.349389 -680.469429
47 -248.994772 -336.349389
48 -125.228034 -248.994772
49 533.446739 -125.228034
50 -66.468149 533.446739
51 -227.655286 -66.468149
52 168.409172 -227.655286
53 -141.775092 168.409172
54 -85.455802 -141.775092
55 311.474527 -85.455802
56 -51.966821 311.474527
57 -40.374171 -51.966821
58 -602.084481 -40.374171
59 566.760764 -602.084481
60 -39.958494 566.760764
61 152.364501 -39.958494
62 -30.573604 152.364501
63 -91.743382 -30.573604
64 174.806092 -91.743382
65 -68.195316 174.806092
66 270.308567 -68.195316
67 410.621730 270.308567
68 -183.612029 410.621730
> 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/77g871323867229.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/8ild61323867229.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/9qrdz1323867229.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/104y271323867229.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/11dkvn1323867229.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/12n84s1323867229.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/13t97r1323867229.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/140gm81323867229.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/15xuff1323867229.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/16ftm61323867229.tab")
+ }
>
> try(system("convert tmp/1o4d61323867229.ps tmp/1o4d61323867229.png",intern=TRUE))
character(0)
> try(system("convert tmp/29dl21323867229.ps tmp/29dl21323867229.png",intern=TRUE))
character(0)
> try(system("convert tmp/32eio1323867229.ps tmp/32eio1323867229.png",intern=TRUE))
character(0)
> try(system("convert tmp/4a2dg1323867229.ps tmp/4a2dg1323867229.png",intern=TRUE))
character(0)
> try(system("convert tmp/55l4f1323867229.ps tmp/55l4f1323867229.png",intern=TRUE))
character(0)
> try(system("convert tmp/6yr8u1323867229.ps tmp/6yr8u1323867229.png",intern=TRUE))
character(0)
> try(system("convert tmp/77g871323867229.ps tmp/77g871323867229.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ild61323867229.ps tmp/8ild61323867229.png",intern=TRUE))
character(0)
> try(system("convert tmp/9qrdz1323867229.ps tmp/9qrdz1323867229.png",intern=TRUE))
character(0)
> try(system("convert tmp/104y271323867229.ps tmp/104y271323867229.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.304 0.586 3.910