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(72772
+ ,26073
+ ,22274
+ ,45104
+ ,18103
+ ,14819
+ ,44525
+ ,15100
+ ,15136
+ ,41169
+ ,14738
+ ,13704
+ ,31118
+ ,22259
+ ,19638
+ ,28435
+ ,10277
+ ,7551
+ ,22162
+ ,6225
+ ,8019
+ ,20202
+ ,7663
+ ,6509
+ ,17773
+ ,6618
+ ,6634
+ ,17094
+ ,9945
+ ,11166
+ ,15153
+ ,7590
+ ,7508
+ ,11218
+ ,4293
+ ,4275
+ ,10796
+ ,4656
+ ,4944
+ ,9594
+ ,5145
+ ,5441
+ ,9309
+ ,2001
+ ,1689
+ ,8556
+ ,1779
+ ,1522
+ ,8041
+ ,1609
+ ,1416
+ ,7639
+ ,2191
+ ,1594
+ ,6884
+ ,1617
+ ,1909
+ ,6642
+ ,2554
+ ,2599
+ ,6321
+ ,2198
+ ,1262
+ ,6216
+ ,1578
+ ,1199
+ ,5865
+ ,3446
+ ,4404
+ ,5799
+ ,1380
+ ,1166
+ ,5695
+ ,1249
+ ,1122
+ ,5644
+ ,1223
+ ,886
+ ,5446
+ ,834
+ ,778
+ ,5395
+ ,3754
+ ,4436
+ ,5363
+ ,2283
+ ,1890
+ ,5338
+ ,3028
+ ,3107
+ ,5160
+ ,1100
+ ,1038
+ ,5091
+ ,457
+ ,300
+ ,5057
+ ,1201
+ ,988
+ ,5039
+ ,2192
+ ,2008
+ ,4880
+ ,1508
+ ,1522
+ ,4735
+ ,1393
+ ,1336
+ ,4693
+ ,952
+ ,976
+ ,4653
+ ,1032
+ ,798
+ ,4586
+ ,1279
+ ,869
+ ,4398
+ ,1370
+ ,1260
+ ,3974
+ ,649
+ ,578
+ ,3858
+ ,1900
+ ,2359
+ ,3826
+ ,666
+ ,736
+ ,3819
+ ,1313
+ ,1690
+ ,3556
+ ,1353
+ ,1201
+ ,3372
+ ,1500
+ ,813
+ ,3193
+ ,877
+ ,778
+ ,3126
+ ,874
+ ,687
+ ,3104
+ ,1133
+ ,1270
+ ,2967
+ ,754
+ ,671
+ ,2848
+ ,695
+ ,1559
+ ,2748
+ ,609
+ ,489
+ ,2649
+ ,696
+ ,773
+ ,2625
+ ,756
+ ,629
+ ,2572
+ ,670
+ ,637
+ ,2548
+ ,301
+ ,277
+ ,2477
+ ,630
+ ,776
+ ,2442
+ ,798
+ ,1651
+ ,2392
+ ,436
+ ,377
+ ,2372
+ ,388
+ ,222
+ ,2
+ ,346
+ ,864
+ ,1
+ ,068
+ ,2
+ ,251
+ ,497
+ ,399
+ ,2
+ ,230
+ ,449
+ ,547
+ ,2
+ ,225
+ ,919
+ ,668
+ ,2
+ ,220
+ ,536
+ ,451
+ ,2
+ ,205
+ ,673
+ ,724
+ ,2
+ ,193
+ ,837
+ ,853
+ ,2
+ ,116
+ ,534
+ ,434
+ ,2
+ ,102
+ ,845
+ ,730
+ ,2
+ ,099
+ ,626
+ ,612
+ ,2
+ ,096
+ ,871
+ ,558
+ ,2
+ ,064
+ ,740
+ ,859
+ ,2
+ ,036
+ ,391
+ ,311
+ ,1
+ ,920
+ ,435
+ ,318
+ ,1
+ ,813
+ ,424
+ ,312
+ ,1
+ ,776
+ ,338
+ ,343
+ ,1
+ ,752
+ ,744
+ ,710
+ ,1
+ ,738
+ ,368
+ ,273
+ ,1
+ ,729
+ ,393
+ ,259
+ ,1
+ ,685
+ ,938
+ ,1
+ ,274)
+ ,dim=c(3
+ ,80)
+ ,dimnames=list(c('weekdag'
+ ,'zaterdag'
+ ,'zondag')
+ ,1:80))
> y <- array(NA,dim=c(3,80),dimnames=list(c('weekdag','zaterdag','zondag'),1:80))
> 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 = '3'
> #'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
zondag weekdag zaterdag
1 22274 72772 26073
2 14819 45104 18103
3 15136 44525 15100
4 13704 41169 14738
5 19638 31118 22259
6 7551 28435 10277
7 8019 22162 6225
8 6509 20202 7663
9 6634 17773 6618
10 11166 17094 9945
11 7508 15153 7590
12 4275 11218 4293
13 4944 10796 4656
14 5441 9594 5145
15 1689 9309 2001
16 1522 8556 1779
17 1416 8041 1609
18 1594 7639 2191
19 1909 6884 1617
20 2599 6642 2554
21 1262 6321 2198
22 1199 6216 1578
23 4404 5865 3446
24 1166 5799 1380
25 1122 5695 1249
26 886 5644 1223
27 778 5446 834
28 4436 5395 3754
29 1890 5363 2283
30 3107 5338 3028
31 1038 5160 1100
32 300 5091 457
33 988 5057 1201
34 2008 5039 2192
35 1522 4880 1508
36 1336 4735 1393
37 976 4693 952
38 798 4653 1032
39 869 4586 1279
40 1260 4398 1370
41 578 3974 649
42 2359 3858 1900
43 736 3826 666
44 1690 3819 1313
45 1201 3556 1353
46 813 3372 1500
47 778 3193 877
48 687 3126 874
49 1270 3104 1133
50 671 2967 754
51 1559 2848 695
52 489 2748 609
53 773 2649 696
54 629 2625 756
55 637 2572 670
56 277 2548 301
57 776 2477 630
58 1651 2442 798
59 377 2392 436
60 222 2372 388
61 864 2 346
62 2 1 68
63 399 251 497
64 449 2 230
65 225 547 2
66 2 919 668
67 451 220 536
68 673 2 205
69 193 724 2
70 2 837 853
71 434 116 534
72 845 2 102
73 99 730 2
74 2 626 612
75 558 96 871
76 740 2 64
77 36 859 2
78 1 391 311
79 318 920 435
80 424 1 813
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) weekdag zaterdag
156.77978 -0.01221 0.91848
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1697.8 -238.2 -116.4 253.7 2415.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 156.77978 84.67367 1.852 0.0679 .
weekdag -0.01221 0.01999 -0.611 0.5433
zaterdag 0.91848 0.04912 18.697 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 637.8 on 77 degrees of freedom
Multiple R-squared: 0.9794, Adjusted R-squared: 0.9789
F-statistic: 1834 on 2 and 77 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.9999982 3.632411e-06 1.816205e-06
[2,] 1.0000000 1.238069e-09 6.190343e-10
[3,] 1.0000000 1.225065e-10 6.125323e-11
[4,] 1.0000000 5.382636e-10 2.691318e-10
[5,] 1.0000000 1.384257e-11 6.921285e-12
[6,] 1.0000000 4.729344e-11 2.364672e-11
[7,] 1.0000000 9.887104e-11 4.943552e-11
[8,] 1.0000000 2.953547e-10 1.476773e-10
[9,] 1.0000000 9.236675e-10 4.618338e-10
[10,] 1.0000000 6.948003e-10 3.474002e-10
[11,] 1.0000000 8.913509e-10 4.456754e-10
[12,] 1.0000000 1.559364e-09 7.796822e-10
[13,] 1.0000000 1.189942e-09 5.949712e-10
[14,] 1.0000000 2.450726e-09 1.225363e-09
[15,] 1.0000000 6.705830e-09 3.352915e-09
[16,] 1.0000000 9.954338e-10 4.977169e-10
[17,] 1.0000000 1.681943e-09 8.409714e-10
[18,] 1.0000000 3.949030e-10 1.974515e-10
[19,] 1.0000000 9.155311e-10 4.577655e-10
[20,] 1.0000000 2.330791e-09 1.165395e-09
[21,] 1.0000000 4.243855e-09 2.121928e-09
[22,] 1.0000000 1.096646e-08 5.483228e-09
[23,] 1.0000000 2.906871e-09 1.453435e-09
[24,] 1.0000000 5.988361e-09 2.994181e-09
[25,] 1.0000000 9.239914e-09 4.619957e-09
[26,] 1.0000000 2.356780e-08 1.178390e-08
[27,] 1.0000000 4.179407e-08 2.089703e-08
[28,] 1.0000000 8.817012e-08 4.408506e-08
[29,] 0.9999999 2.033330e-07 1.016665e-07
[30,] 0.9999998 4.745295e-07 2.372647e-07
[31,] 0.9999994 1.113057e-06 5.565287e-07
[32,] 0.9999987 2.583099e-06 1.291550e-06
[33,] 0.9999977 4.678826e-06 2.339413e-06
[34,] 0.9999967 6.594421e-06 3.297210e-06
[35,] 0.9999929 1.426049e-05 7.130243e-06
[36,] 0.9999870 2.607025e-05 1.303513e-05
[37,] 0.9999951 9.895365e-06 4.947683e-06
[38,] 0.9999894 2.115227e-05 1.057614e-05
[39,] 0.9999910 1.806752e-05 9.033760e-06
[40,] 0.9999824 3.511704e-05 1.755852e-05
[41,] 0.9999742 5.153639e-05 2.576819e-05
[42,] 0.9999464 1.072384e-04 5.361919e-05
[43,] 0.9998979 2.041151e-04 1.020575e-04
[44,] 0.9998421 3.158965e-04 1.579483e-04
[45,] 0.9996903 6.193045e-04 3.096523e-04
[46,] 0.9999218 1.563731e-04 7.818657e-05
[47,] 0.9998410 3.180919e-04 1.590460e-04
[48,] 0.9996929 6.142802e-04 3.071401e-04
[49,] 0.9993890 1.221959e-03 6.109797e-04
[50,] 0.9988150 2.370072e-03 1.185036e-03
[51,] 0.9979374 4.125181e-03 2.062591e-03
[52,] 0.9964871 7.025707e-03 3.512853e-03
[53,] 0.9999825 3.498669e-05 1.749335e-05
[54,] 0.9999845 3.107265e-05 1.553633e-05
[55,] 0.9999980 4.048100e-06 2.024050e-06
[56,] 0.9999978 4.341156e-06 2.170578e-06
[57,] 0.9999999 2.176370e-07 1.088185e-07
[58,] 0.9999995 9.424474e-07 4.712237e-07
[59,] 0.9999987 2.672564e-06 1.336282e-06
[60,] 0.9999949 1.020246e-05 5.101230e-06
[61,] 0.9999825 3.508239e-05 1.754120e-05
[62,] 0.9999325 1.350889e-04 6.754444e-05
[63,] 0.9997577 4.845111e-04 2.422556e-04
[64,] 0.9991435 1.713087e-03 8.565433e-04
[65,] 0.9974801 5.039876e-03 2.519938e-03
[66,] 0.9922653 1.546942e-02 7.734710e-03
[67,] 0.9854110 2.917808e-02 1.458904e-02
[68,] 0.9572557 8.548854e-02 4.274427e-02
[69,] 0.9067746 1.864507e-01 9.322536e-02
> postscript(file="/var/wessaorg/rcomp/tmp/1wvek1322130626.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/2ap6v1322130626.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/3fmy61322130626.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/4e21i1322130626.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/5ffh81322130626.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 = 80
Frequency = 1
1 2 3 4 5 6
-941.839399 -1414.344288 1653.772632 513.292037 -583.270203 -1697.837056
7 8 9 10 11 12
2415.251027 -439.445443 615.710031 2083.649381 564.966553 312.146470
13 14 15 16 17 18
642.587809 675.779075 -192.009736 -164.300349 -120.446307 -481.907186
19 20 21 22 23 24
351.081544 177.514713 -836.426301 -331.252628 1153.748196 -187.484884
25 26 27 28 29 30
-112.434055 -325.176257 -78.305994 897.119777 -298.191810 234.237939
31 32 33 34 35 36
-66.112175 -214.374062 -210.135706 -100.565743 39.731215 -42.414098
37 38 39 40 41 42
2.121353 -249.845083 -406.526717 -101.403139 -126.357594 504.212099
43 44 45 46 47 48
14.221561 373.881755 -155.067940 -680.330219 -145.304488 -234.366976
49 50 51 52 53 54
110.479018 -142.090813 798.646586 -193.585200 9.298773 -190.102808
55 56 57 58 59 60
-103.760831 -125.135955 70.818499 791.087162 -151.034699 -262.191978
61 62 63 64 65 66
389.451735 -217.223982 -211.198507 80.995019 73.060890 -757.103244
67 68 69 70 71 72
-195.397533 327.956933 43.221657 -928.022446 -211.830183 594.560022
73 74 75 76 77 78
-50.705097 -709.245418 -397.600948 524.462132 -112.130300 -436.652780
79 80
-227.085992 -479.489039
> postscript(file="/var/wessaorg/rcomp/tmp/6et881322130626.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 = 80
Frequency = 1
lag(myerror, k = 1) myerror
0 -941.839399 NA
1 -1414.344288 -941.839399
2 1653.772632 -1414.344288
3 513.292037 1653.772632
4 -583.270203 513.292037
5 -1697.837056 -583.270203
6 2415.251027 -1697.837056
7 -439.445443 2415.251027
8 615.710031 -439.445443
9 2083.649381 615.710031
10 564.966553 2083.649381
11 312.146470 564.966553
12 642.587809 312.146470
13 675.779075 642.587809
14 -192.009736 675.779075
15 -164.300349 -192.009736
16 -120.446307 -164.300349
17 -481.907186 -120.446307
18 351.081544 -481.907186
19 177.514713 351.081544
20 -836.426301 177.514713
21 -331.252628 -836.426301
22 1153.748196 -331.252628
23 -187.484884 1153.748196
24 -112.434055 -187.484884
25 -325.176257 -112.434055
26 -78.305994 -325.176257
27 897.119777 -78.305994
28 -298.191810 897.119777
29 234.237939 -298.191810
30 -66.112175 234.237939
31 -214.374062 -66.112175
32 -210.135706 -214.374062
33 -100.565743 -210.135706
34 39.731215 -100.565743
35 -42.414098 39.731215
36 2.121353 -42.414098
37 -249.845083 2.121353
38 -406.526717 -249.845083
39 -101.403139 -406.526717
40 -126.357594 -101.403139
41 504.212099 -126.357594
42 14.221561 504.212099
43 373.881755 14.221561
44 -155.067940 373.881755
45 -680.330219 -155.067940
46 -145.304488 -680.330219
47 -234.366976 -145.304488
48 110.479018 -234.366976
49 -142.090813 110.479018
50 798.646586 -142.090813
51 -193.585200 798.646586
52 9.298773 -193.585200
53 -190.102808 9.298773
54 -103.760831 -190.102808
55 -125.135955 -103.760831
56 70.818499 -125.135955
57 791.087162 70.818499
58 -151.034699 791.087162
59 -262.191978 -151.034699
60 389.451735 -262.191978
61 -217.223982 389.451735
62 -211.198507 -217.223982
63 80.995019 -211.198507
64 73.060890 80.995019
65 -757.103244 73.060890
66 -195.397533 -757.103244
67 327.956933 -195.397533
68 43.221657 327.956933
69 -928.022446 43.221657
70 -211.830183 -928.022446
71 594.560022 -211.830183
72 -50.705097 594.560022
73 -709.245418 -50.705097
74 -397.600948 -709.245418
75 524.462132 -397.600948
76 -112.130300 524.462132
77 -436.652780 -112.130300
78 -227.085992 -436.652780
79 -479.489039 -227.085992
80 NA -479.489039
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1414.344288 -941.839399
[2,] 1653.772632 -1414.344288
[3,] 513.292037 1653.772632
[4,] -583.270203 513.292037
[5,] -1697.837056 -583.270203
[6,] 2415.251027 -1697.837056
[7,] -439.445443 2415.251027
[8,] 615.710031 -439.445443
[9,] 2083.649381 615.710031
[10,] 564.966553 2083.649381
[11,] 312.146470 564.966553
[12,] 642.587809 312.146470
[13,] 675.779075 642.587809
[14,] -192.009736 675.779075
[15,] -164.300349 -192.009736
[16,] -120.446307 -164.300349
[17,] -481.907186 -120.446307
[18,] 351.081544 -481.907186
[19,] 177.514713 351.081544
[20,] -836.426301 177.514713
[21,] -331.252628 -836.426301
[22,] 1153.748196 -331.252628
[23,] -187.484884 1153.748196
[24,] -112.434055 -187.484884
[25,] -325.176257 -112.434055
[26,] -78.305994 -325.176257
[27,] 897.119777 -78.305994
[28,] -298.191810 897.119777
[29,] 234.237939 -298.191810
[30,] -66.112175 234.237939
[31,] -214.374062 -66.112175
[32,] -210.135706 -214.374062
[33,] -100.565743 -210.135706
[34,] 39.731215 -100.565743
[35,] -42.414098 39.731215
[36,] 2.121353 -42.414098
[37,] -249.845083 2.121353
[38,] -406.526717 -249.845083
[39,] -101.403139 -406.526717
[40,] -126.357594 -101.403139
[41,] 504.212099 -126.357594
[42,] 14.221561 504.212099
[43,] 373.881755 14.221561
[44,] -155.067940 373.881755
[45,] -680.330219 -155.067940
[46,] -145.304488 -680.330219
[47,] -234.366976 -145.304488
[48,] 110.479018 -234.366976
[49,] -142.090813 110.479018
[50,] 798.646586 -142.090813
[51,] -193.585200 798.646586
[52,] 9.298773 -193.585200
[53,] -190.102808 9.298773
[54,] -103.760831 -190.102808
[55,] -125.135955 -103.760831
[56,] 70.818499 -125.135955
[57,] 791.087162 70.818499
[58,] -151.034699 791.087162
[59,] -262.191978 -151.034699
[60,] 389.451735 -262.191978
[61,] -217.223982 389.451735
[62,] -211.198507 -217.223982
[63,] 80.995019 -211.198507
[64,] 73.060890 80.995019
[65,] -757.103244 73.060890
[66,] -195.397533 -757.103244
[67,] 327.956933 -195.397533
[68,] 43.221657 327.956933
[69,] -928.022446 43.221657
[70,] -211.830183 -928.022446
[71,] 594.560022 -211.830183
[72,] -50.705097 594.560022
[73,] -709.245418 -50.705097
[74,] -397.600948 -709.245418
[75,] 524.462132 -397.600948
[76,] -112.130300 524.462132
[77,] -436.652780 -112.130300
[78,] -227.085992 -436.652780
[79,] -479.489039 -227.085992
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1414.344288 -941.839399
2 1653.772632 -1414.344288
3 513.292037 1653.772632
4 -583.270203 513.292037
5 -1697.837056 -583.270203
6 2415.251027 -1697.837056
7 -439.445443 2415.251027
8 615.710031 -439.445443
9 2083.649381 615.710031
10 564.966553 2083.649381
11 312.146470 564.966553
12 642.587809 312.146470
13 675.779075 642.587809
14 -192.009736 675.779075
15 -164.300349 -192.009736
16 -120.446307 -164.300349
17 -481.907186 -120.446307
18 351.081544 -481.907186
19 177.514713 351.081544
20 -836.426301 177.514713
21 -331.252628 -836.426301
22 1153.748196 -331.252628
23 -187.484884 1153.748196
24 -112.434055 -187.484884
25 -325.176257 -112.434055
26 -78.305994 -325.176257
27 897.119777 -78.305994
28 -298.191810 897.119777
29 234.237939 -298.191810
30 -66.112175 234.237939
31 -214.374062 -66.112175
32 -210.135706 -214.374062
33 -100.565743 -210.135706
34 39.731215 -100.565743
35 -42.414098 39.731215
36 2.121353 -42.414098
37 -249.845083 2.121353
38 -406.526717 -249.845083
39 -101.403139 -406.526717
40 -126.357594 -101.403139
41 504.212099 -126.357594
42 14.221561 504.212099
43 373.881755 14.221561
44 -155.067940 373.881755
45 -680.330219 -155.067940
46 -145.304488 -680.330219
47 -234.366976 -145.304488
48 110.479018 -234.366976
49 -142.090813 110.479018
50 798.646586 -142.090813
51 -193.585200 798.646586
52 9.298773 -193.585200
53 -190.102808 9.298773
54 -103.760831 -190.102808
55 -125.135955 -103.760831
56 70.818499 -125.135955
57 791.087162 70.818499
58 -151.034699 791.087162
59 -262.191978 -151.034699
60 389.451735 -262.191978
61 -217.223982 389.451735
62 -211.198507 -217.223982
63 80.995019 -211.198507
64 73.060890 80.995019
65 -757.103244 73.060890
66 -195.397533 -757.103244
67 327.956933 -195.397533
68 43.221657 327.956933
69 -928.022446 43.221657
70 -211.830183 -928.022446
71 594.560022 -211.830183
72 -50.705097 594.560022
73 -709.245418 -50.705097
74 -397.600948 -709.245418
75 524.462132 -397.600948
76 -112.130300 524.462132
77 -436.652780 -112.130300
78 -227.085992 -436.652780
79 -479.489039 -227.085992
> 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/7rckf1322130626.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/8b8d71322130626.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/94sse1322130626.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/10izax1322130626.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/11l8401322130626.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/12b6iq1322130626.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/13lzy01322130626.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/14s68k1322130626.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/15jt0v1322130626.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/161d261322130626.tab")
+ }
>
> try(system("convert tmp/1wvek1322130626.ps tmp/1wvek1322130626.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ap6v1322130626.ps tmp/2ap6v1322130626.png",intern=TRUE))
character(0)
> try(system("convert tmp/3fmy61322130626.ps tmp/3fmy61322130626.png",intern=TRUE))
character(0)
> try(system("convert tmp/4e21i1322130626.ps tmp/4e21i1322130626.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ffh81322130626.ps tmp/5ffh81322130626.png",intern=TRUE))
character(0)
> try(system("convert tmp/6et881322130626.ps tmp/6et881322130626.png",intern=TRUE))
character(0)
> try(system("convert tmp/7rckf1322130626.ps tmp/7rckf1322130626.png",intern=TRUE))
character(0)
> try(system("convert tmp/8b8d71322130626.ps tmp/8b8d71322130626.png",intern=TRUE))
character(0)
> try(system("convert tmp/94sse1322130626.ps tmp/94sse1322130626.png",intern=TRUE))
character(0)
> try(system("convert tmp/10izax1322130626.ps tmp/10izax1322130626.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.368 0.504 3.904