R version 2.12.1 (2010-12-16)
Copyright (C) 2010 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(27.25111111
+ ,54
+ ,17
+ ,18
+ ,-3
+ ,32.94777778
+ ,46
+ ,12
+ ,19
+ ,-2
+ ,30.12388889
+ ,60
+ ,17
+ ,19
+ ,0
+ ,27.26277778
+ ,40
+ ,17
+ ,19
+ ,0
+ ,23.3625
+ ,20
+ ,17
+ ,19
+ ,0
+ ,36.27361111
+ ,46
+ ,29
+ ,20
+ ,-4
+ ,18.1875
+ ,18
+ ,13
+ ,20
+ ,1
+ ,33.84666667
+ ,39
+ ,17
+ ,20
+ ,2
+ ,41.82777778
+ ,77
+ ,22
+ ,20
+ ,-4
+ ,62.31388889
+ ,83
+ ,39
+ ,20
+ ,0
+ ,94.88055556
+ ,168
+ ,21
+ ,20
+ ,0
+ ,37.03555556
+ ,55
+ ,20
+ ,20
+ ,-3
+ ,28.20083333
+ ,42
+ ,22
+ ,20
+ ,0
+ ,41.42
+ ,56
+ ,35
+ ,21
+ ,-4
+ ,8.608055556
+ ,14
+ ,17
+ ,21
+ ,-2
+ ,90.22194444
+ ,154
+ ,47
+ ,21
+ ,0
+ ,64.15666667
+ ,53
+ ,30
+ ,21
+ ,-1
+ ,54.59805556
+ ,57
+ ,29
+ ,21
+ ,-3
+ ,39.93222222
+ ,46
+ ,34
+ ,21
+ ,4
+ ,42.30527778
+ ,53
+ ,33
+ ,21
+ ,2
+ ,62.51666667
+ ,93
+ ,41
+ ,21
+ ,1
+ ,42.35388889
+ ,65
+ ,32
+ ,21
+ ,0
+ ,27.1775
+ ,38
+ ,24
+ ,21
+ ,-1
+ ,54.39944444
+ ,67
+ ,31
+ ,21
+ ,-2
+ ,70.69111111
+ ,83
+ ,39
+ ,21
+ ,0
+ ,25.69416667
+ ,32
+ ,18
+ ,21
+ ,-2
+ ,8.826111111
+ ,23
+ ,17
+ ,21
+ ,1
+ ,58.58527778
+ ,56
+ ,30
+ ,21
+ ,2
+ ,41.40583333
+ ,44
+ ,26
+ ,21
+ ,0
+ ,65.8925
+ ,84
+ ,38
+ ,21
+ ,0
+ ,36.98083333
+ ,55
+ ,30
+ ,21
+ ,4
+ ,48.44861111
+ ,100
+ ,31
+ ,21
+ ,3
+ ,81.78444444
+ ,77
+ ,33
+ ,21
+ ,4
+ ,90.3075
+ ,99
+ ,36
+ ,21
+ ,3
+ ,29.55777778
+ ,30
+ ,14
+ ,21
+ ,1
+ ,73.82472222
+ ,146
+ ,32
+ ,21
+ ,-2
+ ,100.6391667
+ ,119
+ ,34
+ ,21
+ ,2
+ ,60.81833333
+ ,41
+ ,29
+ ,21
+ ,2
+ ,31.28083333
+ ,41
+ ,20
+ ,21
+ ,3
+ ,41.235
+ ,91
+ ,37
+ ,21
+ ,3
+ ,50.5775
+ ,63
+ ,33
+ ,21
+ ,2
+ ,36.80194444
+ ,41
+ ,36
+ ,21
+ ,3
+ ,35.67305556
+ ,64
+ ,38
+ ,21
+ ,2
+ ,47.915
+ ,52
+ ,43
+ ,21
+ ,3
+ ,90.16611111
+ ,110
+ ,37
+ ,21
+ ,2
+ ,50.45361111
+ ,70
+ ,30
+ ,21
+ ,6
+ ,21.195
+ ,31
+ ,20
+ ,21
+ ,2
+ ,16.495
+ ,49
+ ,12
+ ,21
+ ,1
+ ,67.79222222
+ ,68
+ ,44
+ ,22
+ ,2
+ ,46.52444444
+ ,45
+ ,28
+ ,22
+ ,0
+ ,95.63805556
+ ,75
+ ,30
+ ,22
+ ,1
+ ,45.2125
+ ,32
+ ,28
+ ,22
+ ,-3
+ ,23.77055556
+ ,34
+ ,21
+ ,22
+ ,0
+ ,88.165
+ ,86
+ ,31
+ ,22
+ ,-3
+ ,39.79055556
+ ,103
+ ,27
+ ,22
+ ,2
+ ,53.70527778
+ ,78
+ ,35
+ ,22
+ ,2
+ ,67.98583333
+ ,95
+ ,33
+ ,22
+ ,0
+ ,63.67833333
+ ,247
+ ,31
+ ,22
+ ,2
+ ,65.77361111
+ ,119
+ ,31
+ ,23
+ ,0
+ ,67.64194444
+ ,71
+ ,42
+ ,23
+ ,0
+ ,62.12
+ ,73
+ ,33
+ ,23
+ ,-1
+ ,27.98611111
+ ,72
+ ,30
+ ,23
+ ,3
+ ,26.45194444
+ ,34
+ ,32
+ ,23
+ ,3
+ ,50.87972222
+ ,66
+ ,39
+ ,24
+ ,-4
+ ,67.51666667
+ ,63
+ ,29
+ ,24
+ ,-1
+ ,42.46416667
+ ,58
+ ,28
+ ,24
+ ,2
+ ,26.82222222
+ ,76
+ ,17
+ ,25
+ ,0
+ ,75.51555556
+ ,103
+ ,37
+ ,25
+ ,-3
+ ,48.53444444
+ ,92
+ ,34
+ ,26
+ ,0)
+ ,dim=c(5
+ ,69)
+ ,dimnames=list(c('time_in_rfc'
+ ,'logins'
+ ,'compendiums_reviewed'
+ ,'What_is_your_age?'
+ ,'Totale_score')
+ ,1:69))
> y <- array(NA,dim=c(5,69),dimnames=list(c('time_in_rfc','logins','compendiums_reviewed','What_is_your_age?','Totale_score'),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 = '5'
> 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
Totale_score time_in_rfc logins compendiums_reviewed What_is_your_age?
1 -3 27.251111 54 17 18
2 -2 32.947778 46 12 19
3 0 30.123889 60 17 19
4 0 27.262778 40 17 19
5 0 23.362500 20 17 19
6 -4 36.273611 46 29 20
7 1 18.187500 18 13 20
8 2 33.846667 39 17 20
9 -4 41.827778 77 22 20
10 0 62.313889 83 39 20
11 0 94.880556 168 21 20
12 -3 37.035556 55 20 20
13 0 28.200833 42 22 20
14 -4 41.420000 56 35 21
15 -2 8.608056 14 17 21
16 0 90.221944 154 47 21
17 -1 64.156667 53 30 21
18 -3 54.598056 57 29 21
19 4 39.932222 46 34 21
20 2 42.305278 53 33 21
21 1 62.516667 93 41 21
22 0 42.353889 65 32 21
23 -1 27.177500 38 24 21
24 -2 54.399444 67 31 21
25 0 70.691111 83 39 21
26 -2 25.694167 32 18 21
27 1 8.826111 23 17 21
28 2 58.585278 56 30 21
29 0 41.405833 44 26 21
30 0 65.892500 84 38 21
31 4 36.980833 55 30 21
32 3 48.448611 100 31 21
33 4 81.784444 77 33 21
34 3 90.307500 99 36 21
35 1 29.557778 30 14 21
36 -2 73.824722 146 32 21
37 2 100.639167 119 34 21
38 2 60.818333 41 29 21
39 3 31.280833 41 20 21
40 3 41.235000 91 37 21
41 2 50.577500 63 33 21
42 3 36.801944 41 36 21
43 2 35.673056 64 38 21
44 3 47.915000 52 43 21
45 2 90.166111 110 37 21
46 6 50.453611 70 30 21
47 2 21.195000 31 20 21
48 1 16.495000 49 12 21
49 2 67.792222 68 44 22
50 0 46.524444 45 28 22
51 1 95.638056 75 30 22
52 -3 45.212500 32 28 22
53 0 23.770556 34 21 22
54 -3 88.165000 86 31 22
55 2 39.790556 103 27 22
56 2 53.705278 78 35 22
57 0 67.985833 95 33 22
58 2 63.678333 247 31 22
59 0 65.773611 119 31 23
60 0 67.641944 71 42 23
61 -1 62.120000 73 33 23
62 3 27.986111 72 30 23
63 3 26.451944 34 32 23
64 -4 50.879722 66 39 24
65 -1 67.516667 63 29 24
66 2 42.464167 58 28 24
67 0 26.822222 76 17 25
68 -3 75.515556 103 37 25
69 0 48.534444 92 34 26
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) time_in_rfc logins
0.980345 -0.020724 0.005173
compendiums_reviewed `What_is_your_age?`
0.077626 -0.100730
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.0132 -1.3155 -0.0281 1.5624 5.4897
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.980345 4.167887 0.235 0.8148
time_in_rfc -0.020724 0.019484 -1.064 0.2915
logins 0.005173 0.009789 0.528 0.5990
compendiums_reviewed 0.077626 0.043863 1.770 0.0815 .
`What_is_your_age?` -0.100730 0.207086 -0.486 0.6283
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.274 on 64 degrees of freedom
Multiple R-squared: 0.05095, Adjusted R-squared: -0.008365
F-statistic: 0.859 on 4 and 64 DF, p-value: 0.4935
> 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.5056880 0.9886241 0.49431204
[2,] 0.4596899 0.9193797 0.54031014
[3,] 0.6274743 0.7450515 0.37252575
[4,] 0.4982067 0.9964133 0.50179334
[5,] 0.5084418 0.9831163 0.49155816
[6,] 0.4725651 0.9451302 0.52743488
[7,] 0.5077655 0.9844690 0.49223452
[8,] 0.4496697 0.8993394 0.55033031
[9,] 0.5446656 0.9106688 0.45533439
[10,] 0.5353800 0.9292400 0.46462002
[11,] 0.6017912 0.7964176 0.39820880
[12,] 0.8800425 0.2399150 0.11995748
[13,] 0.8821398 0.2357205 0.11786023
[14,] 0.8541241 0.2917519 0.14587593
[15,] 0.8186168 0.3627664 0.18138319
[16,] 0.7953388 0.4093224 0.20466121
[17,] 0.8207394 0.3585211 0.17926055
[18,] 0.7843191 0.4313619 0.21568095
[19,] 0.8098789 0.3802422 0.19012109
[20,] 0.7949526 0.4100949 0.20504744
[21,] 0.7777520 0.4444960 0.22224799
[22,] 0.7456509 0.5086981 0.25434905
[23,] 0.7156024 0.5687953 0.28439763
[24,] 0.8038561 0.3922879 0.19614394
[25,] 0.8025773 0.3948453 0.19742267
[26,] 0.8729564 0.2540872 0.12704362
[27,] 0.8736613 0.2526773 0.12633866
[28,] 0.8351223 0.3297554 0.16487771
[29,] 0.8928095 0.2143811 0.10719055
[30,] 0.8696195 0.2607610 0.13038049
[31,] 0.8339945 0.3320111 0.16600553
[32,] 0.8281998 0.3436005 0.17180024
[33,] 0.8243153 0.3513693 0.17568466
[34,] 0.7778740 0.4442520 0.22212599
[35,] 0.7419374 0.5161252 0.25806260
[36,] 0.6984494 0.6031013 0.30155063
[37,] 0.6394509 0.7210982 0.36054909
[38,] 0.5881506 0.8236989 0.41184944
[39,] 0.8547901 0.2904198 0.14520991
[40,] 0.8089469 0.3821063 0.19105313
[41,] 0.7547634 0.4904732 0.24523658
[42,] 0.7367227 0.5265545 0.26327727
[43,] 0.6755183 0.6489634 0.32448169
[44,] 0.8007984 0.3984032 0.19920159
[45,] 0.9042798 0.1914404 0.09572018
[46,] 0.9324155 0.1351689 0.06758447
[47,] 0.9254633 0.1490734 0.07453670
[48,] 0.8867645 0.2264709 0.11323546
[49,] 0.8307530 0.3384939 0.16924696
[50,] 0.7461218 0.5077563 0.25387816
[51,] 0.6340840 0.7318320 0.36591601
[52,] 0.5074396 0.9851209 0.49256043
[53,] 0.4145812 0.8291624 0.58541882
[54,] 0.2733852 0.5467705 0.72661475
> postscript(file="/var/www/rcomp/tmp/1la6b1323935298.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/www/rcomp/tmp/2tt5f1323935298.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/www/rcomp/tmp/39cy01323935298.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/www/rcomp/tmp/49b161323935298.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/www/rcomp/tmp/5h99k1323935298.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
-3.201456985 -1.553154560 -0.072231733 -0.028059683 -0.005423318 -4.703138806
7 8 9 10 11 12
1.308909242 2.214287615 -4.205025633 -1.131149603 0.501296720 -3.035275782
13 14 15 16 17 18
-0.306365386 -5.013243094 -2.078695587 -1.440363181 -1.138399743 -3.279559712
19 20 21 22 23 24
3.085282842 1.175874822 -0.233202128 -0.807571399 -1.361401839 -2.490660148
25 26 27 28 29 30
-0.856810559 -1.895348016 0.879263829 1.730618964 -0.252825205 -0.883804455
31 32 33 34 35 36
3.288061842 2.215296399 3.869882046 2.699824183 1.505571149 -2.574406095
37 38 39 40 41 42
1.965723461 1.932121859 3.018618903 1.646606417 1.295575466 1.891025838
43 44 45 46 47 48
0.593393622 1.521046124 1.562362119 5.489672349 1.861332991 1.291817353
49 50 51 52 53 54
0.873312767 -0.206442456 1.500937670 -3.166378479 -0.077707318 -2.788465820
55 56 57 58 59 60
1.431579706 1.228274201 -0.408469912 0.871173698 -0.322493742 -0.889340347
61 62 63 64 65 66
-1.315491611 2.215168471 2.224707677 -4.877247255 -0.740685598 1.843618717
67 68 69
0.380949154 -3.302124813 -0.470768217
> postscript(file="/var/www/rcomp/tmp/642c71323935298.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 -3.201456985 NA
1 -1.553154560 -3.201456985
2 -0.072231733 -1.553154560
3 -0.028059683 -0.072231733
4 -0.005423318 -0.028059683
5 -4.703138806 -0.005423318
6 1.308909242 -4.703138806
7 2.214287615 1.308909242
8 -4.205025633 2.214287615
9 -1.131149603 -4.205025633
10 0.501296720 -1.131149603
11 -3.035275782 0.501296720
12 -0.306365386 -3.035275782
13 -5.013243094 -0.306365386
14 -2.078695587 -5.013243094
15 -1.440363181 -2.078695587
16 -1.138399743 -1.440363181
17 -3.279559712 -1.138399743
18 3.085282842 -3.279559712
19 1.175874822 3.085282842
20 -0.233202128 1.175874822
21 -0.807571399 -0.233202128
22 -1.361401839 -0.807571399
23 -2.490660148 -1.361401839
24 -0.856810559 -2.490660148
25 -1.895348016 -0.856810559
26 0.879263829 -1.895348016
27 1.730618964 0.879263829
28 -0.252825205 1.730618964
29 -0.883804455 -0.252825205
30 3.288061842 -0.883804455
31 2.215296399 3.288061842
32 3.869882046 2.215296399
33 2.699824183 3.869882046
34 1.505571149 2.699824183
35 -2.574406095 1.505571149
36 1.965723461 -2.574406095
37 1.932121859 1.965723461
38 3.018618903 1.932121859
39 1.646606417 3.018618903
40 1.295575466 1.646606417
41 1.891025838 1.295575466
42 0.593393622 1.891025838
43 1.521046124 0.593393622
44 1.562362119 1.521046124
45 5.489672349 1.562362119
46 1.861332991 5.489672349
47 1.291817353 1.861332991
48 0.873312767 1.291817353
49 -0.206442456 0.873312767
50 1.500937670 -0.206442456
51 -3.166378479 1.500937670
52 -0.077707318 -3.166378479
53 -2.788465820 -0.077707318
54 1.431579706 -2.788465820
55 1.228274201 1.431579706
56 -0.408469912 1.228274201
57 0.871173698 -0.408469912
58 -0.322493742 0.871173698
59 -0.889340347 -0.322493742
60 -1.315491611 -0.889340347
61 2.215168471 -1.315491611
62 2.224707677 2.215168471
63 -4.877247255 2.224707677
64 -0.740685598 -4.877247255
65 1.843618717 -0.740685598
66 0.380949154 1.843618717
67 -3.302124813 0.380949154
68 -0.470768217 -3.302124813
69 NA -0.470768217
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.553154560 -3.201456985
[2,] -0.072231733 -1.553154560
[3,] -0.028059683 -0.072231733
[4,] -0.005423318 -0.028059683
[5,] -4.703138806 -0.005423318
[6,] 1.308909242 -4.703138806
[7,] 2.214287615 1.308909242
[8,] -4.205025633 2.214287615
[9,] -1.131149603 -4.205025633
[10,] 0.501296720 -1.131149603
[11,] -3.035275782 0.501296720
[12,] -0.306365386 -3.035275782
[13,] -5.013243094 -0.306365386
[14,] -2.078695587 -5.013243094
[15,] -1.440363181 -2.078695587
[16,] -1.138399743 -1.440363181
[17,] -3.279559712 -1.138399743
[18,] 3.085282842 -3.279559712
[19,] 1.175874822 3.085282842
[20,] -0.233202128 1.175874822
[21,] -0.807571399 -0.233202128
[22,] -1.361401839 -0.807571399
[23,] -2.490660148 -1.361401839
[24,] -0.856810559 -2.490660148
[25,] -1.895348016 -0.856810559
[26,] 0.879263829 -1.895348016
[27,] 1.730618964 0.879263829
[28,] -0.252825205 1.730618964
[29,] -0.883804455 -0.252825205
[30,] 3.288061842 -0.883804455
[31,] 2.215296399 3.288061842
[32,] 3.869882046 2.215296399
[33,] 2.699824183 3.869882046
[34,] 1.505571149 2.699824183
[35,] -2.574406095 1.505571149
[36,] 1.965723461 -2.574406095
[37,] 1.932121859 1.965723461
[38,] 3.018618903 1.932121859
[39,] 1.646606417 3.018618903
[40,] 1.295575466 1.646606417
[41,] 1.891025838 1.295575466
[42,] 0.593393622 1.891025838
[43,] 1.521046124 0.593393622
[44,] 1.562362119 1.521046124
[45,] 5.489672349 1.562362119
[46,] 1.861332991 5.489672349
[47,] 1.291817353 1.861332991
[48,] 0.873312767 1.291817353
[49,] -0.206442456 0.873312767
[50,] 1.500937670 -0.206442456
[51,] -3.166378479 1.500937670
[52,] -0.077707318 -3.166378479
[53,] -2.788465820 -0.077707318
[54,] 1.431579706 -2.788465820
[55,] 1.228274201 1.431579706
[56,] -0.408469912 1.228274201
[57,] 0.871173698 -0.408469912
[58,] -0.322493742 0.871173698
[59,] -0.889340347 -0.322493742
[60,] -1.315491611 -0.889340347
[61,] 2.215168471 -1.315491611
[62,] 2.224707677 2.215168471
[63,] -4.877247255 2.224707677
[64,] -0.740685598 -4.877247255
[65,] 1.843618717 -0.740685598
[66,] 0.380949154 1.843618717
[67,] -3.302124813 0.380949154
[68,] -0.470768217 -3.302124813
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.553154560 -3.201456985
2 -0.072231733 -1.553154560
3 -0.028059683 -0.072231733
4 -0.005423318 -0.028059683
5 -4.703138806 -0.005423318
6 1.308909242 -4.703138806
7 2.214287615 1.308909242
8 -4.205025633 2.214287615
9 -1.131149603 -4.205025633
10 0.501296720 -1.131149603
11 -3.035275782 0.501296720
12 -0.306365386 -3.035275782
13 -5.013243094 -0.306365386
14 -2.078695587 -5.013243094
15 -1.440363181 -2.078695587
16 -1.138399743 -1.440363181
17 -3.279559712 -1.138399743
18 3.085282842 -3.279559712
19 1.175874822 3.085282842
20 -0.233202128 1.175874822
21 -0.807571399 -0.233202128
22 -1.361401839 -0.807571399
23 -2.490660148 -1.361401839
24 -0.856810559 -2.490660148
25 -1.895348016 -0.856810559
26 0.879263829 -1.895348016
27 1.730618964 0.879263829
28 -0.252825205 1.730618964
29 -0.883804455 -0.252825205
30 3.288061842 -0.883804455
31 2.215296399 3.288061842
32 3.869882046 2.215296399
33 2.699824183 3.869882046
34 1.505571149 2.699824183
35 -2.574406095 1.505571149
36 1.965723461 -2.574406095
37 1.932121859 1.965723461
38 3.018618903 1.932121859
39 1.646606417 3.018618903
40 1.295575466 1.646606417
41 1.891025838 1.295575466
42 0.593393622 1.891025838
43 1.521046124 0.593393622
44 1.562362119 1.521046124
45 5.489672349 1.562362119
46 1.861332991 5.489672349
47 1.291817353 1.861332991
48 0.873312767 1.291817353
49 -0.206442456 0.873312767
50 1.500937670 -0.206442456
51 -3.166378479 1.500937670
52 -0.077707318 -3.166378479
53 -2.788465820 -0.077707318
54 1.431579706 -2.788465820
55 1.228274201 1.431579706
56 -0.408469912 1.228274201
57 0.871173698 -0.408469912
58 -0.322493742 0.871173698
59 -0.889340347 -0.322493742
60 -1.315491611 -0.889340347
61 2.215168471 -1.315491611
62 2.224707677 2.215168471
63 -4.877247255 2.224707677
64 -0.740685598 -4.877247255
65 1.843618717 -0.740685598
66 0.380949154 1.843618717
67 -3.302124813 0.380949154
68 -0.470768217 -3.302124813
> 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/www/rcomp/tmp/7n1jy1323935298.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/www/rcomp/tmp/8otza1323935298.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/www/rcomp/tmp/935mx1323935298.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/www/rcomp/tmp/10zrc81323935298.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/www/rcomp/tmp/113b9c1323935298.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/www/rcomp/tmp/12mbj41323935298.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/www/rcomp/tmp/13ktpp1323935298.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/www/rcomp/tmp/14v1wj1323935298.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/www/rcomp/tmp/152gir1323935298.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/www/rcomp/tmp/169h2d1323935298.tab")
+ }
>
> try(system("convert tmp/1la6b1323935298.ps tmp/1la6b1323935298.png",intern=TRUE))
character(0)
> try(system("convert tmp/2tt5f1323935298.ps tmp/2tt5f1323935298.png",intern=TRUE))
character(0)
> try(system("convert tmp/39cy01323935298.ps tmp/39cy01323935298.png",intern=TRUE))
character(0)
> try(system("convert tmp/49b161323935298.ps tmp/49b161323935298.png",intern=TRUE))
character(0)
> try(system("convert tmp/5h99k1323935298.ps tmp/5h99k1323935298.png",intern=TRUE))
character(0)
> try(system("convert tmp/642c71323935298.ps tmp/642c71323935298.png",intern=TRUE))
character(0)
> try(system("convert tmp/7n1jy1323935298.ps tmp/7n1jy1323935298.png",intern=TRUE))
character(0)
> try(system("convert tmp/8otza1323935298.ps tmp/8otza1323935298.png",intern=TRUE))
character(0)
> try(system("convert tmp/935mx1323935298.ps tmp/935mx1323935298.png",intern=TRUE))
character(0)
> try(system("convert tmp/10zrc81323935298.ps tmp/10zrc81323935298.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.852 0.680 4.623