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)
+ ,dim=c(3
+ ,70)
+ ,dimnames=list(c('weekdag'
+ ,'zaterdag'
+ ,'zondag')
+ ,1:70))
> y <- array(NA,dim=c(3,70),dimnames=list(c('weekdag','zaterdag','zondag'),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 = '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
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) weekdag zaterdag
189.14252 -0.01633 0.92617
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1692.2 -249.2 -130.8 269.9 2426.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 189.14252 96.87310 1.952 0.0551 .
weekdag -0.01633 0.02112 -0.773 0.4421
zaterdag 0.92617 0.05148 17.991 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 663 on 67 degrees of freedom
Multiple R-squared: 0.9798, Adjusted R-squared: 0.9792
F-statistic: 1628 on 2 and 67 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.9999968 6.409407e-06 3.204703e-06
[2,] 1.0000000 8.943439e-10 4.471719e-10
[3,] 1.0000000 8.751355e-11 4.375677e-11
[4,] 1.0000000 4.015335e-10 2.007668e-10
[5,] 1.0000000 1.431640e-11 7.158198e-12
[6,] 1.0000000 4.661200e-11 2.330600e-11
[7,] 1.0000000 1.064383e-10 5.321917e-11
[8,] 1.0000000 3.359736e-10 1.679868e-10
[9,] 1.0000000 1.155604e-09 5.778019e-10
[10,] 1.0000000 1.027479e-09 5.137395e-10
[11,] 1.0000000 1.483038e-09 7.415189e-10
[12,] 1.0000000 2.782665e-09 1.391332e-09
[13,] 1.0000000 2.388310e-09 1.194155e-09
[14,] 1.0000000 4.300365e-09 2.150183e-09
[15,] 1.0000000 1.235582e-08 6.177912e-09
[16,] 1.0000000 1.547261e-09 7.736305e-10
[17,] 1.0000000 2.880257e-09 1.440128e-09
[18,] 1.0000000 8.670794e-10 4.335397e-10
[19,] 1.0000000 2.203537e-09 1.101768e-09
[20,] 1.0000000 6.034924e-09 3.017462e-09
[21,] 1.0000000 1.225684e-08 6.128418e-09
[22,] 1.0000000 3.310484e-08 1.655242e-08
[23,] 1.0000000 1.573387e-08 7.866934e-09
[24,] 1.0000000 3.095890e-08 1.547945e-08
[25,] 1.0000000 6.893895e-08 3.446948e-08
[26,] 0.9999999 1.827843e-07 9.139217e-08
[27,] 0.9999998 3.969889e-07 1.984945e-07
[28,] 0.9999996 8.814658e-07 4.407329e-07
[29,] 0.9999989 2.179127e-06 1.089563e-06
[30,] 0.9999974 5.249780e-06 2.624890e-06
[31,] 0.9999938 1.247743e-05 6.238717e-06
[32,] 0.9999855 2.903870e-05 1.451935e-05
[33,] 0.9999726 5.481128e-05 2.740564e-05
[34,] 0.9999618 7.638382e-05 3.819191e-05
[35,] 0.9999187 1.625281e-04 8.126405e-05
[36,] 0.9998436 3.127077e-04 1.563539e-04
[37,] 0.9999046 1.907062e-04 9.535311e-05
[38,] 0.9997992 4.016158e-04 2.008079e-04
[39,] 0.9998005 3.989032e-04 1.994516e-04
[40,] 0.9996011 7.977236e-04 3.988618e-04
[41,] 0.9994516 1.096892e-03 5.484459e-04
[42,] 0.9989047 2.190677e-03 1.095338e-03
[43,] 0.9980372 3.925512e-03 1.962756e-03
[44,] 0.9966050 6.789951e-03 3.394975e-03
[45,] 0.9937679 1.246422e-02 6.232112e-03
[46,] 0.9976120 4.776067e-03 2.388033e-03
[47,] 0.9954024 9.195119e-03 4.597560e-03
[48,] 0.9911441 1.771183e-02 8.855913e-03
[49,] 0.9833772 3.324553e-02 1.662277e-02
[50,] 0.9696582 6.068356e-02 3.034178e-02
[51,] 0.9522379 9.552416e-02 4.776208e-02
[52,] 0.9200716 1.598568e-01 7.992842e-02
[53,] 0.9968970 6.205985e-03 3.102993e-03
[54,] 0.9945834 1.083318e-02 5.416591e-03
[55,] 0.9963034 7.393166e-03 3.696583e-03
[56,] 0.9965493 6.901373e-03 3.450686e-03
[57,] 0.9999454 1.091665e-04 5.458325e-05
[58,] 0.9994883 1.023483e-03 5.117414e-04
[59,] 0.9995410 9.180165e-04 4.590083e-04
> postscript(file="/var/wessaorg/rcomp/tmp/1u3uk1322130531.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/2p0rd1322130531.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/3wbrw1322130531.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/43rtd1322130531.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/5ptwu1322130531.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
-875.120042 -1400.252884 1688.595480 537.076279 -658.787212 -1692.171098
7 8 9 10 11 12
2426.266730 -447.573591 605.619742 2045.151139 536.600709 292.950010
13 14 15 16 17 18
618.858574 643.333880 -201.426938 -175.110652 -132.069546 -499.666639
19 20 21 22 23 24
334.630396 152.853739 -859.669215 -350.155423 1119.019838 -206.581359
25 26 27 28 29 30
-130.950543 -343.702699 -94.653641 858.084287 -326.035577 200.556279
31 32 33 34 35 36
-85.685646 -229.282059 -230.910974 -129.043741 15.863542 -65.993852
37 38 39 40 41 42
-18.236667 -270.983714 -429.842730 -126.194133 -147.345124 473.116680
43 44 45 46 47 48
-7.506527 347.144316 -183.196737 -710.348600 -171.264503 -260.579908
49 50 51 52 53 54
82.181711 -168.035012 772.666334 -219.315391 -17.508969 -217.471290
55 56 57 58 59 60
-130.685634 -149.319118 43.810253 762.641491 -176.899724 -287.769896
61 62 63 64 65 66
354.433778 -250.106053 -246.353071 46.870014 42.936155 -790.822278
67 68 69 70
-230.980020 294.024375 13.826083 -963.503387
> postscript(file="/var/wessaorg/rcomp/tmp/67v2m1322130531.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 -875.120042 NA
1 -1400.252884 -875.120042
2 1688.595480 -1400.252884
3 537.076279 1688.595480
4 -658.787212 537.076279
5 -1692.171098 -658.787212
6 2426.266730 -1692.171098
7 -447.573591 2426.266730
8 605.619742 -447.573591
9 2045.151139 605.619742
10 536.600709 2045.151139
11 292.950010 536.600709
12 618.858574 292.950010
13 643.333880 618.858574
14 -201.426938 643.333880
15 -175.110652 -201.426938
16 -132.069546 -175.110652
17 -499.666639 -132.069546
18 334.630396 -499.666639
19 152.853739 334.630396
20 -859.669215 152.853739
21 -350.155423 -859.669215
22 1119.019838 -350.155423
23 -206.581359 1119.019838
24 -130.950543 -206.581359
25 -343.702699 -130.950543
26 -94.653641 -343.702699
27 858.084287 -94.653641
28 -326.035577 858.084287
29 200.556279 -326.035577
30 -85.685646 200.556279
31 -229.282059 -85.685646
32 -230.910974 -229.282059
33 -129.043741 -230.910974
34 15.863542 -129.043741
35 -65.993852 15.863542
36 -18.236667 -65.993852
37 -270.983714 -18.236667
38 -429.842730 -270.983714
39 -126.194133 -429.842730
40 -147.345124 -126.194133
41 473.116680 -147.345124
42 -7.506527 473.116680
43 347.144316 -7.506527
44 -183.196737 347.144316
45 -710.348600 -183.196737
46 -171.264503 -710.348600
47 -260.579908 -171.264503
48 82.181711 -260.579908
49 -168.035012 82.181711
50 772.666334 -168.035012
51 -219.315391 772.666334
52 -17.508969 -219.315391
53 -217.471290 -17.508969
54 -130.685634 -217.471290
55 -149.319118 -130.685634
56 43.810253 -149.319118
57 762.641491 43.810253
58 -176.899724 762.641491
59 -287.769896 -176.899724
60 354.433778 -287.769896
61 -250.106053 354.433778
62 -246.353071 -250.106053
63 46.870014 -246.353071
64 42.936155 46.870014
65 -790.822278 42.936155
66 -230.980020 -790.822278
67 294.024375 -230.980020
68 13.826083 294.024375
69 -963.503387 13.826083
70 NA -963.503387
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1400.252884 -875.120042
[2,] 1688.595480 -1400.252884
[3,] 537.076279 1688.595480
[4,] -658.787212 537.076279
[5,] -1692.171098 -658.787212
[6,] 2426.266730 -1692.171098
[7,] -447.573591 2426.266730
[8,] 605.619742 -447.573591
[9,] 2045.151139 605.619742
[10,] 536.600709 2045.151139
[11,] 292.950010 536.600709
[12,] 618.858574 292.950010
[13,] 643.333880 618.858574
[14,] -201.426938 643.333880
[15,] -175.110652 -201.426938
[16,] -132.069546 -175.110652
[17,] -499.666639 -132.069546
[18,] 334.630396 -499.666639
[19,] 152.853739 334.630396
[20,] -859.669215 152.853739
[21,] -350.155423 -859.669215
[22,] 1119.019838 -350.155423
[23,] -206.581359 1119.019838
[24,] -130.950543 -206.581359
[25,] -343.702699 -130.950543
[26,] -94.653641 -343.702699
[27,] 858.084287 -94.653641
[28,] -326.035577 858.084287
[29,] 200.556279 -326.035577
[30,] -85.685646 200.556279
[31,] -229.282059 -85.685646
[32,] -230.910974 -229.282059
[33,] -129.043741 -230.910974
[34,] 15.863542 -129.043741
[35,] -65.993852 15.863542
[36,] -18.236667 -65.993852
[37,] -270.983714 -18.236667
[38,] -429.842730 -270.983714
[39,] -126.194133 -429.842730
[40,] -147.345124 -126.194133
[41,] 473.116680 -147.345124
[42,] -7.506527 473.116680
[43,] 347.144316 -7.506527
[44,] -183.196737 347.144316
[45,] -710.348600 -183.196737
[46,] -171.264503 -710.348600
[47,] -260.579908 -171.264503
[48,] 82.181711 -260.579908
[49,] -168.035012 82.181711
[50,] 772.666334 -168.035012
[51,] -219.315391 772.666334
[52,] -17.508969 -219.315391
[53,] -217.471290 -17.508969
[54,] -130.685634 -217.471290
[55,] -149.319118 -130.685634
[56,] 43.810253 -149.319118
[57,] 762.641491 43.810253
[58,] -176.899724 762.641491
[59,] -287.769896 -176.899724
[60,] 354.433778 -287.769896
[61,] -250.106053 354.433778
[62,] -246.353071 -250.106053
[63,] 46.870014 -246.353071
[64,] 42.936155 46.870014
[65,] -790.822278 42.936155
[66,] -230.980020 -790.822278
[67,] 294.024375 -230.980020
[68,] 13.826083 294.024375
[69,] -963.503387 13.826083
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1400.252884 -875.120042
2 1688.595480 -1400.252884
3 537.076279 1688.595480
4 -658.787212 537.076279
5 -1692.171098 -658.787212
6 2426.266730 -1692.171098
7 -447.573591 2426.266730
8 605.619742 -447.573591
9 2045.151139 605.619742
10 536.600709 2045.151139
11 292.950010 536.600709
12 618.858574 292.950010
13 643.333880 618.858574
14 -201.426938 643.333880
15 -175.110652 -201.426938
16 -132.069546 -175.110652
17 -499.666639 -132.069546
18 334.630396 -499.666639
19 152.853739 334.630396
20 -859.669215 152.853739
21 -350.155423 -859.669215
22 1119.019838 -350.155423
23 -206.581359 1119.019838
24 -130.950543 -206.581359
25 -343.702699 -130.950543
26 -94.653641 -343.702699
27 858.084287 -94.653641
28 -326.035577 858.084287
29 200.556279 -326.035577
30 -85.685646 200.556279
31 -229.282059 -85.685646
32 -230.910974 -229.282059
33 -129.043741 -230.910974
34 15.863542 -129.043741
35 -65.993852 15.863542
36 -18.236667 -65.993852
37 -270.983714 -18.236667
38 -429.842730 -270.983714
39 -126.194133 -429.842730
40 -147.345124 -126.194133
41 473.116680 -147.345124
42 -7.506527 473.116680
43 347.144316 -7.506527
44 -183.196737 347.144316
45 -710.348600 -183.196737
46 -171.264503 -710.348600
47 -260.579908 -171.264503
48 82.181711 -260.579908
49 -168.035012 82.181711
50 772.666334 -168.035012
51 -219.315391 772.666334
52 -17.508969 -219.315391
53 -217.471290 -17.508969
54 -130.685634 -217.471290
55 -149.319118 -130.685634
56 43.810253 -149.319118
57 762.641491 43.810253
58 -176.899724 762.641491
59 -287.769896 -176.899724
60 354.433778 -287.769896
61 -250.106053 354.433778
62 -246.353071 -250.106053
63 46.870014 -246.353071
64 42.936155 46.870014
65 -790.822278 42.936155
66 -230.980020 -790.822278
67 294.024375 -230.980020
68 13.826083 294.024375
69 -963.503387 13.826083
> 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/7sia51322130531.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/8j1xd1322130531.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/94ryt1322130531.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/10rp6o1322130531.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/110yl41322130531.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/12pfif1322130531.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/138e9f1322130531.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/14ac121322130531.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/15vc8e1322130532.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/16d5rh1322130532.tab")
+ }
>
> try(system("convert tmp/1u3uk1322130531.ps tmp/1u3uk1322130531.png",intern=TRUE))
character(0)
> try(system("convert tmp/2p0rd1322130531.ps tmp/2p0rd1322130531.png",intern=TRUE))
character(0)
> try(system("convert tmp/3wbrw1322130531.ps tmp/3wbrw1322130531.png",intern=TRUE))
character(0)
> try(system("convert tmp/43rtd1322130531.ps tmp/43rtd1322130531.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ptwu1322130531.ps tmp/5ptwu1322130531.png",intern=TRUE))
character(0)
> try(system("convert tmp/67v2m1322130531.ps tmp/67v2m1322130531.png",intern=TRUE))
character(0)
> try(system("convert tmp/7sia51322130531.ps tmp/7sia51322130531.png",intern=TRUE))
character(0)
> try(system("convert tmp/8j1xd1322130531.ps tmp/8j1xd1322130531.png",intern=TRUE))
character(0)
> try(system("convert tmp/94ryt1322130531.ps tmp/94ryt1322130531.png",intern=TRUE))
character(0)
> try(system("convert tmp/10rp6o1322130531.ps tmp/10rp6o1322130531.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.151 0.492 3.674