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(869
+ ,58
+ ,28
+ ,103
+ ,84786
+ ,98364
+ ,120982
+ ,2172
+ ,108
+ ,30
+ ,103
+ ,101193
+ ,96933
+ ,179321
+ ,901
+ ,49
+ ,22
+ ,51
+ ,38361
+ ,79234
+ ,123185
+ ,463
+ ,0
+ ,26
+ ,70
+ ,68504
+ ,42551
+ ,52746
+ ,371
+ ,1
+ ,18
+ ,22
+ ,22807
+ ,6853
+ ,33170
+ ,1495
+ ,86
+ ,44
+ ,148
+ ,71701
+ ,75851
+ ,173326
+ ,2187
+ ,104
+ ,40
+ ,124
+ ,80444
+ ,93163
+ ,258873
+ ,1491
+ ,63
+ ,34
+ ,70
+ ,53855
+ ,96037
+ ,180083
+ ,1036
+ ,82
+ ,23
+ ,66
+ ,99645
+ ,94728
+ ,135473
+ ,1882
+ ,115
+ ,36
+ ,134
+ ,114789
+ ,105499
+ ,202925
+ ,1220
+ ,50
+ ,25
+ ,84
+ ,65553
+ ,98958
+ ,153935
+ ,1289
+ ,83
+ ,39
+ ,156
+ ,97500
+ ,77900
+ ,132943
+ ,1812
+ ,105
+ ,33
+ ,110
+ ,77873
+ ,178812
+ ,221698
+ ,1731
+ ,114
+ ,43
+ ,158
+ ,90183
+ ,163253
+ ,260561
+ ,807
+ ,38
+ ,30
+ ,109
+ ,61542
+ ,27032
+ ,84853
+ ,1940
+ ,71
+ ,32
+ ,92
+ ,55813
+ ,86572
+ ,215641
+ ,1499
+ ,59
+ ,28
+ ,70
+ ,55461
+ ,85371
+ ,167542
+ ,2747
+ ,106
+ ,30
+ ,93
+ ,70106
+ ,120642
+ ,269651
+ ,2099
+ ,34
+ ,39
+ ,31
+ ,71570
+ ,78348
+ ,116408
+ ,918
+ ,20
+ ,26
+ ,66
+ ,33032
+ ,56968
+ ,78800
+ ,3373
+ ,115
+ ,39
+ ,133
+ ,139077
+ ,161632
+ ,277965
+ ,1713
+ ,85
+ ,33
+ ,113
+ ,71595
+ ,87850
+ ,150629
+ ,1438
+ ,76
+ ,28
+ ,100
+ ,72260
+ ,127969
+ ,168809
+ ,496
+ ,8
+ ,4
+ ,7
+ ,5950
+ ,15049
+ ,24188
+ ,744
+ ,21
+ ,18
+ ,61
+ ,32551
+ ,25109
+ ,65029
+ ,1161
+ ,30
+ ,14
+ ,41
+ ,31701
+ ,45824
+ ,101097
+ ,2694
+ ,92
+ ,28
+ ,102
+ ,120733
+ ,162647
+ ,233328
+ ,1769
+ ,75
+ ,28
+ ,99
+ ,73107
+ ,60622
+ ,206161
+ ,3148
+ ,128
+ ,38
+ ,129
+ ,132068
+ ,179566
+ ,311473
+ ,2474
+ ,105
+ ,23
+ ,62
+ ,149193
+ ,184301
+ ,235800
+ ,2084
+ ,55
+ ,36
+ ,73
+ ,46821
+ ,75661
+ ,177939
+ ,1954
+ ,56
+ ,32
+ ,114
+ ,87011
+ ,96144
+ ,207176
+ ,1389
+ ,72
+ ,25
+ ,70
+ ,55183
+ ,117286
+ ,174184
+ ,2269
+ ,75
+ ,36
+ ,116
+ ,73511
+ ,109377
+ ,187559
+ ,1268
+ ,118
+ ,23
+ ,74
+ ,78664
+ ,73631
+ ,119016
+ ,1943
+ ,77
+ ,40
+ ,138
+ ,70054
+ ,86767
+ ,182192
+ ,1762
+ ,66
+ ,40
+ ,151
+ ,74011
+ ,93487
+ ,194979
+ ,1857
+ ,116
+ ,33
+ ,115
+ ,93133
+ ,94552
+ ,275541
+ ,1502
+ ,73
+ ,34
+ ,104
+ ,225920
+ ,128754
+ ,182999
+ ,1441
+ ,99
+ ,30
+ ,108
+ ,62133
+ ,66363
+ ,135649
+ ,1416
+ ,53
+ ,22
+ ,69
+ ,43836
+ ,61724
+ ,120221
+ ,1317
+ ,30
+ ,26
+ ,99
+ ,38692
+ ,68580
+ ,145790
+ ,870
+ ,49
+ ,8
+ ,27
+ ,56622
+ ,55792
+ ,80953
+ ,2008
+ ,75
+ ,45
+ ,93
+ ,67267
+ ,129484
+ ,241066
+ ,1885
+ ,68
+ ,33
+ ,69
+ ,41140
+ ,87831
+ ,204713
+ ,1369
+ ,81
+ ,28
+ ,99
+ ,138599
+ ,136048
+ ,182613
+ ,2845
+ ,130
+ ,24
+ ,85
+ ,162901
+ ,186646
+ ,310839
+ ,1391
+ ,39
+ ,32
+ ,50
+ ,37510
+ ,64219
+ ,144966
+ ,602
+ ,13
+ ,19
+ ,64
+ ,43750
+ ,19630
+ ,43287
+ ,1743
+ ,74
+ ,20
+ ,31
+ ,40652
+ ,76825
+ ,155754
+ ,2014
+ ,109
+ ,31
+ ,92
+ ,85872
+ ,109427
+ ,201940
+ ,2143
+ ,151
+ ,32
+ ,106
+ ,89275
+ ,118168
+ ,235454
+ ,1401
+ ,54
+ ,31
+ ,114
+ ,120662
+ ,146304
+ ,224549
+ ,530
+ ,23
+ ,11
+ ,30
+ ,25162
+ ,24610
+ ,61857
+ ,387
+ ,4
+ ,0
+ ,0
+ ,855
+ ,6622
+ ,21054
+ ,1742
+ ,62
+ ,24
+ ,60
+ ,97068
+ ,115814
+ ,209641
+ ,449
+ ,18
+ ,8
+ ,9
+ ,14116
+ ,13155
+ ,31414
+ ,1606
+ ,64
+ ,40
+ ,140
+ ,110681
+ ,68847
+ ,184510
+ ,568
+ ,16
+ ,8
+ ,21
+ ,8773
+ ,13983
+ ,38214
+ ,1459
+ ,48
+ ,35
+ ,124
+ ,83209
+ ,65176
+ ,151101
+ ,1955
+ ,130
+ ,38
+ ,120
+ ,103487
+ ,180759
+ ,250579
+ ,1002
+ ,59
+ ,31
+ ,114
+ ,71220
+ ,100226
+ ,158015
+ ,956
+ ,32
+ ,28
+ ,78
+ ,56926
+ ,54454
+ ,85439
+ ,3604
+ ,95
+ ,40
+ ,141
+ ,115168
+ ,170745
+ ,351619
+ ,1035
+ ,14
+ ,30
+ ,101
+ ,111194
+ ,6940
+ ,84207
+ ,1701
+ ,70
+ ,32
+ ,90
+ ,51633
+ ,86839
+ ,165543
+ ,1249
+ ,19
+ ,27
+ ,36
+ ,75345
+ ,44830
+ ,141722
+ ,3352
+ ,91
+ ,31
+ ,97
+ ,98952
+ ,103300
+ ,299775
+ ,1369
+ ,135
+ ,41
+ ,148
+ ,123969
+ ,58106
+ ,104389
+ ,2201
+ ,87
+ ,32
+ ,105
+ ,135400
+ ,122422
+ ,199476
+ ,207
+ ,4
+ ,0
+ ,0
+ ,6023
+ ,7953
+ ,14688
+ ,151
+ ,7
+ ,0
+ ,0
+ ,1644
+ ,4245
+ ,7199
+ ,474
+ ,12
+ ,5
+ ,13
+ ,6179
+ ,21509
+ ,46660
+ ,141
+ ,0
+ ,1
+ ,4
+ ,3926
+ ,7670
+ ,17547
+ ,1318
+ ,46
+ ,24
+ ,46
+ ,73224
+ ,53608
+ ,152601)
+ ,dim=c(7
+ ,75)
+ ,dimnames=list(c('pageviews'
+ ,'blogged_computations'
+ ,'compendiums_reviewed'
+ ,'feedback_messages'
+ ,'totsize'
+ ,'totseconds'
+ ,'time_in_rfc')
+ ,1:75))
> y <- array(NA,dim=c(7,75),dimnames=list(c('pageviews','blogged_computations','compendiums_reviewed','feedback_messages','totsize','totseconds','time_in_rfc'),1:75))
> 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 = '4'
> #'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
feedback_messages pageviews blogged_computations compendiums_reviewed
1 103 869 58 28
2 103 2172 108 30
3 51 901 49 22
4 70 463 0 26
5 22 371 1 18
6 148 1495 86 44
7 124 2187 104 40
8 70 1491 63 34
9 66 1036 82 23
10 134 1882 115 36
11 84 1220 50 25
12 156 1289 83 39
13 110 1812 105 33
14 158 1731 114 43
15 109 807 38 30
16 92 1940 71 32
17 70 1499 59 28
18 93 2747 106 30
19 31 2099 34 39
20 66 918 20 26
21 133 3373 115 39
22 113 1713 85 33
23 100 1438 76 28
24 7 496 8 4
25 61 744 21 18
26 41 1161 30 14
27 102 2694 92 28
28 99 1769 75 28
29 129 3148 128 38
30 62 2474 105 23
31 73 2084 55 36
32 114 1954 56 32
33 70 1389 72 25
34 116 2269 75 36
35 74 1268 118 23
36 138 1943 77 40
37 151 1762 66 40
38 115 1857 116 33
39 104 1502 73 34
40 108 1441 99 30
41 69 1416 53 22
42 99 1317 30 26
43 27 870 49 8
44 93 2008 75 45
45 69 1885 68 33
46 99 1369 81 28
47 85 2845 130 24
48 50 1391 39 32
49 64 602 13 19
50 31 1743 74 20
51 92 2014 109 31
52 106 2143 151 32
53 114 1401 54 31
54 30 530 23 11
55 0 387 4 0
56 60 1742 62 24
57 9 449 18 8
58 140 1606 64 40
59 21 568 16 8
60 124 1459 48 35
61 120 1955 130 38
62 114 1002 59 31
63 78 956 32 28
64 141 3604 95 40
65 101 1035 14 30
66 90 1701 70 32
67 36 1249 19 27
68 97 3352 91 31
69 148 1369 135 41
70 105 2201 87 32
71 0 207 4 0
72 0 151 7 0
73 13 474 12 5
74 4 141 0 1
75 46 1318 46 24
totsize totseconds time_in_rfc
1 84786 98364 120982
2 101193 96933 179321
3 38361 79234 123185
4 68504 42551 52746
5 22807 6853 33170
6 71701 75851 173326
7 80444 93163 258873
8 53855 96037 180083
9 99645 94728 135473
10 114789 105499 202925
11 65553 98958 153935
12 97500 77900 132943
13 77873 178812 221698
14 90183 163253 260561
15 61542 27032 84853
16 55813 86572 215641
17 55461 85371 167542
18 70106 120642 269651
19 71570 78348 116408
20 33032 56968 78800
21 139077 161632 277965
22 71595 87850 150629
23 72260 127969 168809
24 5950 15049 24188
25 32551 25109 65029
26 31701 45824 101097
27 120733 162647 233328
28 73107 60622 206161
29 132068 179566 311473
30 149193 184301 235800
31 46821 75661 177939
32 87011 96144 207176
33 55183 117286 174184
34 73511 109377 187559
35 78664 73631 119016
36 70054 86767 182192
37 74011 93487 194979
38 93133 94552 275541
39 225920 128754 182999
40 62133 66363 135649
41 43836 61724 120221
42 38692 68580 145790
43 56622 55792 80953
44 67267 129484 241066
45 41140 87831 204713
46 138599 136048 182613
47 162901 186646 310839
48 37510 64219 144966
49 43750 19630 43287
50 40652 76825 155754
51 85872 109427 201940
52 89275 118168 235454
53 120662 146304 224549
54 25162 24610 61857
55 855 6622 21054
56 97068 115814 209641
57 14116 13155 31414
58 110681 68847 184510
59 8773 13983 38214
60 83209 65176 151101
61 103487 180759 250579
62 71220 100226 158015
63 56926 54454 85439
64 115168 170745 351619
65 111194 6940 84207
66 51633 86839 165543
67 75345 44830 141722
68 98952 103300 299775
69 123969 58106 104389
70 135400 122422 199476
71 6023 7953 14688
72 1644 4245 7199
73 6179 21509 46660
74 3926 7670 17547
75 73224 53608 152601
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) pageviews blogged_computations
-6.1210762 -0.0162297 0.2960221
compendiums_reviewed totsize totseconds
2.7089780 0.0001728 -0.0001127
time_in_rfc
0.0001140
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-61.337 -10.310 5.635 12.601 33.340
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.121e+00 5.886e+00 -1.040 0.30203
pageviews -1.623e-02 7.341e-03 -2.211 0.03041 *
blogged_computations 2.960e-01 1.036e-01 2.857 0.00567 **
compendiums_reviewed 2.709e+00 2.971e-01 9.118 2.07e-13 ***
totsize 1.728e-04 8.140e-05 2.123 0.03739 *
totseconds -1.127e-04 1.059e-04 -1.064 0.29104
time_in_rfc 1.140e-04 9.499e-05 1.200 0.23438
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18.81 on 68 degrees of freedom
Multiple R-squared: 0.8167, Adjusted R-squared: 0.8005
F-statistic: 50.48 on 6 and 68 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.45075860 0.90151721 0.54924140
[2,] 0.48980642 0.97961284 0.51019358
[3,] 0.42544357 0.85088714 0.57455643
[4,] 0.29763263 0.59526527 0.70236737
[5,] 0.20808351 0.41616702 0.79191649
[6,] 0.22905697 0.45811393 0.77094303
[7,] 0.18230193 0.36460386 0.81769807
[8,] 0.12181748 0.24363495 0.87818252
[9,] 0.09978827 0.19957654 0.90021173
[10,] 0.62971419 0.74057162 0.37028581
[11,] 0.64219158 0.71561685 0.35780842
[12,] 0.69585378 0.60829244 0.30414622
[13,] 0.67043860 0.65912279 0.32956140
[14,] 0.63382053 0.73235895 0.36617947
[15,] 0.61533122 0.76933756 0.38466878
[16,] 0.60149690 0.79700621 0.39850310
[17,] 0.54434821 0.91130357 0.45565179
[18,] 0.52492514 0.95014972 0.47507486
[19,] 0.45606413 0.91212826 0.54393587
[20,] 0.38846773 0.77693546 0.61153227
[21,] 0.46201390 0.92402781 0.53798610
[22,] 0.46004837 0.92009674 0.53995163
[23,] 0.49409720 0.98819440 0.50590280
[24,] 0.42990967 0.85981933 0.57009033
[25,] 0.44000146 0.88000293 0.55999854
[26,] 0.44314212 0.88628423 0.55685788
[27,] 0.45860317 0.91720634 0.54139683
[28,] 0.59499608 0.81000783 0.40500392
[29,] 0.60067448 0.79865105 0.39932552
[30,] 0.69631331 0.60737338 0.30368669
[31,] 0.66499092 0.67001815 0.33500908
[32,] 0.60799695 0.78400609 0.39200305
[33,] 0.76083601 0.47832798 0.23916399
[34,] 0.70432256 0.59135487 0.29567744
[35,] 0.81000517 0.37998965 0.18999483
[36,] 0.80764743 0.38470514 0.19235257
[37,] 0.75647450 0.48705100 0.24352550
[38,] 0.72770300 0.54459400 0.27229700
[39,] 0.84898796 0.30202409 0.15101204
[40,] 0.81222211 0.37555578 0.18777789
[41,] 0.87255029 0.25489942 0.12744971
[42,] 0.84111734 0.31776532 0.15888266
[43,] 0.79364071 0.41271858 0.20635929
[44,] 0.78875557 0.42248886 0.21124443
[45,] 0.71946406 0.56107188 0.28053594
[46,] 0.65076956 0.69846089 0.34923044
[47,] 0.59863478 0.80273044 0.40136522
[48,] 0.53814369 0.92371263 0.46185631
[49,] 0.55081090 0.89837819 0.44918910
[50,] 0.45086978 0.90173956 0.54913022
[51,] 0.47979552 0.95959105 0.52020448
[52,] 0.41501357 0.83002713 0.58498643
[53,] 0.77083172 0.45833655 0.22916828
[54,] 0.69239084 0.61521832 0.30760916
[55,] 0.90185090 0.19629820 0.09814910
[56,] 0.98089688 0.03820625 0.01910312
> postscript(file="/var/wessaorg/rcomp/tmp/1egwf1324047132.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/2gkun1324047132.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/30rpq1324047132.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/4ozgh1324047132.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/5k9fo1324047132.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 = 75
Frequency = 1
1 2 3 4 5 6
12.8441172 4.1267577 -14.1016066 0.1453587 -21.8651680 10.1310756
7 8 9 10 11 12
-6.4410527 -29.4475843 -19.6341987 8.0193363 9.6708367 29.5955595
13 14 15 16 17 18
6.4693915 15.0914331 18.4392369 -2.5678956 -11.9295730 1.7990599
19 20 21 22 23 24
-61.3372354 2.3942481 16.6646087 12.7207562 13.7992719 5.8772208
25 26 27 28 29 30
14.0097157 7.3189793 19.6237337 6.4770935 7.2862587 -17.0108922
31 32 33 34 35 36
-20.7086367 20.7512135 -6.5498416 17.4627388 -15.4008674 21.4059495
37 38 39 40 41 42
33.3404674 -9.3223167 -24.6115394 8.2108208 8.4916819 31.6047262
43 44 45 46 47 48
-1.6623990 -36.9080883 -24.3584065 -1.9282569 -8.7562785 -35.3054417
49 50 51 52 53 54
14.2896723 -26.7980488 -10.9652484 -13.4365501 12.9324280 -0.5103915
55 56 57 58 59 60
9.4165483 -16.5968601 -9.1298389 12.4809332 5.6351639 20.5181250
61 62 63 64 65 66
-9.6539777 15.9131973 0.8716295 28.3891603 10.4736890 -1.6885947
67 68 69 70 71 72
-40.4977069 6.9774321 -1.4669486 2.0587226 6.4775761 5.8732234
73 74 75
5.7541048 7.8861999 -29.1280082
> postscript(file="/var/wessaorg/rcomp/tmp/6n2wk1324047132.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 = 75
Frequency = 1
lag(myerror, k = 1) myerror
0 12.8441172 NA
1 4.1267577 12.8441172
2 -14.1016066 4.1267577
3 0.1453587 -14.1016066
4 -21.8651680 0.1453587
5 10.1310756 -21.8651680
6 -6.4410527 10.1310756
7 -29.4475843 -6.4410527
8 -19.6341987 -29.4475843
9 8.0193363 -19.6341987
10 9.6708367 8.0193363
11 29.5955595 9.6708367
12 6.4693915 29.5955595
13 15.0914331 6.4693915
14 18.4392369 15.0914331
15 -2.5678956 18.4392369
16 -11.9295730 -2.5678956
17 1.7990599 -11.9295730
18 -61.3372354 1.7990599
19 2.3942481 -61.3372354
20 16.6646087 2.3942481
21 12.7207562 16.6646087
22 13.7992719 12.7207562
23 5.8772208 13.7992719
24 14.0097157 5.8772208
25 7.3189793 14.0097157
26 19.6237337 7.3189793
27 6.4770935 19.6237337
28 7.2862587 6.4770935
29 -17.0108922 7.2862587
30 -20.7086367 -17.0108922
31 20.7512135 -20.7086367
32 -6.5498416 20.7512135
33 17.4627388 -6.5498416
34 -15.4008674 17.4627388
35 21.4059495 -15.4008674
36 33.3404674 21.4059495
37 -9.3223167 33.3404674
38 -24.6115394 -9.3223167
39 8.2108208 -24.6115394
40 8.4916819 8.2108208
41 31.6047262 8.4916819
42 -1.6623990 31.6047262
43 -36.9080883 -1.6623990
44 -24.3584065 -36.9080883
45 -1.9282569 -24.3584065
46 -8.7562785 -1.9282569
47 -35.3054417 -8.7562785
48 14.2896723 -35.3054417
49 -26.7980488 14.2896723
50 -10.9652484 -26.7980488
51 -13.4365501 -10.9652484
52 12.9324280 -13.4365501
53 -0.5103915 12.9324280
54 9.4165483 -0.5103915
55 -16.5968601 9.4165483
56 -9.1298389 -16.5968601
57 12.4809332 -9.1298389
58 5.6351639 12.4809332
59 20.5181250 5.6351639
60 -9.6539777 20.5181250
61 15.9131973 -9.6539777
62 0.8716295 15.9131973
63 28.3891603 0.8716295
64 10.4736890 28.3891603
65 -1.6885947 10.4736890
66 -40.4977069 -1.6885947
67 6.9774321 -40.4977069
68 -1.4669486 6.9774321
69 2.0587226 -1.4669486
70 6.4775761 2.0587226
71 5.8732234 6.4775761
72 5.7541048 5.8732234
73 7.8861999 5.7541048
74 -29.1280082 7.8861999
75 NA -29.1280082
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4.1267577 12.8441172
[2,] -14.1016066 4.1267577
[3,] 0.1453587 -14.1016066
[4,] -21.8651680 0.1453587
[5,] 10.1310756 -21.8651680
[6,] -6.4410527 10.1310756
[7,] -29.4475843 -6.4410527
[8,] -19.6341987 -29.4475843
[9,] 8.0193363 -19.6341987
[10,] 9.6708367 8.0193363
[11,] 29.5955595 9.6708367
[12,] 6.4693915 29.5955595
[13,] 15.0914331 6.4693915
[14,] 18.4392369 15.0914331
[15,] -2.5678956 18.4392369
[16,] -11.9295730 -2.5678956
[17,] 1.7990599 -11.9295730
[18,] -61.3372354 1.7990599
[19,] 2.3942481 -61.3372354
[20,] 16.6646087 2.3942481
[21,] 12.7207562 16.6646087
[22,] 13.7992719 12.7207562
[23,] 5.8772208 13.7992719
[24,] 14.0097157 5.8772208
[25,] 7.3189793 14.0097157
[26,] 19.6237337 7.3189793
[27,] 6.4770935 19.6237337
[28,] 7.2862587 6.4770935
[29,] -17.0108922 7.2862587
[30,] -20.7086367 -17.0108922
[31,] 20.7512135 -20.7086367
[32,] -6.5498416 20.7512135
[33,] 17.4627388 -6.5498416
[34,] -15.4008674 17.4627388
[35,] 21.4059495 -15.4008674
[36,] 33.3404674 21.4059495
[37,] -9.3223167 33.3404674
[38,] -24.6115394 -9.3223167
[39,] 8.2108208 -24.6115394
[40,] 8.4916819 8.2108208
[41,] 31.6047262 8.4916819
[42,] -1.6623990 31.6047262
[43,] -36.9080883 -1.6623990
[44,] -24.3584065 -36.9080883
[45,] -1.9282569 -24.3584065
[46,] -8.7562785 -1.9282569
[47,] -35.3054417 -8.7562785
[48,] 14.2896723 -35.3054417
[49,] -26.7980488 14.2896723
[50,] -10.9652484 -26.7980488
[51,] -13.4365501 -10.9652484
[52,] 12.9324280 -13.4365501
[53,] -0.5103915 12.9324280
[54,] 9.4165483 -0.5103915
[55,] -16.5968601 9.4165483
[56,] -9.1298389 -16.5968601
[57,] 12.4809332 -9.1298389
[58,] 5.6351639 12.4809332
[59,] 20.5181250 5.6351639
[60,] -9.6539777 20.5181250
[61,] 15.9131973 -9.6539777
[62,] 0.8716295 15.9131973
[63,] 28.3891603 0.8716295
[64,] 10.4736890 28.3891603
[65,] -1.6885947 10.4736890
[66,] -40.4977069 -1.6885947
[67,] 6.9774321 -40.4977069
[68,] -1.4669486 6.9774321
[69,] 2.0587226 -1.4669486
[70,] 6.4775761 2.0587226
[71,] 5.8732234 6.4775761
[72,] 5.7541048 5.8732234
[73,] 7.8861999 5.7541048
[74,] -29.1280082 7.8861999
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4.1267577 12.8441172
2 -14.1016066 4.1267577
3 0.1453587 -14.1016066
4 -21.8651680 0.1453587
5 10.1310756 -21.8651680
6 -6.4410527 10.1310756
7 -29.4475843 -6.4410527
8 -19.6341987 -29.4475843
9 8.0193363 -19.6341987
10 9.6708367 8.0193363
11 29.5955595 9.6708367
12 6.4693915 29.5955595
13 15.0914331 6.4693915
14 18.4392369 15.0914331
15 -2.5678956 18.4392369
16 -11.9295730 -2.5678956
17 1.7990599 -11.9295730
18 -61.3372354 1.7990599
19 2.3942481 -61.3372354
20 16.6646087 2.3942481
21 12.7207562 16.6646087
22 13.7992719 12.7207562
23 5.8772208 13.7992719
24 14.0097157 5.8772208
25 7.3189793 14.0097157
26 19.6237337 7.3189793
27 6.4770935 19.6237337
28 7.2862587 6.4770935
29 -17.0108922 7.2862587
30 -20.7086367 -17.0108922
31 20.7512135 -20.7086367
32 -6.5498416 20.7512135
33 17.4627388 -6.5498416
34 -15.4008674 17.4627388
35 21.4059495 -15.4008674
36 33.3404674 21.4059495
37 -9.3223167 33.3404674
38 -24.6115394 -9.3223167
39 8.2108208 -24.6115394
40 8.4916819 8.2108208
41 31.6047262 8.4916819
42 -1.6623990 31.6047262
43 -36.9080883 -1.6623990
44 -24.3584065 -36.9080883
45 -1.9282569 -24.3584065
46 -8.7562785 -1.9282569
47 -35.3054417 -8.7562785
48 14.2896723 -35.3054417
49 -26.7980488 14.2896723
50 -10.9652484 -26.7980488
51 -13.4365501 -10.9652484
52 12.9324280 -13.4365501
53 -0.5103915 12.9324280
54 9.4165483 -0.5103915
55 -16.5968601 9.4165483
56 -9.1298389 -16.5968601
57 12.4809332 -9.1298389
58 5.6351639 12.4809332
59 20.5181250 5.6351639
60 -9.6539777 20.5181250
61 15.9131973 -9.6539777
62 0.8716295 15.9131973
63 28.3891603 0.8716295
64 10.4736890 28.3891603
65 -1.6885947 10.4736890
66 -40.4977069 -1.6885947
67 6.9774321 -40.4977069
68 -1.4669486 6.9774321
69 2.0587226 -1.4669486
70 6.4775761 2.0587226
71 5.8732234 6.4775761
72 5.7541048 5.8732234
73 7.8861999 5.7541048
74 -29.1280082 7.8861999
> 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/7g4041324047132.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/8szd71324047132.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/9g9mw1324047132.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/10thku1324047132.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/118akw1324047132.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/1291et1324047132.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/13auiy1324047132.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/14wfjc1324047132.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/15kjhd1324047132.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/16jukz1324047132.tab")
+ }
>
> try(system("convert tmp/1egwf1324047132.ps tmp/1egwf1324047132.png",intern=TRUE))
character(0)
> try(system("convert tmp/2gkun1324047132.ps tmp/2gkun1324047132.png",intern=TRUE))
character(0)
> try(system("convert tmp/30rpq1324047132.ps tmp/30rpq1324047132.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ozgh1324047132.ps tmp/4ozgh1324047132.png",intern=TRUE))
character(0)
> try(system("convert tmp/5k9fo1324047132.ps tmp/5k9fo1324047132.png",intern=TRUE))
character(0)
> try(system("convert tmp/6n2wk1324047132.ps tmp/6n2wk1324047132.png",intern=TRUE))
character(0)
> try(system("convert tmp/7g4041324047132.ps tmp/7g4041324047132.png",intern=TRUE))
character(0)
> try(system("convert tmp/8szd71324047132.ps tmp/8szd71324047132.png",intern=TRUE))
character(0)
> try(system("convert tmp/9g9mw1324047132.ps tmp/9g9mw1324047132.png",intern=TRUE))
character(0)
> try(system("convert tmp/10thku1324047132.ps tmp/10thku1324047132.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.373 0.632 4.029