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(1
+ ,1910
+ ,61
+ ,17
+ ,56
+ ,84
+ ,4
+ ,21
+ ,51
+ ,2
+ ,2598
+ ,74
+ ,19
+ ,73
+ ,47
+ ,3
+ ,15
+ ,48
+ ,3
+ ,2144
+ ,57
+ ,18
+ ,62
+ ,63
+ ,3
+ ,17
+ ,46
+ ,4
+ ,1331
+ ,50
+ ,15
+ ,42
+ ,28
+ ,3
+ ,20
+ ,42
+ ,5
+ ,1431
+ ,48
+ ,15
+ ,59
+ ,22
+ ,2
+ ,12
+ ,38
+ ,6
+ ,7334
+ ,2
+ ,12
+ ,27
+ ,18
+ ,6
+ ,4
+ ,38
+ ,7
+ ,1133
+ ,41
+ ,15
+ ,59
+ ,20
+ ,5
+ ,9
+ ,36
+ ,8
+ ,1195
+ ,31
+ ,20
+ ,78
+ ,27
+ ,5
+ ,11
+ ,35
+ ,9
+ ,1522
+ ,61
+ ,14
+ ,56
+ ,37
+ ,5
+ ,12
+ ,35
+ ,10
+ ,1551
+ ,12
+ ,12
+ ,47
+ ,23
+ ,6
+ ,7
+ ,35
+ ,11
+ ,2108
+ ,46
+ ,13
+ ,51
+ ,67
+ ,5
+ ,14
+ ,34
+ ,12
+ ,1335
+ ,31
+ ,17
+ ,47
+ ,28
+ ,4
+ ,11
+ ,34
+ ,13
+ ,1065
+ ,33
+ ,12
+ ,48
+ ,28
+ ,5
+ ,9
+ ,31
+ ,14
+ ,842
+ ,49
+ ,10
+ ,35
+ ,45
+ ,3
+ ,14
+ ,31
+ ,15
+ ,1539
+ ,15
+ ,13
+ ,47
+ ,15
+ ,5
+ ,4
+ ,31
+ ,16
+ ,1508
+ ,59
+ ,15
+ ,55
+ ,36
+ ,5
+ ,11
+ ,31
+ ,17
+ ,1598
+ ,28
+ ,12
+ ,42
+ ,12
+ ,2
+ ,10
+ ,30
+ ,18
+ ,1219
+ ,55
+ ,16
+ ,55
+ ,30
+ ,6
+ ,9
+ ,30
+ ,19
+ ,1443
+ ,35
+ ,13
+ ,47
+ ,28
+ ,9
+ ,8
+ ,30
+ ,20
+ ,1546
+ ,44
+ ,15
+ ,54
+ ,27
+ ,2
+ ,14
+ ,30
+ ,21
+ ,914
+ ,41
+ ,15
+ ,60
+ ,43
+ ,5
+ ,13
+ ,30
+ ,22
+ ,1370
+ ,26
+ ,13
+ ,51
+ ,10
+ ,3
+ ,10
+ ,28
+ ,23
+ ,1318
+ ,28
+ ,12
+ ,47
+ ,22
+ ,4
+ ,9
+ ,27
+ ,24
+ ,1313
+ ,40
+ ,15
+ ,52
+ ,27
+ ,4
+ ,11
+ ,27
+ ,25
+ ,1743
+ ,28
+ ,12
+ ,38
+ ,21
+ ,11
+ ,7
+ ,27
+ ,26
+ ,1102
+ ,67
+ ,12
+ ,12
+ ,24
+ ,5
+ ,10
+ ,26
+ ,27
+ ,1275
+ ,56
+ ,12
+ ,48
+ ,52
+ ,3
+ ,15
+ ,26
+ ,28
+ ,1253
+ ,54
+ ,12
+ ,48
+ ,24
+ ,5
+ ,7
+ ,26
+ ,29
+ ,1487
+ ,25
+ ,8
+ ,32
+ ,19
+ ,5
+ ,10
+ ,26
+ ,30
+ ,1098
+ ,19
+ ,9
+ ,27
+ ,12
+ ,0
+ ,4
+ ,26
+ ,31
+ ,1176
+ ,36
+ ,12
+ ,47
+ ,21
+ ,3
+ ,10
+ ,25
+ ,32
+ ,903
+ ,42
+ ,16
+ ,58
+ ,71
+ ,4
+ ,13
+ ,25
+ ,33
+ ,1290
+ ,19
+ ,14
+ ,47
+ ,19
+ ,4
+ ,5
+ ,25
+ ,34
+ ,1050
+ ,57
+ ,13
+ ,46
+ ,24
+ ,5
+ ,10
+ ,25
+ ,35
+ ,930
+ ,28
+ ,15
+ ,60
+ ,12
+ ,2
+ ,10
+ ,25
+ ,36
+ ,821
+ ,32
+ ,15
+ ,56
+ ,29
+ ,5
+ ,11
+ ,24
+ ,37
+ ,826
+ ,10
+ ,12
+ ,41
+ ,13
+ ,3
+ ,7
+ ,24
+ ,38
+ ,1402
+ ,28
+ ,12
+ ,45
+ ,22
+ ,11
+ ,6
+ ,24
+ ,39
+ ,1495
+ ,41
+ ,12
+ ,48
+ ,27
+ ,5
+ ,8
+ ,24
+ ,40
+ ,1064
+ ,48
+ ,15
+ ,60
+ ,36
+ ,5
+ ,10
+ ,24
+ ,41
+ ,1469
+ ,57
+ ,12
+ ,48
+ ,27
+ ,3
+ ,9
+ ,24
+ ,42
+ ,1493
+ ,35
+ ,13
+ ,42
+ ,21
+ ,5
+ ,8
+ ,24
+ ,43
+ ,1239
+ ,30
+ ,12
+ ,47
+ ,28
+ ,4
+ ,11
+ ,24
+ ,44
+ ,1317
+ ,39
+ ,12
+ ,41
+ ,17
+ ,3
+ ,5
+ ,23
+ ,45
+ ,708
+ ,17
+ ,15
+ ,49
+ ,15
+ ,8
+ ,5
+ ,23
+ ,46
+ ,872
+ ,33
+ ,12
+ ,39
+ ,26
+ ,3
+ ,10
+ ,23
+ ,47
+ ,853
+ ,55
+ ,12
+ ,39
+ ,19
+ ,3
+ ,8
+ ,23
+ ,48
+ ,1174
+ ,30
+ ,12
+ ,42
+ ,34
+ ,11
+ ,9
+ ,23
+ ,49
+ ,982
+ ,22
+ ,13
+ ,50
+ ,21
+ ,4
+ ,7
+ ,23
+ ,50
+ ,1202
+ ,42
+ ,12
+ ,41
+ ,32
+ ,6
+ ,8
+ ,23
+ ,51
+ ,873
+ ,49
+ ,15
+ ,52
+ ,14
+ ,14
+ ,5
+ ,23
+ ,52
+ ,1000
+ ,13
+ ,9
+ ,36
+ ,17
+ ,6
+ ,5
+ ,22
+ ,53
+ ,1131
+ ,15
+ ,13
+ ,45
+ ,16
+ ,3
+ ,7
+ ,22
+ ,54
+ ,793
+ ,24
+ ,12
+ ,46
+ ,18
+ ,5
+ ,10
+ ,22
+ ,55
+ ,1106
+ ,3
+ ,13
+ ,55
+ ,8
+ ,8
+ ,2
+ ,22
+ ,56
+ ,1205
+ ,35
+ ,13
+ ,49
+ ,30
+ ,8
+ ,5
+ ,22
+ ,57
+ ,1671
+ ,37
+ ,13
+ ,48
+ ,31
+ ,3
+ ,13
+ ,22
+ ,58
+ ,1374
+ ,28
+ ,13
+ ,39
+ ,19
+ ,3
+ ,10
+ ,21
+ ,59
+ ,775
+ ,19
+ ,12
+ ,48
+ ,10
+ ,3
+ ,5
+ ,21
+ ,60
+ ,804
+ ,38
+ ,15
+ ,45
+ ,24
+ ,5
+ ,10
+ ,21
+ ,61
+ ,1224
+ ,29
+ ,14
+ ,52
+ ,28
+ ,6
+ ,8
+ ,21
+ ,62
+ ,1233
+ ,38
+ ,15
+ ,51
+ ,27
+ ,3
+ ,7
+ ,20
+ ,63
+ ,1170
+ ,35
+ ,14
+ ,41
+ ,16
+ ,3
+ ,10
+ ,20
+ ,64
+ ,913
+ ,23
+ ,9
+ ,32
+ ,17
+ ,3
+ ,5
+ ,20
+ ,65
+ ,613
+ ,27
+ ,14
+ ,52
+ ,30
+ ,3
+ ,9
+ ,20
+ ,66
+ ,1204
+ ,32
+ ,16
+ ,54
+ ,20
+ ,4
+ ,6
+ ,19
+ ,67
+ ,933
+ ,7
+ ,9
+ ,27
+ ,10
+ ,5
+ ,6
+ ,18
+ ,68
+ ,861
+ ,57
+ ,12
+ ,41
+ ,30
+ ,3
+ ,9
+ ,18
+ ,69
+ ,932
+ ,39
+ ,12
+ ,45
+ ,34
+ ,5
+ ,11
+ ,18
+ ,70
+ ,705
+ ,18
+ ,13
+ ,52
+ ,13
+ ,13
+ ,6
+ ,18)
+ ,dim=c(9
+ ,70)
+ ,dimnames=list(c('Rang'
+ ,'Pageviews'
+ ,'Blogs'
+ ,'PR'
+ ,'LFM'
+ ,'KCS'
+ ,'SPR'
+ ,'CH'
+ ,'Hours
')
+ ,1:70))
> y <- array(NA,dim=c(9,70),dimnames=list(c('Rang','Pageviews','Blogs','PR','LFM','KCS','SPR','CH','Hours
'),1:70))
> 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'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal Dummies'
> par1 <- '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
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
Rang Pageviews Blogs PR LFM KCS SPR CH Hours\r
1 1 1910 61 17 56 84 4 21 51
2 2 2598 74 19 73 47 3 15 48
3 3 2144 57 18 62 63 3 17 46
4 4 1331 50 15 42 28 3 20 42
5 5 1431 48 15 59 22 2 12 38
6 6 7334 2 12 27 18 6 4 38
7 7 1133 41 15 59 20 5 9 36
8 8 1195 31 20 78 27 5 11 35
9 9 1522 61 14 56 37 5 12 35
10 10 1551 12 12 47 23 6 7 35
11 11 2108 46 13 51 67 5 14 34
12 12 1335 31 17 47 28 4 11 34
13 13 1065 33 12 48 28 5 9 31
14 14 842 49 10 35 45 3 14 31
15 15 1539 15 13 47 15 5 4 31
16 16 1508 59 15 55 36 5 11 31
17 17 1598 28 12 42 12 2 10 30
18 18 1219 55 16 55 30 6 9 30
19 19 1443 35 13 47 28 9 8 30
20 20 1546 44 15 54 27 2 14 30
21 21 914 41 15 60 43 5 13 30
22 22 1370 26 13 51 10 3 10 28
23 23 1318 28 12 47 22 4 9 27
24 24 1313 40 15 52 27 4 11 27
25 25 1743 28 12 38 21 11 7 27
26 26 1102 67 12 12 24 5 10 26
27 27 1275 56 12 48 52 3 15 26
28 28 1253 54 12 48 24 5 7 26
29 29 1487 25 8 32 19 5 10 26
30 30 1098 19 9 27 12 0 4 26
31 31 1176 36 12 47 21 3 10 25
32 32 903 42 16 58 71 4 13 25
33 33 1290 19 14 47 19 4 5 25
34 34 1050 57 13 46 24 5 10 25
35 35 930 28 15 60 12 2 10 25
36 36 821 32 15 56 29 5 11 24
37 37 826 10 12 41 13 3 7 24
38 38 1402 28 12 45 22 11 6 24
39 39 1495 41 12 48 27 5 8 24
40 40 1064 48 15 60 36 5 10 24
41 41 1469 57 12 48 27 3 9 24
42 42 1493 35 13 42 21 5 8 24
43 43 1239 30 12 47 28 4 11 24
44 44 1317 39 12 41 17 3 5 23
45 45 708 17 15 49 15 8 5 23
46 46 872 33 12 39 26 3 10 23
47 47 853 55 12 39 19 3 8 23
48 48 1174 30 12 42 34 11 9 23
49 49 982 22 13 50 21 4 7 23
50 50 1202 42 12 41 32 6 8 23
51 51 873 49 15 52 14 14 5 23
52 52 1000 13 9 36 17 6 5 22
53 53 1131 15 13 45 16 3 7 22
54 54 793 24 12 46 18 5 10 22
55 55 1106 3 13 55 8 8 2 22
56 56 1205 35 13 49 30 8 5 22
57 57 1671 37 13 48 31 3 13 22
58 58 1374 28 13 39 19 3 10 21
59 59 775 19 12 48 10 3 5 21
60 60 804 38 15 45 24 5 10 21
61 61 1224 29 14 52 28 6 8 21
62 62 1233 38 15 51 27 3 7 20
63 63 1170 35 14 41 16 3 10 20
64 64 913 23 9 32 17 3 5 20
65 65 613 27 14 52 30 3 9 20
66 66 1204 32 16 54 20 4 6 19
67 67 933 7 9 27 10 5 6 18
68 68 861 57 12 41 30 3 9 18
69 69 932 39 12 45 34 5 11 18
70 70 705 18 13 52 13 13 6 18
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Pageviews Blogs PR LFM KCS
95.7151231 0.0008113 -0.0967248 1.5407995 -0.1237261 0.1424668
SPR CH `Hours\\r`
0.0742760 0.2847320 -2.9619334
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-17.927 -7.371 -0.782 6.288 23.186
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 95.7151231 7.1939540 13.305 <2e-16 ***
Pageviews 0.0008113 0.0017391 0.466 0.6425
Blogs -0.0967248 0.0937941 -1.031 0.3065
PR 1.5407995 0.8745562 1.762 0.0831 .
LFM -0.1237261 0.1838818 -0.673 0.5036
KCS 0.1424668 0.1205354 1.182 0.2418
SPR 0.0742760 0.4656257 0.160 0.8738
CH 0.2847320 0.5739759 0.496 0.6216
`Hours\\r` -2.9619334 0.2480927 -11.939 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.951 on 61 degrees of freedom
Multiple R-squared: 0.829, Adjusted R-squared: 0.8065
F-statistic: 36.96 on 8 and 61 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.047787150 9.557430e-02 9.522129e-01
[2,] 0.020565766 4.113153e-02 9.794342e-01
[3,] 0.008413962 1.682792e-02 9.915860e-01
[4,] 0.005849434 1.169887e-02 9.941506e-01
[5,] 0.007792755 1.558551e-02 9.922072e-01
[6,] 0.034390272 6.878054e-02 9.656097e-01
[7,] 0.028209874 5.641975e-02 9.717901e-01
[8,] 0.054022955 1.080459e-01 9.459770e-01
[9,] 0.211695398 4.233908e-01 7.883046e-01
[10,] 0.508256115 9.834878e-01 4.917439e-01
[11,] 0.600158135 7.996837e-01 3.998419e-01
[12,] 0.603908657 7.921827e-01 3.960913e-01
[13,] 0.561668444 8.766631e-01 4.383316e-01
[14,] 0.568746849 8.625063e-01 4.312532e-01
[15,] 0.505105850 9.897883e-01 4.948942e-01
[16,] 0.459261607 9.185232e-01 5.407384e-01
[17,] 0.567216154 8.655677e-01 4.327838e-01
[18,] 0.655524598 6.889508e-01 3.444754e-01
[19,] 0.773824853 4.523503e-01 2.261751e-01
[20,] 0.836710137 3.265797e-01 1.632899e-01
[21,] 0.822124327 3.557513e-01 1.778757e-01
[22,] 0.888482376 2.230352e-01 1.115176e-01
[23,] 0.915640540 1.687189e-01 8.435946e-02
[24,] 0.939450555 1.210989e-01 6.054945e-02
[25,] 0.976730768 4.653846e-02 2.326923e-02
[26,] 0.991465181 1.706964e-02 8.534819e-03
[27,] 0.997886945 4.226110e-03 2.113055e-03
[28,] 0.999082813 1.834373e-03 9.171867e-04
[29,] 0.999713383 5.732341e-04 2.866170e-04
[30,] 0.999822874 3.542515e-04 1.771258e-04
[31,] 0.999852035 2.959297e-04 1.479649e-04
[32,] 0.999900940 1.981199e-04 9.905996e-05
[33,] 0.999977355 4.528984e-05 2.264492e-05
[34,] 0.999991736 1.652785e-05 8.263925e-06
[35,] 0.999994993 1.001366e-05 5.006828e-06
[36,] 0.999995355 9.289648e-06 4.644824e-06
[37,] 0.999994870 1.026075e-05 5.130375e-06
[38,] 0.999995571 8.857344e-06 4.428672e-06
[39,] 0.999989979 2.004126e-05 1.002063e-05
[40,] 0.999968100 6.379990e-05 3.189995e-05
[41,] 0.999962990 7.401959e-05 3.700980e-05
[42,] 0.999961919 7.616124e-05 3.808062e-05
[43,] 0.999962966 7.406817e-05 3.703409e-05
[44,] 0.999893777 2.124453e-04 1.062226e-04
[45,] 0.999628491 7.430189e-04 3.715095e-04
[46,] 0.998019173 3.961654e-03 1.980827e-03
[47,] 0.990802471 1.839506e-02 9.197529e-03
> postscript(file="/var/fisher/rcomp/tmp/14i8w1352130891.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/2182r1352130891.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/3gv2n1352130891.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/4t0my1352130891.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/5cn5j1352130891.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 = 70
Frequency = 1
1 2 3 4 5 6
23.1855667 22.0747266 13.2057294 8.6204973 2.8084638 -2.2159682
7 8 9 10 11 12
-0.6344197 -10.5338340 -2.0839474 -0.5350073 -8.8934521 -8.8906440
13 14 15 16 17 18
-8.0410426 -7.5364176 -6.5554242 -8.3511000 -6.4437932 -9.6562890
19 20 21 22 23 24
-6.7930806 -8.2676264 -8.5202931 -7.5927496 -9.7722985 -12.8931105
25 26 27 28 29 30
-9.0386204 -10.7610845 -11.7754469 -4.8326737 -2.7858099 -1.1329684
31 32 33 34 35 36
-6.8751567 -17.9273584 -8.0592443 -3.9821943 -5.1068676 -10.0178746
37 38 39 40 41 42
-4.8164205 -3.6394331 -1.9224212 -4.8850451 1.5100892 -0.9295026
43 44 45 46 47 48
0.1751755 1.6279341 -2.7249912 0.4552826 5.1653749 0.8420745
49 50 51 52 53 54
3.6145118 4.7973771 6.3043345 8.1892646 4.0226065 6.5441801
55 56 57 58 59 60
8.3114649 7.5955225 6.2382560 5.0970068 11.0726240 5.3265450
61 62 63 64 65 66
7.4474945 5.3342846 7.1116976 15.0311497 8.4409542 5.8536311
67 68 69 70
10.4888136 9.9382523 8.3466265 10.6460990
> postscript(file="/var/fisher/rcomp/tmp/6jk3x1352130891.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 = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 23.1855667 NA
1 22.0747266 23.1855667
2 13.2057294 22.0747266
3 8.6204973 13.2057294
4 2.8084638 8.6204973
5 -2.2159682 2.8084638
6 -0.6344197 -2.2159682
7 -10.5338340 -0.6344197
8 -2.0839474 -10.5338340
9 -0.5350073 -2.0839474
10 -8.8934521 -0.5350073
11 -8.8906440 -8.8934521
12 -8.0410426 -8.8906440
13 -7.5364176 -8.0410426
14 -6.5554242 -7.5364176
15 -8.3511000 -6.5554242
16 -6.4437932 -8.3511000
17 -9.6562890 -6.4437932
18 -6.7930806 -9.6562890
19 -8.2676264 -6.7930806
20 -8.5202931 -8.2676264
21 -7.5927496 -8.5202931
22 -9.7722985 -7.5927496
23 -12.8931105 -9.7722985
24 -9.0386204 -12.8931105
25 -10.7610845 -9.0386204
26 -11.7754469 -10.7610845
27 -4.8326737 -11.7754469
28 -2.7858099 -4.8326737
29 -1.1329684 -2.7858099
30 -6.8751567 -1.1329684
31 -17.9273584 -6.8751567
32 -8.0592443 -17.9273584
33 -3.9821943 -8.0592443
34 -5.1068676 -3.9821943
35 -10.0178746 -5.1068676
36 -4.8164205 -10.0178746
37 -3.6394331 -4.8164205
38 -1.9224212 -3.6394331
39 -4.8850451 -1.9224212
40 1.5100892 -4.8850451
41 -0.9295026 1.5100892
42 0.1751755 -0.9295026
43 1.6279341 0.1751755
44 -2.7249912 1.6279341
45 0.4552826 -2.7249912
46 5.1653749 0.4552826
47 0.8420745 5.1653749
48 3.6145118 0.8420745
49 4.7973771 3.6145118
50 6.3043345 4.7973771
51 8.1892646 6.3043345
52 4.0226065 8.1892646
53 6.5441801 4.0226065
54 8.3114649 6.5441801
55 7.5955225 8.3114649
56 6.2382560 7.5955225
57 5.0970068 6.2382560
58 11.0726240 5.0970068
59 5.3265450 11.0726240
60 7.4474945 5.3265450
61 5.3342846 7.4474945
62 7.1116976 5.3342846
63 15.0311497 7.1116976
64 8.4409542 15.0311497
65 5.8536311 8.4409542
66 10.4888136 5.8536311
67 9.9382523 10.4888136
68 8.3466265 9.9382523
69 10.6460990 8.3466265
70 NA 10.6460990
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 22.0747266 23.1855667
[2,] 13.2057294 22.0747266
[3,] 8.6204973 13.2057294
[4,] 2.8084638 8.6204973
[5,] -2.2159682 2.8084638
[6,] -0.6344197 -2.2159682
[7,] -10.5338340 -0.6344197
[8,] -2.0839474 -10.5338340
[9,] -0.5350073 -2.0839474
[10,] -8.8934521 -0.5350073
[11,] -8.8906440 -8.8934521
[12,] -8.0410426 -8.8906440
[13,] -7.5364176 -8.0410426
[14,] -6.5554242 -7.5364176
[15,] -8.3511000 -6.5554242
[16,] -6.4437932 -8.3511000
[17,] -9.6562890 -6.4437932
[18,] -6.7930806 -9.6562890
[19,] -8.2676264 -6.7930806
[20,] -8.5202931 -8.2676264
[21,] -7.5927496 -8.5202931
[22,] -9.7722985 -7.5927496
[23,] -12.8931105 -9.7722985
[24,] -9.0386204 -12.8931105
[25,] -10.7610845 -9.0386204
[26,] -11.7754469 -10.7610845
[27,] -4.8326737 -11.7754469
[28,] -2.7858099 -4.8326737
[29,] -1.1329684 -2.7858099
[30,] -6.8751567 -1.1329684
[31,] -17.9273584 -6.8751567
[32,] -8.0592443 -17.9273584
[33,] -3.9821943 -8.0592443
[34,] -5.1068676 -3.9821943
[35,] -10.0178746 -5.1068676
[36,] -4.8164205 -10.0178746
[37,] -3.6394331 -4.8164205
[38,] -1.9224212 -3.6394331
[39,] -4.8850451 -1.9224212
[40,] 1.5100892 -4.8850451
[41,] -0.9295026 1.5100892
[42,] 0.1751755 -0.9295026
[43,] 1.6279341 0.1751755
[44,] -2.7249912 1.6279341
[45,] 0.4552826 -2.7249912
[46,] 5.1653749 0.4552826
[47,] 0.8420745 5.1653749
[48,] 3.6145118 0.8420745
[49,] 4.7973771 3.6145118
[50,] 6.3043345 4.7973771
[51,] 8.1892646 6.3043345
[52,] 4.0226065 8.1892646
[53,] 6.5441801 4.0226065
[54,] 8.3114649 6.5441801
[55,] 7.5955225 8.3114649
[56,] 6.2382560 7.5955225
[57,] 5.0970068 6.2382560
[58,] 11.0726240 5.0970068
[59,] 5.3265450 11.0726240
[60,] 7.4474945 5.3265450
[61,] 5.3342846 7.4474945
[62,] 7.1116976 5.3342846
[63,] 15.0311497 7.1116976
[64,] 8.4409542 15.0311497
[65,] 5.8536311 8.4409542
[66,] 10.4888136 5.8536311
[67,] 9.9382523 10.4888136
[68,] 8.3466265 9.9382523
[69,] 10.6460990 8.3466265
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 22.0747266 23.1855667
2 13.2057294 22.0747266
3 8.6204973 13.2057294
4 2.8084638 8.6204973
5 -2.2159682 2.8084638
6 -0.6344197 -2.2159682
7 -10.5338340 -0.6344197
8 -2.0839474 -10.5338340
9 -0.5350073 -2.0839474
10 -8.8934521 -0.5350073
11 -8.8906440 -8.8934521
12 -8.0410426 -8.8906440
13 -7.5364176 -8.0410426
14 -6.5554242 -7.5364176
15 -8.3511000 -6.5554242
16 -6.4437932 -8.3511000
17 -9.6562890 -6.4437932
18 -6.7930806 -9.6562890
19 -8.2676264 -6.7930806
20 -8.5202931 -8.2676264
21 -7.5927496 -8.5202931
22 -9.7722985 -7.5927496
23 -12.8931105 -9.7722985
24 -9.0386204 -12.8931105
25 -10.7610845 -9.0386204
26 -11.7754469 -10.7610845
27 -4.8326737 -11.7754469
28 -2.7858099 -4.8326737
29 -1.1329684 -2.7858099
30 -6.8751567 -1.1329684
31 -17.9273584 -6.8751567
32 -8.0592443 -17.9273584
33 -3.9821943 -8.0592443
34 -5.1068676 -3.9821943
35 -10.0178746 -5.1068676
36 -4.8164205 -10.0178746
37 -3.6394331 -4.8164205
38 -1.9224212 -3.6394331
39 -4.8850451 -1.9224212
40 1.5100892 -4.8850451
41 -0.9295026 1.5100892
42 0.1751755 -0.9295026
43 1.6279341 0.1751755
44 -2.7249912 1.6279341
45 0.4552826 -2.7249912
46 5.1653749 0.4552826
47 0.8420745 5.1653749
48 3.6145118 0.8420745
49 4.7973771 3.6145118
50 6.3043345 4.7973771
51 8.1892646 6.3043345
52 4.0226065 8.1892646
53 6.5441801 4.0226065
54 8.3114649 6.5441801
55 7.5955225 8.3114649
56 6.2382560 7.5955225
57 5.0970068 6.2382560
58 11.0726240 5.0970068
59 5.3265450 11.0726240
60 7.4474945 5.3265450
61 5.3342846 7.4474945
62 7.1116976 5.3342846
63 15.0311497 7.1116976
64 8.4409542 15.0311497
65 5.8536311 8.4409542
66 10.4888136 5.8536311
67 9.9382523 10.4888136
68 8.3466265 9.9382523
69 10.6460990 8.3466265
> 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/7z89q1352130891.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/8tulj1352130891.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/9mzjp1352130891.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/10k8ew1352130891.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/11ud501352130891.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/127hhq1352130891.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/13ujab1352130891.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/14ufwr1352130892.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/15bed21352130892.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/16jt2h1352130892.tab")
+ }
>
> try(system("convert tmp/14i8w1352130891.ps tmp/14i8w1352130891.png",intern=TRUE))
character(0)
> try(system("convert tmp/2182r1352130891.ps tmp/2182r1352130891.png",intern=TRUE))
character(0)
> try(system("convert tmp/3gv2n1352130891.ps tmp/3gv2n1352130891.png",intern=TRUE))
character(0)
> try(system("convert tmp/4t0my1352130891.ps tmp/4t0my1352130891.png",intern=TRUE))
character(0)
> try(system("convert tmp/5cn5j1352130891.ps tmp/5cn5j1352130891.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jk3x1352130891.ps tmp/6jk3x1352130891.png",intern=TRUE))
character(0)
> try(system("convert tmp/7z89q1352130891.ps tmp/7z89q1352130891.png",intern=TRUE))
character(0)
> try(system("convert tmp/8tulj1352130891.ps tmp/8tulj1352130891.png",intern=TRUE))
character(0)
> try(system("convert tmp/9mzjp1352130891.ps tmp/9mzjp1352130891.png",intern=TRUE))
character(0)
> try(system("convert tmp/10k8ew1352130891.ps tmp/10k8ew1352130891.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.126 1.112 7.225