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 = '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
> 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
test tijd login vieuws revieuws size shared blogged intrisieke
1 0 96560 76 129 17 22996 0 78 62
2 3 112611 41 36 20 26706 7 44 56
3 0 98146 40 37 17 27114 1 80 49
4 2 121848 39 45 17 30594 0 73 63
5 1 31774 23 48 17 4143 2 107 67
6 1 65475 18 24 13 69008 0 42 59
7 0 108446 60 90 17 46300 2 76 40
8 2 76302 31 26 20 30976 2 69 34
9 -2 30989 14 35 17 4154 0 62 37
10 -4 150580 77 124 22 45588 0 46 61
11 1 59382 49 49 12 26263 0 133 60
12 0 341570 168 190 21 117105 0 71 57
13 -3 133328 55 79 20 40909 0 46 56
14 0 101523 42 76 22 61056 1 131 67
15 -2 92499 32 57 18 21399 2 47 38
16 -2 118612 46 72 12 30080 0 15 49
17 -3 98104 54 132 17 25568 0 37 32
18 0 84105 20 45 17 20055 1 0 63
19 0 237213 84 74 38 66198 3 79 67
20 4 133131 55 52 30 57793 3 77 43
21 1 344297 75 86 30 67654 5 101 84
22 3 174415 100 63 31 82753 0 105 49
23 4 294424 77 59 33 101494 2 124 58
24 2 362301 119 715 34 100708 0 83 63
25 0 167488 45 77 28 83737 0 106 29
26 2 152299 53 63 33 61370 2 25 58
27 0 243511 71 65 42 101338 2 16 62
28 3 132487 41 97 36 40735 1 22 54
29 3 172494 52 52 43 86687 0 29 53
30 0 224330 83 52 39 130115 0 5 66
31 6 181633 70 48 30 64466 0 27 53
32 2 210907 56 81 30 112285 0 29 26
33 0 236785 119 75 31 101481 5 43 43
34 2 244052 68 66 44 143558 0 158 53
35 4 143756 46 57 34 69094 4 102 54
36 2 182079 63 63 33 102860 0 123 47
37 3 100750 72 67 30 140867 0 105 43
38 0 152474 65 45 32 65567 1 33 57
39 -1 97839 38 42 24 94785 2 96 41
40 0 149061 44 66 26 116174 6 56 37
41 0 324799 154 108 47 97668 5 59 52
42 -1 230964 53 43 30 133824 0 91 52
43 0 174724 92 135 34 69112 1 76 67
44 -1 223632 73 52 33 72654 1 94 70
45 1 106408 30 32 14 31081 2 41 68
46 -2 265769 146 37 32 83122 5 67 43
47 -4 149112 56 65 35 60578 2 100 56
48 2 152871 58 74 28 79892 5 67 74
49 -4 183167 66 66 39 82875 1 135 58
50 2 218946 41 112 29 80670 4 58 63
51 -3 196553 57 50 29 95260 0 56 64
52 2 143246 103 42 27 106671 0 59 53
53 2 193339 78 47 35 84651 1 116 51
54 -4 130585 46 57 29 95364 2 98 54
55 3 148446 91 63 37 126846 8 32 48
56 -1 243060 63 110 29 111813 4 63 50
57 -3 317394 86 53 31 91413 0 113 45
58 0 244749 95 144 33 76643 1 111 61
59 2 128423 64 89 38 92696 10 120 56
60 2 229242 247 128 31 91721 0 25 46
61 2 324598 110 128 37 135777 1 109 51
62 -2 195838 67 50 31 102372 0 37 37
63 0 254488 83 50 39 103772 2 54 42
64 -3 271856 103 91 37 54990 2 55 69
65 3 95227 34 70 32 34777 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 revieuws size
-6.075e-01 -9.729e-06 2.354e-03 3.110e-03 4.835e-02 1.237e-05
shared blogged intrisieke extrisieke
1.243e-01 -1.823e-03 1.976e-02 -1.686e-02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4.9461 -1.4097 -0.0083 1.4743 6.0118
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.075e-01 2.991e+00 -0.203 0.840
tijd -9.729e-06 6.603e-06 -1.473 0.146
login 2.354e-03 1.042e-02 0.226 0.822
vieuws 3.110e-03 3.898e-03 0.798 0.428
revieuws 4.835e-02 5.153e-02 0.938 0.352
size 1.237e-05 1.263e-05 0.980 0.332
shared 1.243e-01 1.368e-01 0.909 0.367
blogged -1.823e-03 8.279e-03 -0.220 0.827
intrisieke 1.976e-02 3.071e-02 0.643 0.523
extrisieke -1.686e-02 4.425e-02 -0.381 0.705
Residual standard error: 2.333 on 55 degrees of freedom
Multiple R-squared: 0.09188, Adjusted R-squared: -0.05672
F-statistic: 0.6183 on 9 and 55 DF, p-value: 0.7763
> 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.40279121 0.80558242 0.5972088
[2,] 0.24451689 0.48903377 0.7554831
[3,] 0.18522125 0.37044251 0.8147787
[4,] 0.11681381 0.23362762 0.8831862
[5,] 0.07587613 0.15175226 0.9241239
[6,] 0.07472314 0.14944628 0.9252769
[7,] 0.04056160 0.08112321 0.9594384
[8,] 0.08738981 0.17477962 0.9126102
[9,] 0.06281498 0.12562996 0.9371850
[10,] 0.04089031 0.08178062 0.9591097
[11,] 0.06281137 0.12562275 0.9371886
[12,] 0.20112450 0.40224901 0.7988755
[13,] 0.16696094 0.33392188 0.8330391
[14,] 0.12864155 0.25728309 0.8713585
[15,] 0.10529946 0.21059891 0.8947005
[16,] 0.12813965 0.25627929 0.8718604
[17,] 0.10027830 0.20055660 0.8997217
[18,] 0.08422503 0.16845006 0.9157750
[19,] 0.47481782 0.94963564 0.5251822
[20,] 0.41764694 0.83529387 0.5823531
[21,] 0.36157825 0.72315650 0.6384218
[22,] 0.36011634 0.72023269 0.6398837
[23,] 0.45420646 0.90841292 0.5457935
[24,] 0.46213206 0.92426411 0.5378679
[25,] 0.42784671 0.85569343 0.5721533
[26,] 0.35229213 0.70458425 0.6477079
[27,] 0.31757543 0.63515085 0.6824246
[28,] 0.27127302 0.54254605 0.7287270
[29,] 0.22015344 0.44030687 0.7798466
[30,] 0.18110886 0.36221772 0.8188911
[31,] 0.15739738 0.31479475 0.8426026
[32,] 0.12600722 0.25201444 0.8739928
[33,] 0.09870729 0.19741458 0.9012927
[34,] 0.07225829 0.14451658 0.9277417
[35,] 0.13747279 0.27494559 0.8625272
[36,] 0.14958543 0.29917086 0.8504146
[37,] 0.21269333 0.42538666 0.7873067
[38,] 0.31384801 0.62769602 0.6861520
[39,] 0.24933937 0.49867874 0.7506606
[40,] 0.44670655 0.89341311 0.5532934
> postscript(file="/var/wessaorg/rcomp/tmp/1n6d61323620161.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/29emc1323620161.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/364h81323620161.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/4pq2y1323620161.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/5yqt91323620161.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
-0.00828895 2.38022125 0.36209036 2.56439618 0.66000471 0.75324262
7 8 9 10 11 12
-0.05883400 2.14727167 -1.91475470 -4.02937616 1.18228975 0.62990450
13 14 15 16 17 18
-2.70151362 -0.53740226 -1.79902247 -1.24975057 -2.58273072 -0.03679742
19 20 21 22 23 24
-0.52416996 3.27450901 1.54518661 2.68656010 4.31842790 0.80851677
25 26 27 28 29 30
0.30233708 1.27887473 -1.05916587 2.21555916 1.86747476 -1.24452655
31 32 33 34 35 36
6.01178976 1.81816946 -0.87239214 1.13418840 2.85542525 1.74366708
37 38 39 40 41 42
1.36996126 -0.57036179 -1.61281933 -1.40974404 -0.85693859 -1.14075699
43 44 45 46 47 48
-1.42926787 -1.28464834 1.15854655 -2.07712008 -4.70981964 0.79037740
49 50 51 52 53 54
-4.70899174 1.40306169 -3.38325957 1.34477044 1.69590652 -4.94611602
55 56 57 58 59 60
0.49494521 -1.52729418 -1.80441389 -0.27456671 -0.50488862 1.47429476
61 62 63 64 65
1.70127551 -1.91418303 -0.28719402 -3.12410902 2.21197238
> postscript(file="/var/wessaorg/rcomp/tmp/6wwdj1323620161.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 -0.00828895 NA
1 2.38022125 -0.00828895
2 0.36209036 2.38022125
3 2.56439618 0.36209036
4 0.66000471 2.56439618
5 0.75324262 0.66000471
6 -0.05883400 0.75324262
7 2.14727167 -0.05883400
8 -1.91475470 2.14727167
9 -4.02937616 -1.91475470
10 1.18228975 -4.02937616
11 0.62990450 1.18228975
12 -2.70151362 0.62990450
13 -0.53740226 -2.70151362
14 -1.79902247 -0.53740226
15 -1.24975057 -1.79902247
16 -2.58273072 -1.24975057
17 -0.03679742 -2.58273072
18 -0.52416996 -0.03679742
19 3.27450901 -0.52416996
20 1.54518661 3.27450901
21 2.68656010 1.54518661
22 4.31842790 2.68656010
23 0.80851677 4.31842790
24 0.30233708 0.80851677
25 1.27887473 0.30233708
26 -1.05916587 1.27887473
27 2.21555916 -1.05916587
28 1.86747476 2.21555916
29 -1.24452655 1.86747476
30 6.01178976 -1.24452655
31 1.81816946 6.01178976
32 -0.87239214 1.81816946
33 1.13418840 -0.87239214
34 2.85542525 1.13418840
35 1.74366708 2.85542525
36 1.36996126 1.74366708
37 -0.57036179 1.36996126
38 -1.61281933 -0.57036179
39 -1.40974404 -1.61281933
40 -0.85693859 -1.40974404
41 -1.14075699 -0.85693859
42 -1.42926787 -1.14075699
43 -1.28464834 -1.42926787
44 1.15854655 -1.28464834
45 -2.07712008 1.15854655
46 -4.70981964 -2.07712008
47 0.79037740 -4.70981964
48 -4.70899174 0.79037740
49 1.40306169 -4.70899174
50 -3.38325957 1.40306169
51 1.34477044 -3.38325957
52 1.69590652 1.34477044
53 -4.94611602 1.69590652
54 0.49494521 -4.94611602
55 -1.52729418 0.49494521
56 -1.80441389 -1.52729418
57 -0.27456671 -1.80441389
58 -0.50488862 -0.27456671
59 1.47429476 -0.50488862
60 1.70127551 1.47429476
61 -1.91418303 1.70127551
62 -0.28719402 -1.91418303
63 -3.12410902 -0.28719402
64 2.21197238 -3.12410902
65 NA 2.21197238
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2.38022125 -0.00828895
[2,] 0.36209036 2.38022125
[3,] 2.56439618 0.36209036
[4,] 0.66000471 2.56439618
[5,] 0.75324262 0.66000471
[6,] -0.05883400 0.75324262
[7,] 2.14727167 -0.05883400
[8,] -1.91475470 2.14727167
[9,] -4.02937616 -1.91475470
[10,] 1.18228975 -4.02937616
[11,] 0.62990450 1.18228975
[12,] -2.70151362 0.62990450
[13,] -0.53740226 -2.70151362
[14,] -1.79902247 -0.53740226
[15,] -1.24975057 -1.79902247
[16,] -2.58273072 -1.24975057
[17,] -0.03679742 -2.58273072
[18,] -0.52416996 -0.03679742
[19,] 3.27450901 -0.52416996
[20,] 1.54518661 3.27450901
[21,] 2.68656010 1.54518661
[22,] 4.31842790 2.68656010
[23,] 0.80851677 4.31842790
[24,] 0.30233708 0.80851677
[25,] 1.27887473 0.30233708
[26,] -1.05916587 1.27887473
[27,] 2.21555916 -1.05916587
[28,] 1.86747476 2.21555916
[29,] -1.24452655 1.86747476
[30,] 6.01178976 -1.24452655
[31,] 1.81816946 6.01178976
[32,] -0.87239214 1.81816946
[33,] 1.13418840 -0.87239214
[34,] 2.85542525 1.13418840
[35,] 1.74366708 2.85542525
[36,] 1.36996126 1.74366708
[37,] -0.57036179 1.36996126
[38,] -1.61281933 -0.57036179
[39,] -1.40974404 -1.61281933
[40,] -0.85693859 -1.40974404
[41,] -1.14075699 -0.85693859
[42,] -1.42926787 -1.14075699
[43,] -1.28464834 -1.42926787
[44,] 1.15854655 -1.28464834
[45,] -2.07712008 1.15854655
[46,] -4.70981964 -2.07712008
[47,] 0.79037740 -4.70981964
[48,] -4.70899174 0.79037740
[49,] 1.40306169 -4.70899174
[50,] -3.38325957 1.40306169
[51,] 1.34477044 -3.38325957
[52,] 1.69590652 1.34477044
[53,] -4.94611602 1.69590652
[54,] 0.49494521 -4.94611602
[55,] -1.52729418 0.49494521
[56,] -1.80441389 -1.52729418
[57,] -0.27456671 -1.80441389
[58,] -0.50488862 -0.27456671
[59,] 1.47429476 -0.50488862
[60,] 1.70127551 1.47429476
[61,] -1.91418303 1.70127551
[62,] -0.28719402 -1.91418303
[63,] -3.12410902 -0.28719402
[64,] 2.21197238 -3.12410902
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2.38022125 -0.00828895
2 0.36209036 2.38022125
3 2.56439618 0.36209036
4 0.66000471 2.56439618
5 0.75324262 0.66000471
6 -0.05883400 0.75324262
7 2.14727167 -0.05883400
8 -1.91475470 2.14727167
9 -4.02937616 -1.91475470
10 1.18228975 -4.02937616
11 0.62990450 1.18228975
12 -2.70151362 0.62990450
13 -0.53740226 -2.70151362
14 -1.79902247 -0.53740226
15 -1.24975057 -1.79902247
16 -2.58273072 -1.24975057
17 -0.03679742 -2.58273072
18 -0.52416996 -0.03679742
19 3.27450901 -0.52416996
20 1.54518661 3.27450901
21 2.68656010 1.54518661
22 4.31842790 2.68656010
23 0.80851677 4.31842790
24 0.30233708 0.80851677
25 1.27887473 0.30233708
26 -1.05916587 1.27887473
27 2.21555916 -1.05916587
28 1.86747476 2.21555916
29 -1.24452655 1.86747476
30 6.01178976 -1.24452655
31 1.81816946 6.01178976
32 -0.87239214 1.81816946
33 1.13418840 -0.87239214
34 2.85542525 1.13418840
35 1.74366708 2.85542525
36 1.36996126 1.74366708
37 -0.57036179 1.36996126
38 -1.61281933 -0.57036179
39 -1.40974404 -1.61281933
40 -0.85693859 -1.40974404
41 -1.14075699 -0.85693859
42 -1.42926787 -1.14075699
43 -1.28464834 -1.42926787
44 1.15854655 -1.28464834
45 -2.07712008 1.15854655
46 -4.70981964 -2.07712008
47 0.79037740 -4.70981964
48 -4.70899174 0.79037740
49 1.40306169 -4.70899174
50 -3.38325957 1.40306169
51 1.34477044 -3.38325957
52 1.69590652 1.34477044
53 -4.94611602 1.69590652
54 0.49494521 -4.94611602
55 -1.52729418 0.49494521
56 -1.80441389 -1.52729418
57 -0.27456671 -1.80441389
58 -0.50488862 -0.27456671
59 1.47429476 -0.50488862
60 1.70127551 1.47429476
61 -1.91418303 1.70127551
62 -0.28719402 -1.91418303
63 -3.12410902 -0.28719402
64 2.21197238 -3.12410902
> 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/7k05e1323620161.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/8n3qb1323620161.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/9a4hk1323620161.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/10i4wp1323620161.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/118eyx1323620161.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/1240w31323620161.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/13rpdn1323620161.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/14mr761323620161.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/15rdle1323620161.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/16hbqr1323620161.tab")
+ }
>
> try(system("convert tmp/1n6d61323620161.ps tmp/1n6d61323620161.png",intern=TRUE))
character(0)
> try(system("convert tmp/29emc1323620161.ps tmp/29emc1323620161.png",intern=TRUE))
character(0)
> try(system("convert tmp/364h81323620161.ps tmp/364h81323620161.png",intern=TRUE))
character(0)
> try(system("convert tmp/4pq2y1323620161.ps tmp/4pq2y1323620161.png",intern=TRUE))
character(0)
> try(system("convert tmp/5yqt91323620161.ps tmp/5yqt91323620161.png",intern=TRUE))
character(0)
> try(system("convert tmp/6wwdj1323620161.ps tmp/6wwdj1323620161.png",intern=TRUE))
character(0)
> try(system("convert tmp/7k05e1323620161.ps tmp/7k05e1323620161.png",intern=TRUE))
character(0)
> try(system("convert tmp/8n3qb1323620161.ps tmp/8n3qb1323620161.png",intern=TRUE))
character(0)
> try(system("convert tmp/9a4hk1323620161.ps tmp/9a4hk1323620161.png",intern=TRUE))
character(0)
> try(system("convert tmp/10i4wp1323620161.ps tmp/10i4wp1323620161.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.455 0.506 4.024