R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(210907
+ ,56
+ ,396
+ ,81
+ ,3
+ ,79
+ ,30
+ ,120982
+ ,56
+ ,297
+ ,55
+ ,4
+ ,58
+ ,28
+ ,176508
+ ,54
+ ,559
+ ,50
+ ,12
+ ,60
+ ,38
+ ,179321
+ ,89
+ ,967
+ ,125
+ ,2
+ ,108
+ ,30
+ ,123185
+ ,40
+ ,270
+ ,40
+ ,1
+ ,49
+ ,22
+ ,52746
+ ,25
+ ,143
+ ,37
+ ,3
+ ,0
+ ,26
+ ,385534
+ ,92
+ ,1562
+ ,63
+ ,0
+ ,121
+ ,25
+ ,33170
+ ,18
+ ,109
+ ,44
+ ,0
+ ,1
+ ,18
+ ,101645
+ ,63
+ ,371
+ ,88
+ ,0
+ ,20
+ ,11
+ ,149061
+ ,44
+ ,656
+ ,66
+ ,5
+ ,43
+ ,26
+ ,165446
+ ,33
+ ,511
+ ,57
+ ,0
+ ,69
+ ,25
+ ,237213
+ ,84
+ ,655
+ ,74
+ ,0
+ ,78
+ ,38
+ ,173326
+ ,88
+ ,465
+ ,49
+ ,7
+ ,86
+ ,44
+ ,133131
+ ,55
+ ,525
+ ,52
+ ,7
+ ,44
+ ,30
+ ,258873
+ ,60
+ ,885
+ ,88
+ ,3
+ ,104
+ ,40
+ ,180083
+ ,66
+ ,497
+ ,36
+ ,9
+ ,63
+ ,34
+ ,324799
+ ,154
+ ,1436
+ ,108
+ ,0
+ ,158
+ ,47
+ ,230964
+ ,53
+ ,612
+ ,43
+ ,4
+ ,102
+ ,30
+ ,236785
+ ,119
+ ,865
+ ,75
+ ,3
+ ,77
+ ,31
+ ,135473
+ ,41
+ ,385
+ ,32
+ ,0
+ ,82
+ ,23
+ ,202925
+ ,61
+ ,567
+ ,44
+ ,7
+ ,115
+ ,36
+ ,215147
+ ,58
+ ,639
+ ,85
+ ,0
+ ,101
+ ,36
+ ,344297
+ ,75
+ ,963
+ ,86
+ ,1
+ ,80
+ ,30
+ ,153935
+ ,33
+ ,398
+ ,56
+ ,5
+ ,50
+ ,25
+ ,132943
+ ,40
+ ,410
+ ,50
+ ,7
+ ,83
+ ,39
+ ,174724
+ ,92
+ ,966
+ ,135
+ ,0
+ ,123
+ ,34
+ ,174415
+ ,100
+ ,801
+ ,63
+ ,0
+ ,73
+ ,31
+ ,225548
+ ,112
+ ,892
+ ,81
+ ,5
+ ,81
+ ,31
+ ,223632
+ ,73
+ ,513
+ ,52
+ ,0
+ ,105
+ ,33
+ ,124817
+ ,40
+ ,469
+ ,44
+ ,0
+ ,47
+ ,25
+ ,221698
+ ,45
+ ,683
+ ,113
+ ,0
+ ,105
+ ,33
+ ,210767
+ ,60
+ ,643
+ ,39
+ ,3
+ ,94
+ ,35
+ ,170266
+ ,62
+ ,535
+ ,73
+ ,4
+ ,44
+ ,42
+ ,260561
+ ,75
+ ,625
+ ,48
+ ,1
+ ,114
+ ,43
+ ,84853
+ ,31
+ ,264
+ ,33
+ ,4
+ ,38
+ ,30
+ ,294424
+ ,77
+ ,992
+ ,59
+ ,2
+ ,107
+ ,33
+ ,101011
+ ,34
+ ,238
+ ,41
+ ,0
+ ,30
+ ,13
+ ,215641
+ ,46
+ ,818
+ ,69
+ ,0
+ ,71
+ ,32
+ ,325107
+ ,99
+ ,937
+ ,64
+ ,0
+ ,84
+ ,36
+ ,7176
+ ,17
+ ,70
+ ,1
+ ,0
+ ,0
+ ,0
+ ,167542
+ ,66
+ ,507
+ ,59
+ ,2
+ ,59
+ ,28
+ ,106408
+ ,30
+ ,260
+ ,32
+ ,1
+ ,33
+ ,14
+ ,96560
+ ,76
+ ,503
+ ,129
+ ,0
+ ,42
+ ,17
+ ,265769
+ ,146
+ ,927
+ ,37
+ ,2
+ ,96
+ ,32
+ ,269651
+ ,67
+ ,1269
+ ,31
+ ,10
+ ,106
+ ,30
+ ,149112
+ ,56
+ ,537
+ ,65
+ ,6
+ ,56
+ ,35
+ ,175824
+ ,107
+ ,910
+ ,107
+ ,0
+ ,57
+ ,20
+ ,152871
+ ,58
+ ,532
+ ,74
+ ,5
+ ,59
+ ,28
+ ,111665
+ ,34
+ ,345
+ ,54
+ ,4
+ ,39
+ ,28
+ ,116408
+ ,61
+ ,918
+ ,76
+ ,1
+ ,34
+ ,39
+ ,362301
+ ,119
+ ,1635
+ ,715
+ ,2
+ ,76
+ ,34
+ ,78800
+ ,42
+ ,330
+ ,57
+ ,2
+ ,20
+ ,26
+ ,183167
+ ,66
+ ,557
+ ,66
+ ,0
+ ,91
+ ,39
+ ,277965
+ ,89
+ ,1178
+ ,106
+ ,8
+ ,115
+ ,39
+ ,150629
+ ,44
+ ,740
+ ,54
+ ,3
+ ,85
+ ,33
+ ,168809
+ ,66
+ ,452
+ ,32
+ ,0
+ ,76
+ ,28
+ ,24188
+ ,24
+ ,218
+ ,20
+ ,0
+ ,8
+ ,4)
+ ,dim=c(7
+ ,57)
+ ,dimnames=list(c('time_in_rfc'
+ ,'logins'
+ ,'compendium_views_info'
+ ,'compendium_views_pr'
+ ,'shared_compendiums'
+ ,'blogged_computations'
+ ,'compendiums_reviewed')
+ ,1:57))
> y <- array(NA,dim=c(7,57),dimnames=list(c('time_in_rfc','logins','compendium_views_info','compendium_views_pr','shared_compendiums','blogged_computations','compendiums_reviewed'),1:57))
> 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'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, as.Date.numeric
> 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
time_in_rfc logins compendium_views_info compendium_views_pr
1 210907 56 396 81
2 120982 56 297 55
3 176508 54 559 50
4 179321 89 967 125
5 123185 40 270 40
6 52746 25 143 37
7 385534 92 1562 63
8 33170 18 109 44
9 101645 63 371 88
10 149061 44 656 66
11 165446 33 511 57
12 237213 84 655 74
13 173326 88 465 49
14 133131 55 525 52
15 258873 60 885 88
16 180083 66 497 36
17 324799 154 1436 108
18 230964 53 612 43
19 236785 119 865 75
20 135473 41 385 32
21 202925 61 567 44
22 215147 58 639 85
23 344297 75 963 86
24 153935 33 398 56
25 132943 40 410 50
26 174724 92 966 135
27 174415 100 801 63
28 225548 112 892 81
29 223632 73 513 52
30 124817 40 469 44
31 221698 45 683 113
32 210767 60 643 39
33 170266 62 535 73
34 260561 75 625 48
35 84853 31 264 33
36 294424 77 992 59
37 101011 34 238 41
38 215641 46 818 69
39 325107 99 937 64
40 7176 17 70 1
41 167542 66 507 59
42 106408 30 260 32
43 96560 76 503 129
44 265769 146 927 37
45 269651 67 1269 31
46 149112 56 537 65
47 175824 107 910 107
48 152871 58 532 74
49 111665 34 345 54
50 116408 61 918 76
51 362301 119 1635 715
52 78800 42 330 57
53 183167 66 557 66
54 277965 89 1178 106
55 150629 44 740 54
56 168809 66 452 32
57 24188 24 218 20
shared_compendiums blogged_computations compendiums_reviewed
1 3 79 30
2 4 58 28
3 12 60 38
4 2 108 30
5 1 49 22
6 3 0 26
7 0 121 25
8 0 1 18
9 0 20 11
10 5 43 26
11 0 69 25
12 0 78 38
13 7 86 44
14 7 44 30
15 3 104 40
16 9 63 34
17 0 158 47
18 4 102 30
19 3 77 31
20 0 82 23
21 7 115 36
22 0 101 36
23 1 80 30
24 5 50 25
25 7 83 39
26 0 123 34
27 0 73 31
28 5 81 31
29 0 105 33
30 0 47 25
31 0 105 33
32 3 94 35
33 4 44 42
34 1 114 43
35 4 38 30
36 2 107 33
37 0 30 13
38 0 71 32
39 0 84 36
40 0 0 0
41 2 59 28
42 1 33 14
43 0 42 17
44 2 96 32
45 10 106 30
46 6 56 35
47 0 57 20
48 5 59 28
49 4 39 28
50 1 34 39
51 2 76 34
52 2 20 26
53 0 91 39
54 8 115 39
55 3 85 33
56 0 76 28
57 0 8 4
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) logins compendium_views_info
12038.80 97.57 121.47
compendium_views_pr shared_compendiums blogged_computations
37.00 -545.19 877.08
compendiums_reviewed
759.61
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-102330 -15061 756 18100 112377
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12038.80 17490.87 0.688 0.49445
logins 97.57 267.12 0.365 0.71644
compendium_views_info 121.47 28.68 4.235 9.78e-05 ***
compendium_views_pr 37.00 69.61 0.532 0.59743
shared_compendiums -545.19 1882.31 -0.290 0.77329
blogged_computations 877.08 260.28 3.370 0.00146 **
compendiums_reviewed 759.61 817.36 0.929 0.35717
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 37540 on 50 degrees of freedom
Multiple R-squared: 0.8149, Adjusted R-squared: 0.7927
F-statistic: 36.69 on 6 and 50 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.7555253 0.48894936 0.244474680
[2,] 0.6146624 0.77067527 0.385337636
[3,] 0.5633366 0.87332675 0.436663373
[4,] 0.5760361 0.84792787 0.423963933
[5,] 0.4593405 0.91868094 0.540659530
[6,] 0.3622323 0.72446452 0.637767741
[7,] 0.2679735 0.53594696 0.732026522
[8,] 0.3101473 0.62029457 0.689852715
[9,] 0.2368533 0.47370669 0.763146653
[10,] 0.1928960 0.38579200 0.807104001
[11,] 0.2184325 0.43686498 0.781567511
[12,] 0.1638437 0.32768745 0.836156273
[13,] 0.1150401 0.23008015 0.884959926
[14,] 0.6708171 0.65836581 0.329182905
[15,] 0.6270878 0.74582448 0.372912240
[16,] 0.6138349 0.77233028 0.386165142
[17,] 0.9321175 0.13576493 0.067882464
[18,] 0.9460560 0.10788790 0.053943952
[19,] 0.9188230 0.16235395 0.081176977
[20,] 0.8924934 0.21501317 0.107506586
[21,] 0.8627809 0.27443819 0.137219094
[22,] 0.8289579 0.34208414 0.171042070
[23,] 0.7684099 0.46318026 0.231590130
[24,] 0.7448796 0.51024089 0.255120445
[25,] 0.6820924 0.63581519 0.317907595
[26,] 0.6265976 0.74680479 0.373402397
[27,] 0.5702701 0.85945976 0.429729880
[28,] 0.5075029 0.98499420 0.492497102
[29,] 0.4611212 0.92224249 0.538878753
[30,] 0.9825781 0.03484372 0.017421862
[31,] 0.9699836 0.06003274 0.030016369
[32,] 0.9577348 0.08453034 0.042265168
[33,] 0.9769293 0.04614145 0.023070726
[34,] 0.9942463 0.01150734 0.005753671
[35,] 0.9842302 0.03153957 0.015769783
[36,] 0.9948028 0.01039436 0.005197178
[37,] 0.9865818 0.02683646 0.013418228
[38,] 0.9723875 0.05522500 0.027612498
> postscript(file="/var/fisher/rcomp/tmp/1yzdu1355155700.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/fisher/rcomp/tmp/20tdc1355155700.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/fisher/rcomp/tmp/3b6kl1355155700.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/fisher/rcomp/tmp/4wram1355155700.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/fisher/rcomp/tmp/54nre1355155700.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 = 57
Frequency = 1
1 2 3 4 5 6
51864.5110 -4590.3153 14503.0877 -79906.7568 13824.2667 1414.9119
7 8 9 10 11 12
47341.0535 -10042.9669 9242.0284 -4132.8897 6500.4613 37402.4152
13 14 15 16 17 18
-10629.4589 -7531.4032 10261.3963 23727.6680 -54967.3791 27755.9425
19 20 21 22 23 24
15844.4991 -17906.3449 -9958.8697 756.1881 112377.1691 28142.5729
25 26 27 28 29 30
-33255.9462 -102330.0619 -34581.0590 -628.9763 23073.6909 -9933.2444
31 32 33 34 35 36
965.7411 5931.8636 16177.9964 31406.7090 -17435.2246 34370.6320
37 38 39 40 41 42
19041.4965 10621.6218 86206.2370 -15061.1591 13370.7126 19643.8737
43 44 45 46 47 48
-38515.3031 18099.7050 -14519.1521 -8454.7842 -26333.6942 -2475.7247
49 50 51 52 53 54
-889.4755 -74800.6493 22206.1124 -15730.6996 -14849.2989 -15893.9604
55 56 57
-55569.1132 6316.5078 -27467.1603
> postscript(file="/var/fisher/rcomp/tmp/6uysl1355155700.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 = 57
Frequency = 1
lag(myerror, k = 1) myerror
0 51864.5110 NA
1 -4590.3153 51864.5110
2 14503.0877 -4590.3153
3 -79906.7568 14503.0877
4 13824.2667 -79906.7568
5 1414.9119 13824.2667
6 47341.0535 1414.9119
7 -10042.9669 47341.0535
8 9242.0284 -10042.9669
9 -4132.8897 9242.0284
10 6500.4613 -4132.8897
11 37402.4152 6500.4613
12 -10629.4589 37402.4152
13 -7531.4032 -10629.4589
14 10261.3963 -7531.4032
15 23727.6680 10261.3963
16 -54967.3791 23727.6680
17 27755.9425 -54967.3791
18 15844.4991 27755.9425
19 -17906.3449 15844.4991
20 -9958.8697 -17906.3449
21 756.1881 -9958.8697
22 112377.1691 756.1881
23 28142.5729 112377.1691
24 -33255.9462 28142.5729
25 -102330.0619 -33255.9462
26 -34581.0590 -102330.0619
27 -628.9763 -34581.0590
28 23073.6909 -628.9763
29 -9933.2444 23073.6909
30 965.7411 -9933.2444
31 5931.8636 965.7411
32 16177.9964 5931.8636
33 31406.7090 16177.9964
34 -17435.2246 31406.7090
35 34370.6320 -17435.2246
36 19041.4965 34370.6320
37 10621.6218 19041.4965
38 86206.2370 10621.6218
39 -15061.1591 86206.2370
40 13370.7126 -15061.1591
41 19643.8737 13370.7126
42 -38515.3031 19643.8737
43 18099.7050 -38515.3031
44 -14519.1521 18099.7050
45 -8454.7842 -14519.1521
46 -26333.6942 -8454.7842
47 -2475.7247 -26333.6942
48 -889.4755 -2475.7247
49 -74800.6493 -889.4755
50 22206.1124 -74800.6493
51 -15730.6996 22206.1124
52 -14849.2989 -15730.6996
53 -15893.9604 -14849.2989
54 -55569.1132 -15893.9604
55 6316.5078 -55569.1132
56 -27467.1603 6316.5078
57 NA -27467.1603
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -4590.3153 51864.5110
[2,] 14503.0877 -4590.3153
[3,] -79906.7568 14503.0877
[4,] 13824.2667 -79906.7568
[5,] 1414.9119 13824.2667
[6,] 47341.0535 1414.9119
[7,] -10042.9669 47341.0535
[8,] 9242.0284 -10042.9669
[9,] -4132.8897 9242.0284
[10,] 6500.4613 -4132.8897
[11,] 37402.4152 6500.4613
[12,] -10629.4589 37402.4152
[13,] -7531.4032 -10629.4589
[14,] 10261.3963 -7531.4032
[15,] 23727.6680 10261.3963
[16,] -54967.3791 23727.6680
[17,] 27755.9425 -54967.3791
[18,] 15844.4991 27755.9425
[19,] -17906.3449 15844.4991
[20,] -9958.8697 -17906.3449
[21,] 756.1881 -9958.8697
[22,] 112377.1691 756.1881
[23,] 28142.5729 112377.1691
[24,] -33255.9462 28142.5729
[25,] -102330.0619 -33255.9462
[26,] -34581.0590 -102330.0619
[27,] -628.9763 -34581.0590
[28,] 23073.6909 -628.9763
[29,] -9933.2444 23073.6909
[30,] 965.7411 -9933.2444
[31,] 5931.8636 965.7411
[32,] 16177.9964 5931.8636
[33,] 31406.7090 16177.9964
[34,] -17435.2246 31406.7090
[35,] 34370.6320 -17435.2246
[36,] 19041.4965 34370.6320
[37,] 10621.6218 19041.4965
[38,] 86206.2370 10621.6218
[39,] -15061.1591 86206.2370
[40,] 13370.7126 -15061.1591
[41,] 19643.8737 13370.7126
[42,] -38515.3031 19643.8737
[43,] 18099.7050 -38515.3031
[44,] -14519.1521 18099.7050
[45,] -8454.7842 -14519.1521
[46,] -26333.6942 -8454.7842
[47,] -2475.7247 -26333.6942
[48,] -889.4755 -2475.7247
[49,] -74800.6493 -889.4755
[50,] 22206.1124 -74800.6493
[51,] -15730.6996 22206.1124
[52,] -14849.2989 -15730.6996
[53,] -15893.9604 -14849.2989
[54,] -55569.1132 -15893.9604
[55,] 6316.5078 -55569.1132
[56,] -27467.1603 6316.5078
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -4590.3153 51864.5110
2 14503.0877 -4590.3153
3 -79906.7568 14503.0877
4 13824.2667 -79906.7568
5 1414.9119 13824.2667
6 47341.0535 1414.9119
7 -10042.9669 47341.0535
8 9242.0284 -10042.9669
9 -4132.8897 9242.0284
10 6500.4613 -4132.8897
11 37402.4152 6500.4613
12 -10629.4589 37402.4152
13 -7531.4032 -10629.4589
14 10261.3963 -7531.4032
15 23727.6680 10261.3963
16 -54967.3791 23727.6680
17 27755.9425 -54967.3791
18 15844.4991 27755.9425
19 -17906.3449 15844.4991
20 -9958.8697 -17906.3449
21 756.1881 -9958.8697
22 112377.1691 756.1881
23 28142.5729 112377.1691
24 -33255.9462 28142.5729
25 -102330.0619 -33255.9462
26 -34581.0590 -102330.0619
27 -628.9763 -34581.0590
28 23073.6909 -628.9763
29 -9933.2444 23073.6909
30 965.7411 -9933.2444
31 5931.8636 965.7411
32 16177.9964 5931.8636
33 31406.7090 16177.9964
34 -17435.2246 31406.7090
35 34370.6320 -17435.2246
36 19041.4965 34370.6320
37 10621.6218 19041.4965
38 86206.2370 10621.6218
39 -15061.1591 86206.2370
40 13370.7126 -15061.1591
41 19643.8737 13370.7126
42 -38515.3031 19643.8737
43 18099.7050 -38515.3031
44 -14519.1521 18099.7050
45 -8454.7842 -14519.1521
46 -26333.6942 -8454.7842
47 -2475.7247 -26333.6942
48 -889.4755 -2475.7247
49 -74800.6493 -889.4755
50 22206.1124 -74800.6493
51 -15730.6996 22206.1124
52 -14849.2989 -15730.6996
53 -15893.9604 -14849.2989
54 -55569.1132 -15893.9604
55 6316.5078 -55569.1132
56 -27467.1603 6316.5078
> 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/fisher/rcomp/tmp/76pzq1355155700.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/fisher/rcomp/tmp/8omk61355155700.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/fisher/rcomp/tmp/96a531355155700.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/fisher/rcomp/tmp/101trv1355155700.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/11s91f1355155700.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/fisher/rcomp/tmp/12jdht1355155700.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/fisher/rcomp/tmp/13tzqc1355155701.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/fisher/rcomp/tmp/14c0rm1355155701.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/fisher/rcomp/tmp/15tvch1355155701.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/fisher/rcomp/tmp/16tcty1355155701.tab")
+ }
>
> try(system("convert tmp/1yzdu1355155700.ps tmp/1yzdu1355155700.png",intern=TRUE))
character(0)
> try(system("convert tmp/20tdc1355155700.ps tmp/20tdc1355155700.png",intern=TRUE))
character(0)
> try(system("convert tmp/3b6kl1355155700.ps tmp/3b6kl1355155700.png",intern=TRUE))
character(0)
> try(system("convert tmp/4wram1355155700.ps tmp/4wram1355155700.png",intern=TRUE))
character(0)
> try(system("convert tmp/54nre1355155700.ps tmp/54nre1355155700.png",intern=TRUE))
character(0)
> try(system("convert tmp/6uysl1355155700.ps tmp/6uysl1355155700.png",intern=TRUE))
character(0)
> try(system("convert tmp/76pzq1355155700.ps tmp/76pzq1355155700.png",intern=TRUE))
character(0)
> try(system("convert tmp/8omk61355155700.ps tmp/8omk61355155700.png",intern=TRUE))
character(0)
> try(system("convert tmp/96a531355155700.ps tmp/96a531355155700.png",intern=TRUE))
character(0)
> try(system("convert tmp/101trv1355155700.ps tmp/101trv1355155700.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
5.898 1.521 7.417