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(27.25111111
+ ,54
+ ,17
+ ,1
+ ,18
+ ,-3
+ ,32.94777778
+ ,46
+ ,12
+ ,1
+ ,19
+ ,-2
+ ,30.12388889
+ ,60
+ ,17
+ ,1
+ ,19
+ ,0
+ ,27.26277778
+ ,40
+ ,17
+ ,1
+ ,19
+ ,0
+ ,23.3625
+ ,20
+ ,17
+ ,1
+ ,19
+ ,0
+ ,36.27361111
+ ,46
+ ,29
+ ,1
+ ,20
+ ,-4
+ ,18.1875
+ ,18
+ ,13
+ ,0
+ ,20
+ ,1
+ ,33.84666667
+ ,39
+ ,17
+ ,1
+ ,20
+ ,2
+ ,41.82777778
+ ,77
+ ,22
+ ,0
+ ,20
+ ,-4
+ ,62.31388889
+ ,83
+ ,39
+ ,1
+ ,20
+ ,0
+ ,94.88055556
+ ,168
+ ,21
+ ,0
+ ,20
+ ,0
+ ,37.03555556
+ ,55
+ ,20
+ ,0
+ ,20
+ ,-3
+ ,28.20083333
+ ,42
+ ,22
+ ,0
+ ,20
+ ,0
+ ,41.42
+ ,56
+ ,35
+ ,1
+ ,21
+ ,-4
+ ,8.608055556
+ ,14
+ ,17
+ ,0
+ ,21
+ ,-2
+ ,90.22194444
+ ,154
+ ,47
+ ,0
+ ,21
+ ,0
+ ,64.15666667
+ ,53
+ ,30
+ ,1
+ ,21
+ ,-1
+ ,54.59805556
+ ,57
+ ,29
+ ,0
+ ,21
+ ,-3
+ ,39.93222222
+ ,46
+ ,34
+ ,1
+ ,21
+ ,4
+ ,42.30527778
+ ,53
+ ,33
+ ,0
+ ,21
+ ,2
+ ,62.51666667
+ ,93
+ ,41
+ ,1
+ ,21
+ ,1
+ ,42.35388889
+ ,65
+ ,32
+ ,0
+ ,21
+ ,0
+ ,27.1775
+ ,38
+ ,24
+ ,1
+ ,21
+ ,-1
+ ,54.39944444
+ ,67
+ ,31
+ ,1
+ ,21
+ ,-2
+ ,70.69111111
+ ,83
+ ,39
+ ,1
+ ,21
+ ,0
+ ,25.69416667
+ ,32
+ ,18
+ ,0
+ ,21
+ ,-2
+ ,8.826111111
+ ,23
+ ,17
+ ,1
+ ,21
+ ,1
+ ,58.58527778
+ ,56
+ ,30
+ ,1
+ ,21
+ ,2
+ ,41.40583333
+ ,44
+ ,26
+ ,1
+ ,21
+ ,0
+ ,65.8925
+ ,84
+ ,38
+ ,0
+ ,21
+ ,0
+ ,36.98083333
+ ,55
+ ,30
+ ,0
+ ,21
+ ,4
+ ,48.44861111
+ ,100
+ ,31
+ ,0
+ ,21
+ ,3
+ ,81.78444444
+ ,77
+ ,33
+ ,1
+ ,21
+ ,4
+ ,90.3075
+ ,99
+ ,36
+ ,0
+ ,21
+ ,3
+ ,29.55777778
+ ,30
+ ,14
+ ,1
+ ,21
+ ,1
+ ,73.82472222
+ ,146
+ ,32
+ ,0
+ ,21
+ ,-2
+ ,100.6391667
+ ,119
+ ,34
+ ,0
+ ,21
+ ,2
+ ,60.81833333
+ ,41
+ ,29
+ ,0
+ ,21
+ ,2
+ ,31.28083333
+ ,41
+ ,20
+ ,0
+ ,21
+ ,3
+ ,41.235
+ ,91
+ ,37
+ ,0
+ ,21
+ ,3
+ ,50.5775
+ ,63
+ ,33
+ ,1
+ ,21
+ ,2
+ ,36.80194444
+ ,41
+ ,36
+ ,0
+ ,21
+ ,3
+ ,35.67305556
+ ,64
+ ,38
+ ,1
+ ,21
+ ,2
+ ,47.915
+ ,52
+ ,43
+ ,0
+ ,21
+ ,3
+ ,90.16611111
+ ,110
+ ,37
+ ,1
+ ,21
+ ,2
+ ,50.45361111
+ ,70
+ ,30
+ ,0
+ ,21
+ ,6
+ ,21.195
+ ,31
+ ,20
+ ,0
+ ,21
+ ,2
+ ,16.495
+ ,49
+ ,12
+ ,1
+ ,21
+ ,1
+ ,67.79222222
+ ,68
+ ,44
+ ,0
+ ,22
+ ,2
+ ,46.52444444
+ ,45
+ ,28
+ ,1
+ ,22
+ ,0
+ ,95.63805556
+ ,75
+ ,30
+ ,0
+ ,22
+ ,1
+ ,45.2125
+ ,32
+ ,28
+ ,0
+ ,22
+ ,-3
+ ,23.77055556
+ ,34
+ ,21
+ ,0
+ ,22
+ ,0
+ ,88.165
+ ,86
+ ,31
+ ,1
+ ,22
+ ,-3
+ ,39.79055556
+ ,103
+ ,27
+ ,1
+ ,22
+ ,2
+ ,53.70527778
+ ,78
+ ,35
+ ,0
+ ,22
+ ,2
+ ,67.98583333
+ ,95
+ ,33
+ ,0
+ ,22
+ ,0
+ ,63.67833333
+ ,247
+ ,31
+ ,0
+ ,22
+ ,2
+ ,65.77361111
+ ,119
+ ,31
+ ,0
+ ,23
+ ,0
+ ,67.64194444
+ ,71
+ ,42
+ ,0
+ ,23
+ ,0
+ ,62.12
+ ,73
+ ,33
+ ,0
+ ,23
+ ,-1
+ ,27.98611111
+ ,72
+ ,30
+ ,0
+ ,23
+ ,3
+ ,26.45194444
+ ,34
+ ,32
+ ,0
+ ,23
+ ,3
+ ,50.87972222
+ ,66
+ ,39
+ ,1
+ ,24
+ ,-4
+ ,67.51666667
+ ,63
+ ,29
+ ,0
+ ,24
+ ,-1
+ ,42.46416667
+ ,58
+ ,28
+ ,1
+ ,24
+ ,2
+ ,26.82222222
+ ,76
+ ,17
+ ,1
+ ,25
+ ,0
+ ,75.51555556
+ ,103
+ ,37
+ ,0
+ ,25
+ ,-3
+ ,48.53444444
+ ,92
+ ,34
+ ,0
+ ,26
+ ,0)
+ ,dim=c(6
+ ,69)
+ ,dimnames=list(c('time_in_rfc'
+ ,'logins'
+ ,'compendiums_reviewed'
+ ,'What_is_your_gender?'
+ ,'What_is_your_age?'
+ ,'Totale_score')
+ ,1:69))
> y <- array(NA,dim=c(6,69),dimnames=list(c('time_in_rfc','logins','compendiums_reviewed','What_is_your_gender?','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 = '6'
> 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_gender?
1 -3 27.251111 54 17 1
2 -2 32.947778 46 12 1
3 0 30.123889 60 17 1
4 0 27.262778 40 17 1
5 0 23.362500 20 17 1
6 -4 36.273611 46 29 1
7 1 18.187500 18 13 0
8 2 33.846667 39 17 1
9 -4 41.827778 77 22 0
10 0 62.313889 83 39 1
11 0 94.880556 168 21 0
12 -3 37.035556 55 20 0
13 0 28.200833 42 22 0
14 -4 41.420000 56 35 1
15 -2 8.608056 14 17 0
16 0 90.221944 154 47 0
17 -1 64.156667 53 30 1
18 -3 54.598056 57 29 0
19 4 39.932222 46 34 1
20 2 42.305278 53 33 0
21 1 62.516667 93 41 1
22 0 42.353889 65 32 0
23 -1 27.177500 38 24 1
24 -2 54.399444 67 31 1
25 0 70.691111 83 39 1
26 -2 25.694167 32 18 0
27 1 8.826111 23 17 1
28 2 58.585278 56 30 1
29 0 41.405833 44 26 1
30 0 65.892500 84 38 0
31 4 36.980833 55 30 0
32 3 48.448611 100 31 0
33 4 81.784444 77 33 1
34 3 90.307500 99 36 0
35 1 29.557778 30 14 1
36 -2 73.824722 146 32 0
37 2 100.639167 119 34 0
38 2 60.818333 41 29 0
39 3 31.280833 41 20 0
40 3 41.235000 91 37 0
41 2 50.577500 63 33 1
42 3 36.801944 41 36 0
43 2 35.673056 64 38 1
44 3 47.915000 52 43 0
45 2 90.166111 110 37 1
46 6 50.453611 70 30 0
47 2 21.195000 31 20 0
48 1 16.495000 49 12 1
49 2 67.792222 68 44 0
50 0 46.524444 45 28 1
51 1 95.638056 75 30 0
52 -3 45.212500 32 28 0
53 0 23.770556 34 21 0
54 -3 88.165000 86 31 1
55 2 39.790556 103 27 1
56 2 53.705278 78 35 0
57 0 67.985833 95 33 0
58 2 63.678333 247 31 0
59 0 65.773611 119 31 0
60 0 67.641944 71 42 0
61 -1 62.120000 73 33 0
62 3 27.986111 72 30 0
63 3 26.451944 34 32 0
64 -4 50.879722 66 39 1
65 -1 67.516667 63 29 0
66 2 42.464167 58 28 1
67 0 26.822222 76 17 1
68 -3 75.515556 103 37 0
69 0 48.534444 92 34 0
What_is_your_age?
1 18
2 19
3 19
4 19
5 19
6 20
7 20
8 20
9 20
10 20
11 20
12 20
13 20
14 21
15 21
16 21
17 21
18 21
19 21
20 21
21 21
22 21
23 21
24 21
25 21
26 21
27 21
28 21
29 21
30 21
31 21
32 21
33 21
34 21
35 21
36 21
37 21
38 21
39 21
40 21
41 21
42 21
43 21
44 21
45 21
46 21
47 21
48 21
49 22
50 22
51 22
52 22
53 22
54 22
55 22
56 22
57 22
58 22
59 23
60 23
61 23
62 23
63 23
64 24
65 24
66 24
67 25
68 25
69 26
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) time_in_rfc logins
2.022138 -0.020563 0.004021
compendiums_reviewed `What_is_your_gender?` `What_is_your_age?`
0.075462 -0.531116 -0.132478
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4.724 -1.114 0.143 1.635 5.252
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.022138 4.324452 0.468 0.6417
time_in_rfc -0.020563 0.019508 -1.054 0.2959
logins 0.004021 0.009881 0.407 0.6855
compendiums_reviewed 0.075462 0.043980 1.716 0.0911 .
`What_is_your_gender?` -0.531116 0.578352 -0.918 0.3620
`What_is_your_age?` -0.132478 0.210202 -0.630 0.5308
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.277 on 63 degrees of freedom
Multiple R-squared: 0.06349, Adjusted R-squared: -0.01084
F-statistic: 0.8542 on 5 and 63 DF, p-value: 0.5169
> 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.2592762 0.51855250 0.74072375
[2,] 0.7223342 0.55533155 0.27766577
[3,] 0.5956367 0.80872666 0.40436333
[4,] 0.5559049 0.88819016 0.44409508
[5,] 0.5548230 0.89035396 0.44517698
[6,] 0.5885197 0.82296055 0.41148028
[7,] 0.5319760 0.93604792 0.46802396
[8,] 0.6092904 0.78141916 0.39070958
[9,] 0.5938942 0.81221163 0.40610582
[10,] 0.6687339 0.66253212 0.33126606
[11,] 0.8927454 0.21450929 0.10725465
[12,] 0.9080237 0.18395265 0.09197633
[13,] 0.8750874 0.24982528 0.12491264
[14,] 0.8478376 0.30432471 0.15216236
[15,] 0.8249036 0.35019272 0.17509636
[16,] 0.8419887 0.31602266 0.15801133
[17,] 0.8014758 0.39704847 0.19852423
[18,] 0.8349093 0.33018132 0.16509066
[19,] 0.8047151 0.39056990 0.19528495
[20,] 0.7832232 0.43355363 0.21677682
[21,] 0.7416656 0.51666874 0.25833437
[22,] 0.7167164 0.56656712 0.28328356
[23,] 0.8137318 0.37253644 0.18626822
[24,] 0.8097988 0.38040237 0.19020119
[25,] 0.8902216 0.21955681 0.10977841
[26,] 0.8904131 0.21917372 0.10958686
[27,] 0.8522084 0.29558325 0.14779163
[28,] 0.9146073 0.17078531 0.08539265
[29,] 0.8911329 0.21773426 0.10886713
[30,] 0.8555837 0.28883257 0.14441628
[31,] 0.8431984 0.31360316 0.15680158
[32,] 0.8372696 0.32546073 0.16273037
[33,] 0.7953692 0.40926159 0.20463080
[34,] 0.7551306 0.48973887 0.24486944
[35,] 0.7016275 0.59674502 0.29837251
[36,] 0.6382658 0.72346843 0.36173421
[37,] 0.6218335 0.75633293 0.37816646
[38,] 0.8580864 0.28382719 0.14191360
[39,] 0.8082966 0.38340678 0.19170339
[40,] 0.7517479 0.49650426 0.24825213
[41,] 0.7353284 0.52934318 0.26467159
[42,] 0.6666044 0.66679129 0.33339565
[43,] 0.7712800 0.45743993 0.22871997
[44,] 0.9073844 0.18523130 0.09261565
[45,] 0.9789382 0.04212364 0.02106182
[46,] 0.9659521 0.06809588 0.03404794
[47,] 0.9405734 0.11885311 0.05942655
[48,] 0.8990827 0.20183459 0.10091730
[49,] 0.8281536 0.34369272 0.17184636
[50,] 0.7350126 0.52997474 0.26498737
[51,] 0.6122744 0.77545115 0.38772558
[52,] 0.5315571 0.93688582 0.46844291
> postscript(file="/var/wessaorg/rcomp/tmp/119qi1323935115.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/2jzk21323935115.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/3zf8d1323935115.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/4cj7i1323935115.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/5w91y1323935115.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.046014439 -1.386917553 0.121414059 0.142991617 0.143200456 -4.468909368
7 8 9 10 11 12
0.948038193 2.414876920 -4.482214322 -0.836819604 0.318317068 -3.341380957
13 14 15 16 17 18
-0.621709438 -4.723584323 -2.402236056 -1.550734625 -0.866669384 -3.534961881
19 20 21 22 23 24
3.361490135 0.926490628 0.108697802 -1.045294314 -1.114001149 -2.199060712
25 26 27 28 29 30
-0.532077984 -2.198721552 1.097178456 2.006702620 0.003532495 -1.090427337
31 32 33 34 35 36
3.035348281 2.014775643 4.172935415 2.502242278 1.721733461 -2.723816021
37 38 39 40 41 42
1.785208883 1.657276909 2.729048982 1.449850857 1.587504995 1.635183638
43 44 45 46 47 48
0.899687723 1.291242149 1.910761908 5.252084729 1.561856489 1.527653124
49 50 51 52 53 54
0.692670416 0.086320949 1.293602270 -3.419505251 -0.340227522 -2.448641530
55 56 57 58 59 60
1.790118795 1.041951498 -0.581817893 0.869402831 -0.440399535 -1.039078711
61 62 63 64 65 66
-1.481508163 2.046993306 2.017302894 -4.473681909 -0.896001481 2.215516851
67 68 69
0.784060172 -3.363561541 -0.515291344
> postscript(file="/var/wessaorg/rcomp/tmp/63goh1323935115.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.046014439 NA
1 -1.386917553 -3.046014439
2 0.121414059 -1.386917553
3 0.142991617 0.121414059
4 0.143200456 0.142991617
5 -4.468909368 0.143200456
6 0.948038193 -4.468909368
7 2.414876920 0.948038193
8 -4.482214322 2.414876920
9 -0.836819604 -4.482214322
10 0.318317068 -0.836819604
11 -3.341380957 0.318317068
12 -0.621709438 -3.341380957
13 -4.723584323 -0.621709438
14 -2.402236056 -4.723584323
15 -1.550734625 -2.402236056
16 -0.866669384 -1.550734625
17 -3.534961881 -0.866669384
18 3.361490135 -3.534961881
19 0.926490628 3.361490135
20 0.108697802 0.926490628
21 -1.045294314 0.108697802
22 -1.114001149 -1.045294314
23 -2.199060712 -1.114001149
24 -0.532077984 -2.199060712
25 -2.198721552 -0.532077984
26 1.097178456 -2.198721552
27 2.006702620 1.097178456
28 0.003532495 2.006702620
29 -1.090427337 0.003532495
30 3.035348281 -1.090427337
31 2.014775643 3.035348281
32 4.172935415 2.014775643
33 2.502242278 4.172935415
34 1.721733461 2.502242278
35 -2.723816021 1.721733461
36 1.785208883 -2.723816021
37 1.657276909 1.785208883
38 2.729048982 1.657276909
39 1.449850857 2.729048982
40 1.587504995 1.449850857
41 1.635183638 1.587504995
42 0.899687723 1.635183638
43 1.291242149 0.899687723
44 1.910761908 1.291242149
45 5.252084729 1.910761908
46 1.561856489 5.252084729
47 1.527653124 1.561856489
48 0.692670416 1.527653124
49 0.086320949 0.692670416
50 1.293602270 0.086320949
51 -3.419505251 1.293602270
52 -0.340227522 -3.419505251
53 -2.448641530 -0.340227522
54 1.790118795 -2.448641530
55 1.041951498 1.790118795
56 -0.581817893 1.041951498
57 0.869402831 -0.581817893
58 -0.440399535 0.869402831
59 -1.039078711 -0.440399535
60 -1.481508163 -1.039078711
61 2.046993306 -1.481508163
62 2.017302894 2.046993306
63 -4.473681909 2.017302894
64 -0.896001481 -4.473681909
65 2.215516851 -0.896001481
66 0.784060172 2.215516851
67 -3.363561541 0.784060172
68 -0.515291344 -3.363561541
69 NA -0.515291344
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.386917553 -3.046014439
[2,] 0.121414059 -1.386917553
[3,] 0.142991617 0.121414059
[4,] 0.143200456 0.142991617
[5,] -4.468909368 0.143200456
[6,] 0.948038193 -4.468909368
[7,] 2.414876920 0.948038193
[8,] -4.482214322 2.414876920
[9,] -0.836819604 -4.482214322
[10,] 0.318317068 -0.836819604
[11,] -3.341380957 0.318317068
[12,] -0.621709438 -3.341380957
[13,] -4.723584323 -0.621709438
[14,] -2.402236056 -4.723584323
[15,] -1.550734625 -2.402236056
[16,] -0.866669384 -1.550734625
[17,] -3.534961881 -0.866669384
[18,] 3.361490135 -3.534961881
[19,] 0.926490628 3.361490135
[20,] 0.108697802 0.926490628
[21,] -1.045294314 0.108697802
[22,] -1.114001149 -1.045294314
[23,] -2.199060712 -1.114001149
[24,] -0.532077984 -2.199060712
[25,] -2.198721552 -0.532077984
[26,] 1.097178456 -2.198721552
[27,] 2.006702620 1.097178456
[28,] 0.003532495 2.006702620
[29,] -1.090427337 0.003532495
[30,] 3.035348281 -1.090427337
[31,] 2.014775643 3.035348281
[32,] 4.172935415 2.014775643
[33,] 2.502242278 4.172935415
[34,] 1.721733461 2.502242278
[35,] -2.723816021 1.721733461
[36,] 1.785208883 -2.723816021
[37,] 1.657276909 1.785208883
[38,] 2.729048982 1.657276909
[39,] 1.449850857 2.729048982
[40,] 1.587504995 1.449850857
[41,] 1.635183638 1.587504995
[42,] 0.899687723 1.635183638
[43,] 1.291242149 0.899687723
[44,] 1.910761908 1.291242149
[45,] 5.252084729 1.910761908
[46,] 1.561856489 5.252084729
[47,] 1.527653124 1.561856489
[48,] 0.692670416 1.527653124
[49,] 0.086320949 0.692670416
[50,] 1.293602270 0.086320949
[51,] -3.419505251 1.293602270
[52,] -0.340227522 -3.419505251
[53,] -2.448641530 -0.340227522
[54,] 1.790118795 -2.448641530
[55,] 1.041951498 1.790118795
[56,] -0.581817893 1.041951498
[57,] 0.869402831 -0.581817893
[58,] -0.440399535 0.869402831
[59,] -1.039078711 -0.440399535
[60,] -1.481508163 -1.039078711
[61,] 2.046993306 -1.481508163
[62,] 2.017302894 2.046993306
[63,] -4.473681909 2.017302894
[64,] -0.896001481 -4.473681909
[65,] 2.215516851 -0.896001481
[66,] 0.784060172 2.215516851
[67,] -3.363561541 0.784060172
[68,] -0.515291344 -3.363561541
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.386917553 -3.046014439
2 0.121414059 -1.386917553
3 0.142991617 0.121414059
4 0.143200456 0.142991617
5 -4.468909368 0.143200456
6 0.948038193 -4.468909368
7 2.414876920 0.948038193
8 -4.482214322 2.414876920
9 -0.836819604 -4.482214322
10 0.318317068 -0.836819604
11 -3.341380957 0.318317068
12 -0.621709438 -3.341380957
13 -4.723584323 -0.621709438
14 -2.402236056 -4.723584323
15 -1.550734625 -2.402236056
16 -0.866669384 -1.550734625
17 -3.534961881 -0.866669384
18 3.361490135 -3.534961881
19 0.926490628 3.361490135
20 0.108697802 0.926490628
21 -1.045294314 0.108697802
22 -1.114001149 -1.045294314
23 -2.199060712 -1.114001149
24 -0.532077984 -2.199060712
25 -2.198721552 -0.532077984
26 1.097178456 -2.198721552
27 2.006702620 1.097178456
28 0.003532495 2.006702620
29 -1.090427337 0.003532495
30 3.035348281 -1.090427337
31 2.014775643 3.035348281
32 4.172935415 2.014775643
33 2.502242278 4.172935415
34 1.721733461 2.502242278
35 -2.723816021 1.721733461
36 1.785208883 -2.723816021
37 1.657276909 1.785208883
38 2.729048982 1.657276909
39 1.449850857 2.729048982
40 1.587504995 1.449850857
41 1.635183638 1.587504995
42 0.899687723 1.635183638
43 1.291242149 0.899687723
44 1.910761908 1.291242149
45 5.252084729 1.910761908
46 1.561856489 5.252084729
47 1.527653124 1.561856489
48 0.692670416 1.527653124
49 0.086320949 0.692670416
50 1.293602270 0.086320949
51 -3.419505251 1.293602270
52 -0.340227522 -3.419505251
53 -2.448641530 -0.340227522
54 1.790118795 -2.448641530
55 1.041951498 1.790118795
56 -0.581817893 1.041951498
57 0.869402831 -0.581817893
58 -0.440399535 0.869402831
59 -1.039078711 -0.440399535
60 -1.481508163 -1.039078711
61 2.046993306 -1.481508163
62 2.017302894 2.046993306
63 -4.473681909 2.017302894
64 -0.896001481 -4.473681909
65 2.215516851 -0.896001481
66 0.784060172 2.215516851
67 -3.363561541 0.784060172
68 -0.515291344 -3.363561541
> 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/7apqp1323935115.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/8bl9c1323935115.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/9v0ve1323935115.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/10zfhz1323935115.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/11upc71323935115.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/129kq51323935115.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/13ts3w1323935115.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/149sek1323935115.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/15h59t1323935115.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/16ri6m1323935115.tab")
+ }
>
> try(system("convert tmp/119qi1323935115.ps tmp/119qi1323935115.png",intern=TRUE))
character(0)
> try(system("convert tmp/2jzk21323935115.ps tmp/2jzk21323935115.png",intern=TRUE))
character(0)
> try(system("convert tmp/3zf8d1323935115.ps tmp/3zf8d1323935115.png",intern=TRUE))
character(0)
> try(system("convert tmp/4cj7i1323935115.ps tmp/4cj7i1323935115.png",intern=TRUE))
character(0)
> try(system("convert tmp/5w91y1323935115.ps tmp/5w91y1323935115.png",intern=TRUE))
character(0)
> try(system("convert tmp/63goh1323935115.ps tmp/63goh1323935115.png",intern=TRUE))
character(0)
> try(system("convert tmp/7apqp1323935115.ps tmp/7apqp1323935115.png",intern=TRUE))
character(0)
> try(system("convert tmp/8bl9c1323935115.ps tmp/8bl9c1323935115.png",intern=TRUE))
character(0)
> try(system("convert tmp/9v0ve1323935115.ps tmp/9v0ve1323935115.png",intern=TRUE))
character(0)
> try(system("convert tmp/10zfhz1323935115.ps tmp/10zfhz1323935115.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.371 0.550 3.947