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(1
+ ,13534
+ ,80
+ ,15
+ ,54
+ ,2
+ ,21392
+ ,36
+ ,12
+ ,41
+ ,3
+ ,31809
+ ,39
+ ,13
+ ,48
+ ,4
+ ,41947
+ ,45
+ ,15
+ ,59
+ ,5
+ ,51294
+ ,20
+ ,14
+ ,46
+ ,6
+ ,61020
+ ,23
+ ,11
+ ,43
+ ,7
+ ,71094
+ ,47
+ ,18
+ ,71
+ ,8
+ ,81574
+ ,39
+ ,9
+ ,26
+ ,9
+ ,91312
+ ,21
+ ,17
+ ,62
+ ,10
+ ,101055
+ ,39
+ ,10
+ ,30
+ ,11
+ ,111244
+ ,51
+ ,12
+ ,40
+ ,12
+ ,12890
+ ,29
+ ,10
+ ,34
+ ,13
+ ,131221
+ ,46
+ ,15
+ ,40
+ ,14
+ ,14592
+ ,28
+ ,13
+ ,49
+ ,15
+ ,15960
+ ,30
+ ,8
+ ,32
+ ,16
+ ,161190
+ ,23
+ ,13
+ ,46
+ ,17
+ ,171116
+ ,29
+ ,11
+ ,37
+ ,18
+ ,18916
+ ,40
+ ,12
+ ,48
+ ,19
+ ,191189
+ ,19
+ ,7
+ ,23
+ ,20
+ ,201681
+ ,24
+ ,12
+ ,46
+ ,21
+ ,211061
+ ,42
+ ,9
+ ,32
+ ,22
+ ,22650
+ ,19
+ ,10
+ ,29
+ ,23
+ ,23797
+ ,41
+ ,9
+ ,34
+ ,24
+ ,241405
+ ,37
+ ,12
+ ,28
+ ,25
+ ,25623
+ ,31
+ ,15
+ ,39
+ ,26
+ ,261079
+ ,39
+ ,10
+ ,31
+ ,27
+ ,27744
+ ,26
+ ,8
+ ,23
+ ,28
+ ,28947
+ ,22
+ ,6
+ ,20
+ ,29
+ ,291157
+ ,35
+ ,10
+ ,34
+ ,30
+ ,30854
+ ,15
+ ,10
+ ,35
+ ,31
+ ,31718
+ ,33
+ ,13
+ ,43
+ ,32
+ ,321151
+ ,35
+ ,15
+ ,29
+ ,33
+ ,33558
+ ,10
+ ,11
+ ,37
+ ,34
+ ,34609
+ ,36
+ ,10
+ ,38
+ ,35
+ ,351602
+ ,38
+ ,20
+ ,62
+ ,36
+ ,36880
+ ,11
+ ,13
+ ,45
+ ,37
+ ,371684
+ ,39
+ ,12
+ ,31
+ ,38
+ ,381062
+ ,26
+ ,12
+ ,32
+ ,39
+ ,39688
+ ,18
+ ,9
+ ,35
+ ,40
+ ,40788
+ ,28
+ ,8
+ ,24
+ ,41
+ ,41923
+ ,48
+ ,16
+ ,63
+ ,42
+ ,42538
+ ,26
+ ,7
+ ,21
+ ,43
+ ,43611
+ ,19
+ ,7
+ ,22
+ ,44
+ ,44815
+ ,21
+ ,10
+ ,30
+ ,45
+ ,451140
+ ,14
+ ,16
+ ,55
+ ,46
+ ,46834
+ ,35
+ ,6
+ ,24
+ ,47
+ ,47965
+ ,9
+ ,13
+ ,31
+ ,48
+ ,481183
+ ,35
+ ,8
+ ,24
+ ,49
+ ,49588
+ ,11
+ ,9
+ ,21
+ ,50
+ ,50566
+ ,27
+ ,13
+ ,49
+ ,51
+ ,51626
+ ,8
+ ,10
+ ,30
+ ,52
+ ,521108
+ ,24
+ ,9
+ ,26
+ ,53
+ ,531066
+ ,24
+ ,11
+ ,31
+ ,54
+ ,54929
+ ,31
+ ,12
+ ,41
+ ,55
+ ,56753
+ ,20
+ ,14
+ ,41
+ ,56
+ ,57793
+ ,25
+ ,7
+ ,17
+ ,57
+ ,58853
+ ,11
+ ,12
+ ,25
+ ,58
+ ,59830
+ ,13
+ ,5
+ ,20
+ ,59
+ ,60769
+ ,11
+ ,8
+ ,32
+ ,60
+ ,61771
+ ,24
+ ,12
+ ,44)
+ ,dim=c(5
+ ,60)
+ ,dimnames=list(c('ranking'
+ ,'pages'
+ ,'blogs'
+ ,'pr'
+ ,'lfm')
+ ,1:60))
> y <- array(NA,dim=c(5,60),dimnames=list(c('ranking','pages','blogs','pr','lfm'),1:60))
> 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'
> #'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
> 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
ranking pages blogs pr lfm
1 1 13534 80 15 54
2 2 21392 36 12 41
3 3 31809 39 13 48
4 4 41947 45 15 59
5 5 51294 20 14 46
6 6 61020 23 11 43
7 7 71094 47 18 71
8 8 81574 39 9 26
9 9 91312 21 17 62
10 10 101055 39 10 30
11 11 111244 51 12 40
12 12 12890 29 10 34
13 13 131221 46 15 40
14 14 14592 28 13 49
15 15 15960 30 8 32
16 16 161190 23 13 46
17 17 171116 29 11 37
18 18 18916 40 12 48
19 19 191189 19 7 23
20 20 201681 24 12 46
21 21 211061 42 9 32
22 22 22650 19 10 29
23 23 23797 41 9 34
24 24 241405 37 12 28
25 25 25623 31 15 39
26 26 261079 39 10 31
27 27 27744 26 8 23
28 28 28947 22 6 20
29 29 291157 35 10 34
30 30 30854 15 10 35
31 31 31718 33 13 43
32 32 321151 35 15 29
33 33 33558 10 11 37
34 34 34609 36 10 38
35 35 351602 38 20 62
36 36 36880 11 13 45
37 37 371684 39 12 31
38 38 381062 26 12 32
39 39 39688 18 9 35
40 40 40788 28 8 24
41 41 41923 48 16 63
42 42 42538 26 7 21
43 43 43611 19 7 22
44 44 44815 21 10 30
45 45 451140 14 16 55
46 46 46834 35 6 24
47 47 47965 9 13 31
48 48 481183 35 8 24
49 49 49588 11 9 21
50 50 50566 27 13 49
51 51 51626 8 10 30
52 52 521108 24 9 26
53 53 531066 24 11 31
54 54 54929 31 12 41
55 55 56753 20 14 41
56 56 57793 25 7 17
57 57 58853 11 12 25
58 58 59830 13 5 20
59 59 60769 11 8 32
60 60 61771 24 12 44
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) pages blogs pr lfm
5.600e+01 2.649e-05 -6.135e-01 1.795e-01 -3.558e-01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-26.233 -9.449 -0.282 7.522 32.884
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.600e+01 7.355e+00 7.614 3.73e-10 ***
pages 2.649e-05 1.435e-05 1.846 0.070342 .
blogs -6.135e-01 1.555e-01 -3.945 0.000228 ***
pr 1.795e-01 1.226e+00 0.146 0.884184
lfm -3.558e-01 3.187e-01 -1.117 0.269049
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.34 on 55 degrees of freedom
Multiple R-squared: 0.3719, Adjusted R-squared: 0.3262
F-statistic: 8.14 on 4 and 55 DF, p-value: 3.139e-05
> 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,] 3.763663e-08 7.527326e-08 0.999999962
[2,] 2.807597e-10 5.615195e-10 1.000000000
[3,] 1.242715e-12 2.485431e-12 1.000000000
[4,] 4.109035e-15 8.218070e-15 1.000000000
[5,] 4.051446e-04 8.102892e-04 0.999594855
[6,] 1.165864e-04 2.331727e-04 0.999883414
[7,] 5.806665e-04 1.161333e-03 0.999419333
[8,] 6.808593e-04 1.361719e-03 0.999319141
[9,] 3.762376e-04 7.524752e-04 0.999623762
[10,] 1.852028e-04 3.704056e-04 0.999814797
[11,] 6.975583e-04 1.395117e-03 0.999302442
[12,] 5.674823e-04 1.134965e-03 0.999432518
[13,] 5.064480e-04 1.012896e-03 0.999493552
[14,] 2.818017e-04 5.636034e-04 0.999718198
[15,] 2.108698e-03 4.217396e-03 0.997891302
[16,] 5.281689e-03 1.056338e-02 0.994718311
[17,] 5.372666e-03 1.074533e-02 0.994627334
[18,] 1.722794e-02 3.445588e-02 0.982772060
[19,] 1.693524e-02 3.387048e-02 0.983064761
[20,] 2.837393e-02 5.674787e-02 0.971626066
[21,] 5.011170e-02 1.002234e-01 0.949888299
[22,] 6.645159e-02 1.329032e-01 0.933548409
[23,] 1.646720e-01 3.293441e-01 0.835327965
[24,] 3.282537e-01 6.565074e-01 0.671746298
[25,] 3.050757e-01 6.101513e-01 0.694924350
[26,] 4.951370e-01 9.902740e-01 0.504862979
[27,] 6.851767e-01 6.296466e-01 0.314823300
[28,] 7.318015e-01 5.363970e-01 0.268198481
[29,] 8.700864e-01 2.598272e-01 0.129913605
[30,] 8.542111e-01 2.915778e-01 0.145788889
[31,] 8.708958e-01 2.582085e-01 0.129104239
[32,] 9.424798e-01 1.150405e-01 0.057520235
[33,] 9.629735e-01 7.405292e-02 0.037026460
[34,] 9.861438e-01 2.771241e-02 0.013856204
[35,] 9.904450e-01 1.911002e-02 0.009555008
[36,] 9.946416e-01 1.071684e-02 0.005358420
[37,] 9.967523e-01 6.495410e-03 0.003247705
[38,] 9.970647e-01 5.870525e-03 0.002935263
[39,] 9.981835e-01 3.632946e-03 0.001816473
[40,] 9.978637e-01 4.272504e-03 0.002136252
[41,] 9.959529e-01 8.094134e-03 0.004047067
[42,] 9.960985e-01 7.803065e-03 0.003901533
[43,] 9.965290e-01 6.942074e-03 0.003471037
[44,] 9.988785e-01 2.243010e-03 0.001121505
[45,] 9.928517e-01 1.429662e-02 0.007148311
> postscript(file="/var/wessaorg/rcomp/tmp/1hy431322149621.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/29xa41322149621.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/30n0u1322149621.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/4oad41322149621.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/5vyxr1322149621.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 = 60
Frequency = 1
1 2 3 4 5 6
10.24507907 -20.04521034 -15.16919194 -7.20132714 -26.23326871 -24.17955394
7 8 9 10 11 12
-0.01472066 -18.59834894 -17.52479244 -15.87052966 -4.57879275 -16.24654666
13 14 15 16 17 18
-6.71397494 -10.10585016 -13.06713602 -16.12469859 -14.55028172 0.96534905
19 20 21 22 23 24
-23.48117286 -12.40443591 -5.05318285 -14.41945558 2.00609753 -7.88639506
25 26 27 28 29 30
-1.47493171 -3.75412811 -7.03597895 -9.23050938 -2.93748050 -6.95577362
31 32 33 34 35 36
7.37292303 -3.40863072 -6.56273692 10.89602350 11.47077289 -0.54939243
37 38 39 40 41 42
3.95674117 -2.91150931 3.83018662 7.20132055 32.88376367 7.03986049
43 44 45 46 47 48
4.07269930 8.57620537 2.33640907 17.69464890 3.94807865 7.82871394
49 50 51 52 53 54
4.29149414 24.32759539 7.42011583 4.55460275 6.71109493 27.99875177
55 56 57 58 59 60
21.84288128 18.59882344 12.93103912 14.60917928 18.08904184 30.59044594
> postscript(file="/var/wessaorg/rcomp/tmp/6lh3g1322149621.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 10.24507907 NA
1 -20.04521034 10.24507907
2 -15.16919194 -20.04521034
3 -7.20132714 -15.16919194
4 -26.23326871 -7.20132714
5 -24.17955394 -26.23326871
6 -0.01472066 -24.17955394
7 -18.59834894 -0.01472066
8 -17.52479244 -18.59834894
9 -15.87052966 -17.52479244
10 -4.57879275 -15.87052966
11 -16.24654666 -4.57879275
12 -6.71397494 -16.24654666
13 -10.10585016 -6.71397494
14 -13.06713602 -10.10585016
15 -16.12469859 -13.06713602
16 -14.55028172 -16.12469859
17 0.96534905 -14.55028172
18 -23.48117286 0.96534905
19 -12.40443591 -23.48117286
20 -5.05318285 -12.40443591
21 -14.41945558 -5.05318285
22 2.00609753 -14.41945558
23 -7.88639506 2.00609753
24 -1.47493171 -7.88639506
25 -3.75412811 -1.47493171
26 -7.03597895 -3.75412811
27 -9.23050938 -7.03597895
28 -2.93748050 -9.23050938
29 -6.95577362 -2.93748050
30 7.37292303 -6.95577362
31 -3.40863072 7.37292303
32 -6.56273692 -3.40863072
33 10.89602350 -6.56273692
34 11.47077289 10.89602350
35 -0.54939243 11.47077289
36 3.95674117 -0.54939243
37 -2.91150931 3.95674117
38 3.83018662 -2.91150931
39 7.20132055 3.83018662
40 32.88376367 7.20132055
41 7.03986049 32.88376367
42 4.07269930 7.03986049
43 8.57620537 4.07269930
44 2.33640907 8.57620537
45 17.69464890 2.33640907
46 3.94807865 17.69464890
47 7.82871394 3.94807865
48 4.29149414 7.82871394
49 24.32759539 4.29149414
50 7.42011583 24.32759539
51 4.55460275 7.42011583
52 6.71109493 4.55460275
53 27.99875177 6.71109493
54 21.84288128 27.99875177
55 18.59882344 21.84288128
56 12.93103912 18.59882344
57 14.60917928 12.93103912
58 18.08904184 14.60917928
59 30.59044594 18.08904184
60 NA 30.59044594
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -20.04521034 10.24507907
[2,] -15.16919194 -20.04521034
[3,] -7.20132714 -15.16919194
[4,] -26.23326871 -7.20132714
[5,] -24.17955394 -26.23326871
[6,] -0.01472066 -24.17955394
[7,] -18.59834894 -0.01472066
[8,] -17.52479244 -18.59834894
[9,] -15.87052966 -17.52479244
[10,] -4.57879275 -15.87052966
[11,] -16.24654666 -4.57879275
[12,] -6.71397494 -16.24654666
[13,] -10.10585016 -6.71397494
[14,] -13.06713602 -10.10585016
[15,] -16.12469859 -13.06713602
[16,] -14.55028172 -16.12469859
[17,] 0.96534905 -14.55028172
[18,] -23.48117286 0.96534905
[19,] -12.40443591 -23.48117286
[20,] -5.05318285 -12.40443591
[21,] -14.41945558 -5.05318285
[22,] 2.00609753 -14.41945558
[23,] -7.88639506 2.00609753
[24,] -1.47493171 -7.88639506
[25,] -3.75412811 -1.47493171
[26,] -7.03597895 -3.75412811
[27,] -9.23050938 -7.03597895
[28,] -2.93748050 -9.23050938
[29,] -6.95577362 -2.93748050
[30,] 7.37292303 -6.95577362
[31,] -3.40863072 7.37292303
[32,] -6.56273692 -3.40863072
[33,] 10.89602350 -6.56273692
[34,] 11.47077289 10.89602350
[35,] -0.54939243 11.47077289
[36,] 3.95674117 -0.54939243
[37,] -2.91150931 3.95674117
[38,] 3.83018662 -2.91150931
[39,] 7.20132055 3.83018662
[40,] 32.88376367 7.20132055
[41,] 7.03986049 32.88376367
[42,] 4.07269930 7.03986049
[43,] 8.57620537 4.07269930
[44,] 2.33640907 8.57620537
[45,] 17.69464890 2.33640907
[46,] 3.94807865 17.69464890
[47,] 7.82871394 3.94807865
[48,] 4.29149414 7.82871394
[49,] 24.32759539 4.29149414
[50,] 7.42011583 24.32759539
[51,] 4.55460275 7.42011583
[52,] 6.71109493 4.55460275
[53,] 27.99875177 6.71109493
[54,] 21.84288128 27.99875177
[55,] 18.59882344 21.84288128
[56,] 12.93103912 18.59882344
[57,] 14.60917928 12.93103912
[58,] 18.08904184 14.60917928
[59,] 30.59044594 18.08904184
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -20.04521034 10.24507907
2 -15.16919194 -20.04521034
3 -7.20132714 -15.16919194
4 -26.23326871 -7.20132714
5 -24.17955394 -26.23326871
6 -0.01472066 -24.17955394
7 -18.59834894 -0.01472066
8 -17.52479244 -18.59834894
9 -15.87052966 -17.52479244
10 -4.57879275 -15.87052966
11 -16.24654666 -4.57879275
12 -6.71397494 -16.24654666
13 -10.10585016 -6.71397494
14 -13.06713602 -10.10585016
15 -16.12469859 -13.06713602
16 -14.55028172 -16.12469859
17 0.96534905 -14.55028172
18 -23.48117286 0.96534905
19 -12.40443591 -23.48117286
20 -5.05318285 -12.40443591
21 -14.41945558 -5.05318285
22 2.00609753 -14.41945558
23 -7.88639506 2.00609753
24 -1.47493171 -7.88639506
25 -3.75412811 -1.47493171
26 -7.03597895 -3.75412811
27 -9.23050938 -7.03597895
28 -2.93748050 -9.23050938
29 -6.95577362 -2.93748050
30 7.37292303 -6.95577362
31 -3.40863072 7.37292303
32 -6.56273692 -3.40863072
33 10.89602350 -6.56273692
34 11.47077289 10.89602350
35 -0.54939243 11.47077289
36 3.95674117 -0.54939243
37 -2.91150931 3.95674117
38 3.83018662 -2.91150931
39 7.20132055 3.83018662
40 32.88376367 7.20132055
41 7.03986049 32.88376367
42 4.07269930 7.03986049
43 8.57620537 4.07269930
44 2.33640907 8.57620537
45 17.69464890 2.33640907
46 3.94807865 17.69464890
47 7.82871394 3.94807865
48 4.29149414 7.82871394
49 24.32759539 4.29149414
50 7.42011583 24.32759539
51 4.55460275 7.42011583
52 6.71109493 4.55460275
53 27.99875177 6.71109493
54 21.84288128 27.99875177
55 18.59882344 21.84288128
56 12.93103912 18.59882344
57 14.60917928 12.93103912
58 18.08904184 14.60917928
59 30.59044594 18.08904184
> 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/7jggx1322149621.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/8hr9y1322149621.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/9qd101322149621.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/10bc4c1322149621.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/11046o1322149621.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/12fjta1322149621.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/13iq131322149621.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/1475jr1322149621.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/15l0o21322149621.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/16ry1i1322149621.tab")
+ }
>
> try(system("convert tmp/1hy431322149621.ps tmp/1hy431322149621.png",intern=TRUE))
character(0)
> try(system("convert tmp/29xa41322149621.ps tmp/29xa41322149621.png",intern=TRUE))
character(0)
> try(system("convert tmp/30n0u1322149621.ps tmp/30n0u1322149621.png",intern=TRUE))
character(0)
> try(system("convert tmp/4oad41322149621.ps tmp/4oad41322149621.png",intern=TRUE))
character(0)
> try(system("convert tmp/5vyxr1322149621.ps tmp/5vyxr1322149621.png",intern=TRUE))
character(0)
> try(system("convert tmp/6lh3g1322149621.ps tmp/6lh3g1322149621.png",intern=TRUE))
character(0)
> try(system("convert tmp/7jggx1322149621.ps tmp/7jggx1322149621.png",intern=TRUE))
character(0)
> try(system("convert tmp/8hr9y1322149621.ps tmp/8hr9y1322149621.png",intern=TRUE))
character(0)
> try(system("convert tmp/9qd101322149621.ps tmp/9qd101322149621.png",intern=TRUE))
character(0)
> try(system("convert tmp/10bc4c1322149621.ps tmp/10bc4c1322149621.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.284 0.498 3.927