R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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
+ ,2346
+ ,864
+ ,1068
+ ,2251
+ ,497
+ ,399
+ ,2230
+ ,449
+ ,547
+ ,2225
+ ,919
+ ,668
+ ,2220
+ ,536
+ ,451
+ ,2205
+ ,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
+ ,1684
+ ,804
+ ,625
+ ,1549
+ ,456
+ ,245
+ ,1533
+ ,267
+ ,235
+ ,1528
+ ,338
+ ,250
+ ,1489
+ ,NA
+ ,NA
+ ,1456
+ ,635
+ ,303
+ ,1393
+ ,402
+ ,250
+ ,1370
+ ,462
+ ,403
+ ,1357
+ ,525
+ ,441
+ ,1292
+ ,1069
+ ,1507
+ ,1278
+ ,487
+ ,388
+ ,1256
+ ,756
+ ,530
+ ,1219
+ ,360
+ ,107
+ ,1217
+ ,481
+ ,868
+ ,1187
+ ,NA
+ ,NA
+ ,1187
+ ,480
+ ,753
+ ,1182
+ ,68
+ ,62
+ ,1157
+ ,262
+ ,346
+ ,1156
+ ,520
+ ,478
+ ,1155
+ ,298
+ ,292)
+ ,dim=c(3
+ ,100)
+ ,dimnames=list(c('weekdag'
+ ,'zaterdag'
+ ,'zondag')
+ ,1:100))
> y <- array(NA,dim=c(3,100),dimnames=list(c('weekdag','zaterdag','zondag'),1:100))
> 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'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, as.Date.numeric
> 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 1068 2346 864
62 399 2251 497
63 547 2230 449
64 668 2225 919
65 451 2220 536
66 724 2205 673
67 853 2193 837
68 434 2116 534
69 730 2102 845
70 612 2099 626
71 558 2096 871
72 859 2064 740
73 311 2036 391
74 318 1920 435
75 312 1813 424
76 343 1776 338
77 710 1752 744
78 273 1738 368
79 259 1729 393
80 1274 1685 938
81 625 1684 804
82 245 1549 456
83 235 1533 267
84 250 1528 338
85 NA 1489 NA
86 303 1456 635
87 250 1393 402
88 403 1370 462
89 441 1357 525
90 1507 1292 1069
91 388 1278 487
92 530 1256 756
93 107 1219 360
94 868 1217 481
95 NA 1187 NA
96 753 1187 480
97 62 1182 68
98 346 1157 262
99 478 1156 520
100 292 1155 298
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) weekdag zaterdag
136.51725 -0.01725 0.93285
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1682.03 -190.59 -117.64 87.82 2457.70
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 136.51725 66.74894 2.045 0.0436 *
weekdag -0.01725 0.01772 -0.973 0.3328
zaterdag 0.93285 0.04323 21.578 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 556.5 on 95 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.9814, Adjusted R-squared: 0.981
F-statistic: 2504 on 2 and 95 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,] 1.0000000 4.188328e-09 2.094164e-09
[2,] 1.0000000 8.240749e-16 4.120375e-16
[3,] 1.0000000 6.932995e-18 3.466498e-18
[4,] 1.0000000 4.090118e-17 2.045059e-17
[5,] 1.0000000 1.448017e-19 7.240086e-20
[6,] 1.0000000 3.734062e-19 1.867031e-19
[7,] 1.0000000 7.823405e-19 3.911702e-19
[8,] 1.0000000 3.089296e-18 1.544648e-18
[9,] 1.0000000 1.417801e-17 7.089005e-18
[10,] 1.0000000 6.869389e-18 3.434695e-18
[11,] 1.0000000 7.371597e-18 3.685799e-18
[12,] 1.0000000 1.206371e-17 6.031853e-18
[13,] 1.0000000 4.948583e-18 2.474291e-18
[14,] 1.0000000 6.207916e-18 3.103958e-18
[15,] 1.0000000 2.324176e-17 1.162088e-17
[16,] 1.0000000 1.702917e-19 8.514587e-20
[17,] 1.0000000 2.842954e-19 1.421477e-19
[18,] 1.0000000 2.673118e-20 1.336559e-20
[19,] 1.0000000 8.019474e-20 4.009737e-20
[20,] 1.0000000 2.786151e-19 1.393076e-19
[21,] 1.0000000 5.918065e-19 2.959032e-19
[22,] 1.0000000 1.962418e-18 9.812088e-19
[23,] 1.0000000 8.806314e-19 4.403157e-19
[24,] 1.0000000 1.059563e-18 5.297816e-19
[25,] 1.0000000 4.210892e-18 2.105446e-18
[26,] 1.0000000 1.487302e-17 7.436510e-18
[27,] 1.0000000 4.727350e-17 2.363675e-17
[28,] 1.0000000 1.340015e-16 6.700077e-17
[29,] 1.0000000 3.374921e-16 1.687461e-16
[30,] 1.0000000 1.181346e-15 5.906728e-16
[31,] 1.0000000 3.939284e-15 1.969642e-15
[32,] 1.0000000 1.238614e-14 6.193070e-15
[33,] 1.0000000 3.030279e-14 1.515140e-14
[34,] 1.0000000 2.681275e-14 1.340638e-14
[35,] 1.0000000 6.982054e-14 3.491027e-14
[36,] 1.0000000 2.087322e-13 1.043661e-13
[37,] 1.0000000 2.621888e-13 1.310944e-13
[38,] 1.0000000 8.096958e-13 4.048479e-13
[39,] 1.0000000 9.514323e-13 4.757161e-13
[40,] 1.0000000 2.583403e-12 1.291701e-12
[41,] 1.0000000 1.039387e-13 5.196937e-14
[42,] 1.0000000 2.679356e-13 1.339678e-13
[43,] 1.0000000 4.538452e-13 2.269226e-13
[44,] 1.0000000 1.465150e-12 7.325751e-13
[45,] 1.0000000 3.622569e-12 1.811284e-12
[46,] 1.0000000 2.938253e-14 1.469127e-14
[47,] 1.0000000 9.028102e-14 4.514051e-14
[48,] 1.0000000 3.063607e-13 1.531803e-13
[49,] 1.0000000 8.122934e-13 4.061467e-13
[50,] 1.0000000 2.672150e-12 1.336075e-12
[51,] 1.0000000 8.610833e-12 4.305417e-12
[52,] 1.0000000 2.347280e-11 1.173640e-11
[53,] 1.0000000 1.295625e-14 6.478123e-15
[54,] 1.0000000 4.875339e-14 2.437670e-14
[55,] 1.0000000 1.826111e-13 9.130553e-14
[56,] 1.0000000 2.869330e-13 1.434665e-13
[57,] 1.0000000 1.101783e-12 5.508914e-13
[58,] 1.0000000 1.914259e-12 9.571296e-13
[59,] 1.0000000 3.688273e-12 1.844136e-12
[60,] 1.0000000 1.388449e-11 6.942246e-12
[61,] 1.0000000 3.908031e-11 1.954015e-11
[62,] 1.0000000 1.334745e-10 6.673723e-11
[63,] 1.0000000 4.758865e-10 2.379432e-10
[64,] 1.0000000 1.622019e-09 8.110093e-10
[65,] 1.0000000 5.041850e-09 2.520925e-09
[66,] 1.0000000 5.871163e-09 2.935582e-09
[67,] 1.0000000 1.468412e-08 7.342058e-09
[68,] 1.0000000 4.574582e-08 2.287291e-08
[69,] 0.9999999 1.524614e-07 7.623068e-08
[70,] 0.9999998 4.987989e-07 2.493994e-07
[71,] 0.9999994 1.127577e-06 5.637883e-07
[72,] 0.9999982 3.602352e-06 1.801176e-06
[73,] 0.9999949 1.013476e-05 5.067382e-06
[74,] 0.9999855 2.899327e-05 1.449663e-05
[75,] 0.9999913 1.741487e-05 8.707435e-06
[76,] 0.9999726 5.472693e-05 2.736347e-05
[77,] 0.9999190 1.620801e-04 8.104005e-05
[78,] 0.9998441 3.118383e-04 1.559191e-04
[79,] 0.9997435 5.130837e-04 2.565418e-04
[80,] 0.9994420 1.115931e-03 5.579657e-04
[81,] 0.9983701 3.259865e-03 1.629933e-03
[82,] 0.9956041 8.791780e-03 4.395890e-03
[83,] 0.9884534 2.309327e-02 1.154663e-02
[84,] 0.9867528 2.649443e-02 1.324721e-02
[85,] 0.9647006 7.059884e-02 3.529942e-02
[86,] 0.9489018 1.021964e-01 5.109820e-02
[87,] 0.9935372 1.292561e-02 6.462807e-03
[88,] 0.9581294 8.374123e-02 4.187061e-02
[89,] 0.9144659 1.710682e-01 8.553411e-02
> postscript(file="/var/wessaorg/rcomp/tmp/1bgx41352118421.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)
Warning message:
In x[, 1] - mysum$resid :
longer object length is not a multiple of shorter object length
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/23j0p1352118421.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/3enyu1352118421.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/4kha11352118421.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/5pb4v1352118421.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 = 98
Frequency = 1
1 2 3 4 5 6
-929.692356 -1427.045158 1681.326765 529.138385 -726.200546 -1682.028127
7 8 9 10 11 12
2457.701841 -427.545264 630.393101 2047.079926 552.472484 327.222014
13 14 15 16 17 18
650.318054 670.421872 -153.603269 -126.496926 -82.794141 -454.647998
19 20 21 22 23 24
382.788157 194.530968 -815.909628 -302.351624 1154.024980 -157.838736
25 26 27 28 29 30
-81.428671 -294.054091 -42.589141 890.600083 -283.724852 237.868380
31 32 33 34 35 36
-35.660743 -175.026222 -181.655360 -86.423272 62.905991 -18.316725
37 38 39 40 41 42
32.347133 -220.971002 -381.541282 -78.673370 -95.399019 516.601039
43 44 45 46 47 48
44.189909 394.513197 -136.336912 -664.639788 -121.559509 -210.916504
49 50 51 52 53 54
130.095096 -117.716413 823.269519 -168.229819 32.904496 -167.480623
55 56 57 58 59 60
-80.169348 -96.360469 94.506305 812.183329 -124.986188 -235.554179
61 62 63 64 65 66
165.959301 -162.322069 30.092693 -287.434519 -147.238001 -2.297588
67 68 69 70 71 72
-26.492468 -164.165991 -158.524778 -72.281681 -354.882442 67.769413
73 74 75 76 77 78
-155.147759 -191.193959 -188.778012 -78.190784 -90.343090 -176.831768
79 80 81 82 83 84
-214.308320 291.527846 -232.487080 -290.182543 -124.149253 -175.468061
86 87 88 89 90 91
-400.767235 -237.499019 -140.866890 -161.860850 395.545980 -180.774951
92 93 94 96 97 98
-290.091882 -344.320180 303.770096 189.185536 -117.565204 -14.969891
99 100
-123.663249 -102.587098
> postscript(file="/var/wessaorg/rcomp/tmp/6dyur1352118421.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 = 98
Frequency = 1
lag(myerror, k = 1) myerror
0 -929.692356 NA
1 -1427.045158 -929.692356
2 1681.326765 -1427.045158
3 529.138385 1681.326765
4 -726.200546 529.138385
5 -1682.028127 -726.200546
6 2457.701841 -1682.028127
7 -427.545264 2457.701841
8 630.393101 -427.545264
9 2047.079926 630.393101
10 552.472484 2047.079926
11 327.222014 552.472484
12 650.318054 327.222014
13 670.421872 650.318054
14 -153.603269 670.421872
15 -126.496926 -153.603269
16 -82.794141 -126.496926
17 -454.647998 -82.794141
18 382.788157 -454.647998
19 194.530968 382.788157
20 -815.909628 194.530968
21 -302.351624 -815.909628
22 1154.024980 -302.351624
23 -157.838736 1154.024980
24 -81.428671 -157.838736
25 -294.054091 -81.428671
26 -42.589141 -294.054091
27 890.600083 -42.589141
28 -283.724852 890.600083
29 237.868380 -283.724852
30 -35.660743 237.868380
31 -175.026222 -35.660743
32 -181.655360 -175.026222
33 -86.423272 -181.655360
34 62.905991 -86.423272
35 -18.316725 62.905991
36 32.347133 -18.316725
37 -220.971002 32.347133
38 -381.541282 -220.971002
39 -78.673370 -381.541282
40 -95.399019 -78.673370
41 516.601039 -95.399019
42 44.189909 516.601039
43 394.513197 44.189909
44 -136.336912 394.513197
45 -664.639788 -136.336912
46 -121.559509 -664.639788
47 -210.916504 -121.559509
48 130.095096 -210.916504
49 -117.716413 130.095096
50 823.269519 -117.716413
51 -168.229819 823.269519
52 32.904496 -168.229819
53 -167.480623 32.904496
54 -80.169348 -167.480623
55 -96.360469 -80.169348
56 94.506305 -96.360469
57 812.183329 94.506305
58 -124.986188 812.183329
59 -235.554179 -124.986188
60 165.959301 -235.554179
61 -162.322069 165.959301
62 30.092693 -162.322069
63 -287.434519 30.092693
64 -147.238001 -287.434519
65 -2.297588 -147.238001
66 -26.492468 -2.297588
67 -164.165991 -26.492468
68 -158.524778 -164.165991
69 -72.281681 -158.524778
70 -354.882442 -72.281681
71 67.769413 -354.882442
72 -155.147759 67.769413
73 -191.193959 -155.147759
74 -188.778012 -191.193959
75 -78.190784 -188.778012
76 -90.343090 -78.190784
77 -176.831768 -90.343090
78 -214.308320 -176.831768
79 291.527846 -214.308320
80 -232.487080 291.527846
81 -290.182543 -232.487080
82 -124.149253 -290.182543
83 -175.468061 -124.149253
84 -400.767235 -175.468061
85 -237.499019 -400.767235
86 -140.866890 -237.499019
87 -161.860850 -140.866890
88 395.545980 -161.860850
89 -180.774951 395.545980
90 -290.091882 -180.774951
91 -344.320180 -290.091882
92 303.770096 -344.320180
93 189.185536 303.770096
94 -117.565204 189.185536
95 -14.969891 -117.565204
96 -123.663249 -14.969891
97 -102.587098 -123.663249
98 NA -102.587098
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1427.045158 -929.692356
[2,] 1681.326765 -1427.045158
[3,] 529.138385 1681.326765
[4,] -726.200546 529.138385
[5,] -1682.028127 -726.200546
[6,] 2457.701841 -1682.028127
[7,] -427.545264 2457.701841
[8,] 630.393101 -427.545264
[9,] 2047.079926 630.393101
[10,] 552.472484 2047.079926
[11,] 327.222014 552.472484
[12,] 650.318054 327.222014
[13,] 670.421872 650.318054
[14,] -153.603269 670.421872
[15,] -126.496926 -153.603269
[16,] -82.794141 -126.496926
[17,] -454.647998 -82.794141
[18,] 382.788157 -454.647998
[19,] 194.530968 382.788157
[20,] -815.909628 194.530968
[21,] -302.351624 -815.909628
[22,] 1154.024980 -302.351624
[23,] -157.838736 1154.024980
[24,] -81.428671 -157.838736
[25,] -294.054091 -81.428671
[26,] -42.589141 -294.054091
[27,] 890.600083 -42.589141
[28,] -283.724852 890.600083
[29,] 237.868380 -283.724852
[30,] -35.660743 237.868380
[31,] -175.026222 -35.660743
[32,] -181.655360 -175.026222
[33,] -86.423272 -181.655360
[34,] 62.905991 -86.423272
[35,] -18.316725 62.905991
[36,] 32.347133 -18.316725
[37,] -220.971002 32.347133
[38,] -381.541282 -220.971002
[39,] -78.673370 -381.541282
[40,] -95.399019 -78.673370
[41,] 516.601039 -95.399019
[42,] 44.189909 516.601039
[43,] 394.513197 44.189909
[44,] -136.336912 394.513197
[45,] -664.639788 -136.336912
[46,] -121.559509 -664.639788
[47,] -210.916504 -121.559509
[48,] 130.095096 -210.916504
[49,] -117.716413 130.095096
[50,] 823.269519 -117.716413
[51,] -168.229819 823.269519
[52,] 32.904496 -168.229819
[53,] -167.480623 32.904496
[54,] -80.169348 -167.480623
[55,] -96.360469 -80.169348
[56,] 94.506305 -96.360469
[57,] 812.183329 94.506305
[58,] -124.986188 812.183329
[59,] -235.554179 -124.986188
[60,] 165.959301 -235.554179
[61,] -162.322069 165.959301
[62,] 30.092693 -162.322069
[63,] -287.434519 30.092693
[64,] -147.238001 -287.434519
[65,] -2.297588 -147.238001
[66,] -26.492468 -2.297588
[67,] -164.165991 -26.492468
[68,] -158.524778 -164.165991
[69,] -72.281681 -158.524778
[70,] -354.882442 -72.281681
[71,] 67.769413 -354.882442
[72,] -155.147759 67.769413
[73,] -191.193959 -155.147759
[74,] -188.778012 -191.193959
[75,] -78.190784 -188.778012
[76,] -90.343090 -78.190784
[77,] -176.831768 -90.343090
[78,] -214.308320 -176.831768
[79,] 291.527846 -214.308320
[80,] -232.487080 291.527846
[81,] -290.182543 -232.487080
[82,] -124.149253 -290.182543
[83,] -175.468061 -124.149253
[84,] -400.767235 -175.468061
[85,] -237.499019 -400.767235
[86,] -140.866890 -237.499019
[87,] -161.860850 -140.866890
[88,] 395.545980 -161.860850
[89,] -180.774951 395.545980
[90,] -290.091882 -180.774951
[91,] -344.320180 -290.091882
[92,] 303.770096 -344.320180
[93,] 189.185536 303.770096
[94,] -117.565204 189.185536
[95,] -14.969891 -117.565204
[96,] -123.663249 -14.969891
[97,] -102.587098 -123.663249
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1427.045158 -929.692356
2 1681.326765 -1427.045158
3 529.138385 1681.326765
4 -726.200546 529.138385
5 -1682.028127 -726.200546
6 2457.701841 -1682.028127
7 -427.545264 2457.701841
8 630.393101 -427.545264
9 2047.079926 630.393101
10 552.472484 2047.079926
11 327.222014 552.472484
12 650.318054 327.222014
13 670.421872 650.318054
14 -153.603269 670.421872
15 -126.496926 -153.603269
16 -82.794141 -126.496926
17 -454.647998 -82.794141
18 382.788157 -454.647998
19 194.530968 382.788157
20 -815.909628 194.530968
21 -302.351624 -815.909628
22 1154.024980 -302.351624
23 -157.838736 1154.024980
24 -81.428671 -157.838736
25 -294.054091 -81.428671
26 -42.589141 -294.054091
27 890.600083 -42.589141
28 -283.724852 890.600083
29 237.868380 -283.724852
30 -35.660743 237.868380
31 -175.026222 -35.660743
32 -181.655360 -175.026222
33 -86.423272 -181.655360
34 62.905991 -86.423272
35 -18.316725 62.905991
36 32.347133 -18.316725
37 -220.971002 32.347133
38 -381.541282 -220.971002
39 -78.673370 -381.541282
40 -95.399019 -78.673370
41 516.601039 -95.399019
42 44.189909 516.601039
43 394.513197 44.189909
44 -136.336912 394.513197
45 -664.639788 -136.336912
46 -121.559509 -664.639788
47 -210.916504 -121.559509
48 130.095096 -210.916504
49 -117.716413 130.095096
50 823.269519 -117.716413
51 -168.229819 823.269519
52 32.904496 -168.229819
53 -167.480623 32.904496
54 -80.169348 -167.480623
55 -96.360469 -80.169348
56 94.506305 -96.360469
57 812.183329 94.506305
58 -124.986188 812.183329
59 -235.554179 -124.986188
60 165.959301 -235.554179
61 -162.322069 165.959301
62 30.092693 -162.322069
63 -287.434519 30.092693
64 -147.238001 -287.434519
65 -2.297588 -147.238001
66 -26.492468 -2.297588
67 -164.165991 -26.492468
68 -158.524778 -164.165991
69 -72.281681 -158.524778
70 -354.882442 -72.281681
71 67.769413 -354.882442
72 -155.147759 67.769413
73 -191.193959 -155.147759
74 -188.778012 -191.193959
75 -78.190784 -188.778012
76 -90.343090 -78.190784
77 -176.831768 -90.343090
78 -214.308320 -176.831768
79 291.527846 -214.308320
80 -232.487080 291.527846
81 -290.182543 -232.487080
82 -124.149253 -290.182543
83 -175.468061 -124.149253
84 -400.767235 -175.468061
85 -237.499019 -400.767235
86 -140.866890 -237.499019
87 -161.860850 -140.866890
88 395.545980 -161.860850
89 -180.774951 395.545980
90 -290.091882 -180.774951
91 -344.320180 -290.091882
92 303.770096 -344.320180
93 189.185536 303.770096
94 -117.565204 189.185536
95 -14.969891 -117.565204
96 -123.663249 -14.969891
97 -102.587098 -123.663249
> 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/74r4j1352118421.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/870fy1352118421.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/98xph1352118421.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/10o1aj1352118421.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/11u9b51352118421.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/12n2hp1352118421.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/13587e1352118421.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/14jq251352118421.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/150n1u1352118421.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/16qica1352118421.tab")
+ }
>
> try(system("convert tmp/1bgx41352118421.ps tmp/1bgx41352118421.png",intern=TRUE))
character(0)
> try(system("convert tmp/23j0p1352118421.ps tmp/23j0p1352118421.png",intern=TRUE))
character(0)
> try(system("convert tmp/3enyu1352118421.ps tmp/3enyu1352118421.png",intern=TRUE))
character(0)
> try(system("convert tmp/4kha11352118421.ps tmp/4kha11352118421.png",intern=TRUE))
character(0)
> try(system("convert tmp/5pb4v1352118421.ps tmp/5pb4v1352118421.png",intern=TRUE))
character(0)
> try(system("convert tmp/6dyur1352118421.ps tmp/6dyur1352118421.png",intern=TRUE))
character(0)
> try(system("convert tmp/74r4j1352118421.ps tmp/74r4j1352118421.png",intern=TRUE))
character(0)
> try(system("convert tmp/870fy1352118421.ps tmp/870fy1352118421.png",intern=TRUE))
character(0)
> try(system("convert tmp/98xph1352118421.ps tmp/98xph1352118421.png",intern=TRUE))
character(0)
> try(system("convert tmp/10o1aj1352118421.ps tmp/10o1aj1352118421.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
7.138 0.893 8.025