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
+ ,2193
+ ,837
+ ,853
+ ,2116
+ ,534
+ ,434
+ ,2102
+ ,845
+ ,730
+ ,2099
+ ,626
+ ,612
+ ,2096
+ ,871
+ ,558
+ ,2064
+ ,740
+ ,859
+ ,2036
+ ,391
+ ,311
+ ,1920
+ ,435
+ ,318
+ ,1813
+ ,424
+ ,312
+ ,1776
+ ,338
+ ,343
+ ,1752
+ ,744
+ ,710
+ ,1738
+ ,368
+ ,273
+ ,1729
+ ,393
+ ,259
+ ,1685
+ ,938
+ ,1274)
+ ,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 = 'Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '3'
> 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 t
1 22274 72772 26073 1
2 14819 45104 18103 2
3 15136 44525 15100 3
4 13704 41169 14738 4
5 19638 31118 22259 5
6 7551 28435 10277 6
7 8019 22162 6225 7
8 6509 20202 7663 8
9 6634 17773 6618 9
10 11166 17094 9945 10
11 7508 15153 7590 11
12 4275 11218 4293 12
13 4944 10796 4656 13
14 5441 9594 5145 14
15 1689 9309 2001 15
16 1522 8556 1779 16
17 1416 8041 1609 17
18 1594 7639 2191 18
19 1909 6884 1617 19
20 2599 6642 2554 20
21 1262 6321 2198 21
22 1199 6216 1578 22
23 4404 5865 3446 23
24 1166 5799 1380 24
25 1122 5695 1249 25
26 886 5644 1223 26
27 778 5446 834 27
28 4436 5395 3754 28
29 1890 5363 2283 29
30 3107 5338 3028 30
31 1038 5160 1100 31
32 300 5091 457 32
33 988 5057 1201 33
34 2008 5039 2192 34
35 1522 4880 1508 35
36 1336 4735 1393 36
37 976 4693 952 37
38 798 4653 1032 38
39 869 4586 1279 39
40 1260 4398 1370 40
41 578 3974 649 41
42 2359 3858 1900 42
43 736 3826 666 43
44 1690 3819 1313 44
45 1201 3556 1353 45
46 813 3372 1500 46
47 778 3193 877 47
48 687 3126 874 48
49 1270 3104 1133 49
50 671 2967 754 50
51 1559 2848 695 51
52 489 2748 609 52
53 773 2649 696 53
54 629 2625 756 54
55 637 2572 670 55
56 277 2548 301 56
57 776 2477 630 57
58 1651 2442 798 58
59 377 2392 436 59
60 222 2372 388 60
61 864 2 346 61
62 2 1 68 62
63 399 251 497 63
64 449 2 230 64
65 225 547 2 65
66 2 919 668 66
67 451 220 536 67
68 673 2 205 68
69 837 724 2193 69
70 534 853 2116 70
71 845 434 2102 71
72 626 730 2099 72
73 871 612 2096 73
74 740 558 2064 74
75 391 859 2036 75
76 435 311 1920 76
77 424 318 1813 77
78 338 312 1776 78
79 744 343 1752 79
80 368 710 1738 80
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) weekdag zaterdag t
1.056e+03 -7.683e-03 8.513e-01 -2.226e+01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1902.09 -421.85 -40.18 324.36 1997.46
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.056e+03 2.467e+02 4.280 5.38e-05 ***
weekdag -7.683e-03 2.156e-02 -0.356 0.723
zaterdag 8.513e-01 4.873e-02 17.471 < 2e-16 ***
t -2.226e+01 4.589e+00 -4.851 6.39e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 673.5 on 76 degrees of freedom
Multiple R-squared: 0.9771, Adjusted R-squared: 0.9762
F-statistic: 1082 on 3 and 76 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.9999971 5.782119e-06 2.891059e-06
[2,] 0.9999982 3.685634e-06 1.842817e-06
[3,] 0.9999940 1.208762e-05 6.043811e-06
[4,] 0.9999997 5.282040e-07 2.641020e-07
[5,] 0.9999998 4.774745e-07 2.387373e-07
[6,] 0.9999997 5.327464e-07 2.663732e-07
[7,] 0.9999998 4.210953e-07 2.105477e-07
[8,] 0.9999999 1.404914e-07 7.024569e-08
[9,] 0.9999999 1.377630e-07 6.888150e-08
[10,] 0.9999999 2.729676e-07 1.364838e-07
[11,] 0.9999997 6.799480e-07 3.399740e-07
[12,] 0.9999994 1.295636e-06 6.478178e-07
[13,] 0.9999987 2.587995e-06 1.293997e-06
[14,] 0.9999973 5.497919e-06 2.748959e-06
[15,] 0.9999986 2.842435e-06 1.421217e-06
[16,] 0.9999982 3.657009e-06 1.828504e-06
[17,] 0.9999999 2.016109e-07 1.008055e-07
[18,] 0.9999998 3.165122e-07 1.582561e-07
[19,] 0.9999997 5.058032e-07 2.529016e-07
[20,] 0.9999998 4.966992e-07 2.483496e-07
[21,] 0.9999997 5.746098e-07 2.873049e-07
[22,] 1.0000000 1.302733e-08 6.513666e-09
[23,] 1.0000000 3.244602e-08 1.622301e-08
[24,] 1.0000000 1.462817e-08 7.314087e-09
[25,] 1.0000000 3.124110e-08 1.562055e-08
[26,] 1.0000000 1.570313e-08 7.851563e-09
[27,] 1.0000000 2.574516e-08 1.287258e-08
[28,] 1.0000000 5.661408e-08 2.830704e-08
[29,] 0.9999999 1.367712e-07 6.838561e-08
[30,] 0.9999998 3.365318e-07 1.682659e-07
[31,] 0.9999996 7.131467e-07 3.565733e-07
[32,] 0.9999995 1.071875e-06 5.359374e-07
[33,] 0.9999993 1.392547e-06 6.962735e-07
[34,] 0.9999985 3.099227e-06 1.549614e-06
[35,] 0.9999985 3.097706e-06 1.548853e-06
[36,] 0.9999997 6.078627e-07 3.039313e-07
[37,] 0.9999994 1.115966e-06 5.579832e-07
[38,] 0.9999996 8.692442e-07 4.346221e-07
[39,] 0.9999989 2.142837e-06 1.071418e-06
[40,] 0.9999985 3.023832e-06 1.511916e-06
[41,] 0.9999970 6.068599e-06 3.034299e-06
[42,] 0.9999952 9.634614e-06 4.817307e-06
[43,] 0.9999902 1.951971e-05 9.759856e-06
[44,] 0.9999822 3.557964e-05 1.778982e-05
[45,] 0.9999962 7.666766e-06 3.833383e-06
[46,] 0.9999928 1.443456e-05 7.217282e-06
[47,] 0.9999822 3.559889e-05 1.779944e-05
[48,] 0.9999606 7.878379e-05 3.939189e-05
[49,] 0.9999114 1.772563e-04 8.862816e-05
[50,] 0.9998755 2.490143e-04 1.245071e-04
[51,] 0.9997199 5.601370e-04 2.800685e-04
[52,] 0.9999991 1.727438e-06 8.637192e-07
[53,] 0.9999980 4.011797e-06 2.005899e-06
[54,] 0.9999977 4.603965e-06 2.301983e-06
[55,] 0.9999974 5.126337e-06 2.563168e-06
[56,] 0.9999994 1.274692e-06 6.373460e-07
[57,] 0.9999980 3.916903e-06 1.958451e-06
[58,] 0.9999945 1.104876e-05 5.524378e-06
[59,] 0.9999787 4.264075e-05 2.132038e-05
[60,] 0.9999657 6.860981e-05 3.430491e-05
[61,] 0.9999152 1.696919e-04 8.484593e-05
[62,] 0.9997663 4.674621e-04 2.337310e-04
[63,] 0.9993198 1.360377e-03 6.801884e-04
[64,] 0.9978045 4.391028e-03 2.195514e-03
[65,] 0.9932524 1.349522e-02 6.747611e-03
[66,] 0.9779860 4.402801e-02 2.201401e-02
[67,] 0.9694519 6.109620e-02 3.054810e-02
> postscript(file="/var/wessaorg/rcomp/tmp/1t6j41322140929.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/2ghuz1322140929.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/3pdc91322140929.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/4h9c21322140929.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/5ioy41322140929.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
-397.258820 -1257.526126 1633.808108 506.465706 -17.284726 -1902.094076
7 8 9 10 11 12
1989.531380 -737.465631 280.766866 1997.464065 351.677672 -82.482748
13 14 15 16 17 18
296.507563 390.238738 -665.131883 -626.661151 -569.630669 -867.926239
19 20 21 22 23 24
-47.805500 -135.091166 -1149.224322 -662.948973 971.347085 -486.066213
25 26 27 28 29 30
-397.079678 -589.074774 -345.169322 848.840737 -422.847984 181.987349
31 32 33 34 35 36
-224.768865 -393.636529 -317.019019 -118.555139 -1.210046 -68.159657
37 38 39 40 41 42
-30.786947 -254.937673 -372.466622 -38.118940 -87.310862 650.056367
43 44 45 46 47 48
99.604314 525.007500 22.196401 -470.099223 46.161568 -20.536861
49 50 51 52 53 54
364.064050 108.924930 1068.501034 93.208806 324.645542 151.644181
55 56 57 58 59 60
254.713044 230.928854 471.560765 1225.532109 281.588913 189.561068
61 62 63 64 65 66
871.370802 268.293007 324.258876 621.911213 618.462102 -146.398066
67 68 69 70 71 72
431.868549 956.243665 -544.375324 -758.570098 -416.608316 -608.517902
73 74 75 76 77 78
-339.608152 -421.518365 -722.106485 -561.300927 -458.893335 -491.178165
79 80
-42.245920 -381.245487
> postscript(file="/var/wessaorg/rcomp/tmp/6x0dj1322140929.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 -397.258820 NA
1 -1257.526126 -397.258820
2 1633.808108 -1257.526126
3 506.465706 1633.808108
4 -17.284726 506.465706
5 -1902.094076 -17.284726
6 1989.531380 -1902.094076
7 -737.465631 1989.531380
8 280.766866 -737.465631
9 1997.464065 280.766866
10 351.677672 1997.464065
11 -82.482748 351.677672
12 296.507563 -82.482748
13 390.238738 296.507563
14 -665.131883 390.238738
15 -626.661151 -665.131883
16 -569.630669 -626.661151
17 -867.926239 -569.630669
18 -47.805500 -867.926239
19 -135.091166 -47.805500
20 -1149.224322 -135.091166
21 -662.948973 -1149.224322
22 971.347085 -662.948973
23 -486.066213 971.347085
24 -397.079678 -486.066213
25 -589.074774 -397.079678
26 -345.169322 -589.074774
27 848.840737 -345.169322
28 -422.847984 848.840737
29 181.987349 -422.847984
30 -224.768865 181.987349
31 -393.636529 -224.768865
32 -317.019019 -393.636529
33 -118.555139 -317.019019
34 -1.210046 -118.555139
35 -68.159657 -1.210046
36 -30.786947 -68.159657
37 -254.937673 -30.786947
38 -372.466622 -254.937673
39 -38.118940 -372.466622
40 -87.310862 -38.118940
41 650.056367 -87.310862
42 99.604314 650.056367
43 525.007500 99.604314
44 22.196401 525.007500
45 -470.099223 22.196401
46 46.161568 -470.099223
47 -20.536861 46.161568
48 364.064050 -20.536861
49 108.924930 364.064050
50 1068.501034 108.924930
51 93.208806 1068.501034
52 324.645542 93.208806
53 151.644181 324.645542
54 254.713044 151.644181
55 230.928854 254.713044
56 471.560765 230.928854
57 1225.532109 471.560765
58 281.588913 1225.532109
59 189.561068 281.588913
60 871.370802 189.561068
61 268.293007 871.370802
62 324.258876 268.293007
63 621.911213 324.258876
64 618.462102 621.911213
65 -146.398066 618.462102
66 431.868549 -146.398066
67 956.243665 431.868549
68 -544.375324 956.243665
69 -758.570098 -544.375324
70 -416.608316 -758.570098
71 -608.517902 -416.608316
72 -339.608152 -608.517902
73 -421.518365 -339.608152
74 -722.106485 -421.518365
75 -561.300927 -722.106485
76 -458.893335 -561.300927
77 -491.178165 -458.893335
78 -42.245920 -491.178165
79 -381.245487 -42.245920
80 NA -381.245487
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1257.526126 -397.258820
[2,] 1633.808108 -1257.526126
[3,] 506.465706 1633.808108
[4,] -17.284726 506.465706
[5,] -1902.094076 -17.284726
[6,] 1989.531380 -1902.094076
[7,] -737.465631 1989.531380
[8,] 280.766866 -737.465631
[9,] 1997.464065 280.766866
[10,] 351.677672 1997.464065
[11,] -82.482748 351.677672
[12,] 296.507563 -82.482748
[13,] 390.238738 296.507563
[14,] -665.131883 390.238738
[15,] -626.661151 -665.131883
[16,] -569.630669 -626.661151
[17,] -867.926239 -569.630669
[18,] -47.805500 -867.926239
[19,] -135.091166 -47.805500
[20,] -1149.224322 -135.091166
[21,] -662.948973 -1149.224322
[22,] 971.347085 -662.948973
[23,] -486.066213 971.347085
[24,] -397.079678 -486.066213
[25,] -589.074774 -397.079678
[26,] -345.169322 -589.074774
[27,] 848.840737 -345.169322
[28,] -422.847984 848.840737
[29,] 181.987349 -422.847984
[30,] -224.768865 181.987349
[31,] -393.636529 -224.768865
[32,] -317.019019 -393.636529
[33,] -118.555139 -317.019019
[34,] -1.210046 -118.555139
[35,] -68.159657 -1.210046
[36,] -30.786947 -68.159657
[37,] -254.937673 -30.786947
[38,] -372.466622 -254.937673
[39,] -38.118940 -372.466622
[40,] -87.310862 -38.118940
[41,] 650.056367 -87.310862
[42,] 99.604314 650.056367
[43,] 525.007500 99.604314
[44,] 22.196401 525.007500
[45,] -470.099223 22.196401
[46,] 46.161568 -470.099223
[47,] -20.536861 46.161568
[48,] 364.064050 -20.536861
[49,] 108.924930 364.064050
[50,] 1068.501034 108.924930
[51,] 93.208806 1068.501034
[52,] 324.645542 93.208806
[53,] 151.644181 324.645542
[54,] 254.713044 151.644181
[55,] 230.928854 254.713044
[56,] 471.560765 230.928854
[57,] 1225.532109 471.560765
[58,] 281.588913 1225.532109
[59,] 189.561068 281.588913
[60,] 871.370802 189.561068
[61,] 268.293007 871.370802
[62,] 324.258876 268.293007
[63,] 621.911213 324.258876
[64,] 618.462102 621.911213
[65,] -146.398066 618.462102
[66,] 431.868549 -146.398066
[67,] 956.243665 431.868549
[68,] -544.375324 956.243665
[69,] -758.570098 -544.375324
[70,] -416.608316 -758.570098
[71,] -608.517902 -416.608316
[72,] -339.608152 -608.517902
[73,] -421.518365 -339.608152
[74,] -722.106485 -421.518365
[75,] -561.300927 -722.106485
[76,] -458.893335 -561.300927
[77,] -491.178165 -458.893335
[78,] -42.245920 -491.178165
[79,] -381.245487 -42.245920
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1257.526126 -397.258820
2 1633.808108 -1257.526126
3 506.465706 1633.808108
4 -17.284726 506.465706
5 -1902.094076 -17.284726
6 1989.531380 -1902.094076
7 -737.465631 1989.531380
8 280.766866 -737.465631
9 1997.464065 280.766866
10 351.677672 1997.464065
11 -82.482748 351.677672
12 296.507563 -82.482748
13 390.238738 296.507563
14 -665.131883 390.238738
15 -626.661151 -665.131883
16 -569.630669 -626.661151
17 -867.926239 -569.630669
18 -47.805500 -867.926239
19 -135.091166 -47.805500
20 -1149.224322 -135.091166
21 -662.948973 -1149.224322
22 971.347085 -662.948973
23 -486.066213 971.347085
24 -397.079678 -486.066213
25 -589.074774 -397.079678
26 -345.169322 -589.074774
27 848.840737 -345.169322
28 -422.847984 848.840737
29 181.987349 -422.847984
30 -224.768865 181.987349
31 -393.636529 -224.768865
32 -317.019019 -393.636529
33 -118.555139 -317.019019
34 -1.210046 -118.555139
35 -68.159657 -1.210046
36 -30.786947 -68.159657
37 -254.937673 -30.786947
38 -372.466622 -254.937673
39 -38.118940 -372.466622
40 -87.310862 -38.118940
41 650.056367 -87.310862
42 99.604314 650.056367
43 525.007500 99.604314
44 22.196401 525.007500
45 -470.099223 22.196401
46 46.161568 -470.099223
47 -20.536861 46.161568
48 364.064050 -20.536861
49 108.924930 364.064050
50 1068.501034 108.924930
51 93.208806 1068.501034
52 324.645542 93.208806
53 151.644181 324.645542
54 254.713044 151.644181
55 230.928854 254.713044
56 471.560765 230.928854
57 1225.532109 471.560765
58 281.588913 1225.532109
59 189.561068 281.588913
60 871.370802 189.561068
61 268.293007 871.370802
62 324.258876 268.293007
63 621.911213 324.258876
64 618.462102 621.911213
65 -146.398066 618.462102
66 431.868549 -146.398066
67 956.243665 431.868549
68 -544.375324 956.243665
69 -758.570098 -544.375324
70 -416.608316 -758.570098
71 -608.517902 -416.608316
72 -339.608152 -608.517902
73 -421.518365 -339.608152
74 -722.106485 -421.518365
75 -561.300927 -722.106485
76 -458.893335 -561.300927
77 -491.178165 -458.893335
78 -42.245920 -491.178165
79 -381.245487 -42.245920
> 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/7lmrx1322140929.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/8mo6g1322140929.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/9upka1322140929.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/107eaw1322140929.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/11sm591322140929.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/12wice1322140929.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/136aj41322140929.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/14bwow1322140929.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/15ifu51322140929.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/1630ly1322140929.tab")
+ }
>
> try(system("convert tmp/1t6j41322140929.ps tmp/1t6j41322140929.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ghuz1322140929.ps tmp/2ghuz1322140929.png",intern=TRUE))
character(0)
> try(system("convert tmp/3pdc91322140929.ps tmp/3pdc91322140929.png",intern=TRUE))
character(0)
> try(system("convert tmp/4h9c21322140929.ps tmp/4h9c21322140929.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ioy41322140929.ps tmp/5ioy41322140929.png",intern=TRUE))
character(0)
> try(system("convert tmp/6x0dj1322140929.ps tmp/6x0dj1322140929.png",intern=TRUE))
character(0)
> try(system("convert tmp/7lmrx1322140929.ps tmp/7lmrx1322140929.png",intern=TRUE))
character(0)
> try(system("convert tmp/8mo6g1322140929.ps tmp/8mo6g1322140929.png",intern=TRUE))
character(0)
> try(system("convert tmp/9upka1322140929.ps tmp/9upka1322140929.png",intern=TRUE))
character(0)
> try(system("convert tmp/107eaw1322140929.ps tmp/107eaw1322140929.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.663 0.518 4.236