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