R version 2.12.0 (2010-10-15)
Copyright (C) 2010 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(921365
+ ,18919
+ ,48873
+ ,137852
+ ,1
+ ,987921
+ ,19147
+ ,52118
+ ,145224
+ ,2
+ ,1132614
+ ,21518
+ ,60530
+ ,163575
+ ,3
+ ,1332224
+ ,20941
+ ,55644
+ ,190761
+ ,4
+ ,1418133
+ ,22401
+ ,57121
+ ,196562
+ ,5
+ ,1411549
+ ,22181
+ ,55697
+ ,204493
+ ,6
+ ,1695920
+ ,22494
+ ,56483
+ ,259479
+ ,7
+ ,1636173
+ ,21479
+ ,51541
+ ,259479
+ ,8
+ ,1539653
+ ,22322
+ ,56328
+ ,223164
+ ,9
+ ,1395314
+ ,21829
+ ,54349
+ ,194886
+ ,10
+ ,1127575
+ ,20370
+ ,59885
+ ,160407
+ ,11
+ ,1036076
+ ,18467
+ ,55806
+ ,151747
+ ,12
+ ,989236
+ ,18780
+ ,54559
+ ,152448
+ ,1
+ ,1008380
+ ,18815
+ ,55590
+ ,148388
+ ,2
+ ,1207763
+ ,20881
+ ,63442
+ ,168510
+ ,3
+ ,1368839
+ ,21443
+ ,61258
+ ,188041
+ ,4
+ ,1469798
+ ,22333
+ ,55829
+ ,192020
+ ,5
+ ,1498721
+ ,22944
+ ,58023
+ ,205250
+ ,6
+ ,1761769
+ ,22536
+ ,58887
+ ,261642
+ ,7
+ ,1653214
+ ,21658
+ ,51510
+ ,251614
+ ,8
+ ,1599104
+ ,23035
+ ,60006
+ ,222726
+ ,9
+ ,1421179
+ ,21969
+ ,60831
+ ,179039
+ ,10
+ ,1163995
+ ,20297
+ ,61559
+ ,151462
+ ,11
+ ,1037735
+ ,18564
+ ,61325
+ ,143653
+ ,12
+ ,1015407
+ ,18844
+ ,55222
+ ,143762
+ ,1
+ ,1039210
+ ,18762
+ ,56370
+ ,134580
+ ,2
+ ,1258049
+ ,21757
+ ,66063
+ ,165273
+ ,3
+ ,1469445
+ ,20501
+ ,60864
+ ,181016
+ ,4
+ ,1552346
+ ,23181
+ ,57596
+ ,189079
+ ,5
+ ,1549144
+ ,23015
+ ,57650
+ ,199266
+ ,6
+ ,1785895
+ ,22828
+ ,55324
+ ,248742
+ ,7
+ ,1662335
+ ,21597
+ ,54203
+ ,244139
+ ,8
+ ,1629440
+ ,23005
+ ,61155
+ ,219777
+ ,9
+ ,1467430
+ ,22243
+ ,63908
+ ,180679
+ ,10
+ ,1202209
+ ,20729
+ ,67466
+ ,156369
+ ,11
+ ,1076982
+ ,18310
+ ,63739
+ ,149176
+ ,12
+ ,1039367
+ ,19427
+ ,56602
+ ,147247
+ ,1
+ ,1063449
+ ,18849
+ ,57640
+ ,142026
+ ,2
+ ,1335135
+ ,21817
+ ,70025
+ ,174119
+ ,3
+ ,1491602
+ ,21101
+ ,61068
+ ,190271
+ ,4
+ ,1591972
+ ,23546
+ ,60467
+ ,202998
+ ,5
+ ,1641248
+ ,23456
+ ,65297
+ ,219097
+ ,6
+ ,1898849
+ ,23649
+ ,64505
+ ,266542
+ ,7
+ ,1798580
+ ,22432
+ ,62517
+ ,257522
+ ,8
+ ,1762444
+ ,23745
+ ,67403
+ ,226187
+ ,9
+ ,1622044
+ ,23874
+ ,70508
+ ,196827
+ ,10
+ ,1368955
+ ,22327
+ ,75601
+ ,174065
+ ,11
+ ,1262973
+ ,20143
+ ,72094
+ ,165891
+ ,12
+ ,1195650
+ ,21252
+ ,66527
+ ,153950
+ ,1
+ ,1269530
+ ,21094
+ ,69324
+ ,154796
+ ,2
+ ,1479279
+ ,21800
+ ,75423
+ ,179944
+ ,3
+ ,1607819
+ ,22480
+ ,57761
+ ,195820
+ ,4
+ ,1712466
+ ,23055
+ ,55801
+ ,203015
+ ,5
+ ,1721766
+ ,23352
+ ,52949
+ ,214055
+ ,6
+ ,1949843
+ ,23171
+ ,45719
+ ,256871
+ ,7
+ ,1821326
+ ,20691
+ ,46610
+ ,235046
+ ,8
+ ,1757802
+ ,23183
+ ,48713
+ ,214295
+ ,9
+ ,1590367
+ ,22412
+ ,50018
+ ,191605
+ ,10
+ ,1260647
+ ,18958
+ ,49123
+ ,159512
+ ,11
+ ,1149235
+ ,17347
+ ,43157
+ ,149715
+ ,12
+ ,1016367
+ ,17353
+ ,36613
+ ,131871
+ ,1
+ ,1027885
+ ,17153
+ ,38355
+ ,130864
+ ,2
+ ,1262159
+ ,20141
+ ,42107
+ ,154383
+ ,3
+ ,1520854
+ ,19699
+ ,36495
+ ,178030
+ ,4
+ ,1544144
+ ,20780
+ ,35589
+ ,183488
+ ,5
+ ,1564709
+ ,21101
+ ,36864
+ ,204119
+ ,6
+ ,1821776
+ ,20871
+ ,36068
+ ,237511
+ ,7
+ ,1741365
+ ,19574
+ ,25131
+ ,228871
+ ,8
+ ,1623386
+ ,21002
+ ,35198
+ ,196125
+ ,9
+ ,1498658
+ ,20105
+ ,38749
+ ,177142
+ ,10
+ ,1241822
+ ,17772
+ ,39385
+ ,151338
+ ,11
+ ,1136029
+ ,16117
+ ,38579
+ ,144732
+ ,12)
+ ,dim=c(5
+ ,72)
+ ,dimnames=list(c('passagiers'
+ ,'bewegingen'
+ ,'cargo'
+ ,'auto'
+ ,'maand')
+ ,1:72))
> y <- array(NA,dim=c(5,72),dimnames=list(c('passagiers','bewegingen','cargo','auto','maand'),1:72))
> 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 = '2'
> #'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
bewegingen passagiers cargo auto maand
1 18919 921365 48873 137852 1
2 19147 987921 52118 145224 2
3 21518 1132614 60530 163575 3
4 20941 1332224 55644 190761 4
5 22401 1418133 57121 196562 5
6 22181 1411549 55697 204493 6
7 22494 1695920 56483 259479 7
8 21479 1636173 51541 259479 8
9 22322 1539653 56328 223164 9
10 21829 1395314 54349 194886 10
11 20370 1127575 59885 160407 11
12 18467 1036076 55806 151747 12
13 18780 989236 54559 152448 1
14 18815 1008380 55590 148388 2
15 20881 1207763 63442 168510 3
16 21443 1368839 61258 188041 4
17 22333 1469798 55829 192020 5
18 22944 1498721 58023 205250 6
19 22536 1761769 58887 261642 7
20 21658 1653214 51510 251614 8
21 23035 1599104 60006 222726 9
22 21969 1421179 60831 179039 10
23 20297 1163995 61559 151462 11
24 18564 1037735 61325 143653 12
25 18844 1015407 55222 143762 1
26 18762 1039210 56370 134580 2
27 21757 1258049 66063 165273 3
28 20501 1469445 60864 181016 4
29 23181 1552346 57596 189079 5
30 23015 1549144 57650 199266 6
31 22828 1785895 55324 248742 7
32 21597 1662335 54203 244139 8
33 23005 1629440 61155 219777 9
34 22243 1467430 63908 180679 10
35 20729 1202209 67466 156369 11
36 18310 1076982 63739 149176 12
37 19427 1039367 56602 147247 1
38 18849 1063449 57640 142026 2
39 21817 1335135 70025 174119 3
40 21101 1491602 61068 190271 4
41 23546 1591972 60467 202998 5
42 23456 1641248 65297 219097 6
43 23649 1898849 64505 266542 7
44 22432 1798580 62517 257522 8
45 23745 1762444 67403 226187 9
46 23874 1622044 70508 196827 10
47 22327 1368955 75601 174065 11
48 20143 1262973 72094 165891 12
49 21252 1195650 66527 153950 1
50 21094 1269530 69324 154796 2
51 21800 1479279 75423 179944 3
52 22480 1607819 57761 195820 4
53 23055 1712466 55801 203015 5
54 23352 1721766 52949 214055 6
55 23171 1949843 45719 256871 7
56 20691 1821326 46610 235046 8
57 23183 1757802 48713 214295 9
58 22412 1590367 50018 191605 10
59 18958 1260647 49123 159512 11
60 17347 1149235 43157 149715 12
61 17353 1016367 36613 131871 1
62 17153 1027885 38355 130864 2
63 20141 1262159 42107 154383 3
64 19699 1520854 36495 178030 4
65 20780 1544144 35589 183488 5
66 21101 1564709 36864 204119 6
67 20871 1821776 36068 237511 7
68 19574 1741365 25131 228871 8
69 21002 1623386 35198 196125 9
70 20105 1498658 38749 177142 10
71 17772 1241822 39385 151338 11
72 16117 1136029 38579 144732 12
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) passagiers cargo auto maand
8.775e+03 6.111e-03 8.621e-02 -3.291e-03 -7.975e+01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1819.87 -610.65 -3.86 657.96 1381.58
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.775e+03 7.240e+02 12.120 < 2e-16 ***
passagiers 6.110e-03 9.032e-04 6.766 3.97e-09 ***
cargo 8.621e-02 8.944e-03 9.639 2.79e-14 ***
auto -3.291e-03 6.473e-03 -0.508 0.6128
maand -7.975e+01 2.792e+01 -2.856 0.0057 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 790.8 on 67 degrees of freedom
Multiple R-squared: 0.8295, Adjusted R-squared: 0.8193
F-statistic: 81.47 on 4 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.17040942 0.3408188 0.8295906
[2,] 0.11744853 0.2348971 0.8825515
[3,] 0.05970083 0.1194017 0.9402992
[4,] 0.06990929 0.1398186 0.9300907
[5,] 0.08514360 0.1702872 0.9148564
[6,] 0.07248104 0.1449621 0.9275190
[7,] 0.09771858 0.1954372 0.9022814
[8,] 0.10392990 0.2078598 0.8960701
[9,] 0.13257048 0.2651410 0.8674295
[10,] 0.12195586 0.2439117 0.8780441
[11,] 0.11320636 0.2264127 0.8867936
[12,] 0.11950236 0.2390047 0.8804976
[13,] 0.09065374 0.1813075 0.9093463
[14,] 0.07464264 0.1492853 0.9253574
[15,] 0.08136615 0.1627323 0.9186339
[16,] 0.07096288 0.1419258 0.9290371
[17,] 0.07904758 0.1580952 0.9209524
[18,] 0.09080894 0.1816179 0.9091911
[19,] 0.16577847 0.3315569 0.8342215
[20,] 0.13603385 0.2720677 0.8639661
[21,] 0.56003689 0.8799262 0.4399631
[22,] 0.54572413 0.9085517 0.4542759
[23,] 0.55534406 0.8893119 0.4446559
[24,] 0.52216562 0.9556688 0.4778344
[25,] 0.51008926 0.9798215 0.4899107
[26,] 0.49206926 0.9841385 0.5079307
[27,] 0.43563918 0.8712784 0.5643608
[28,] 0.40302971 0.8060594 0.5969703
[29,] 0.45521524 0.9104305 0.5447848
[30,] 0.43990014 0.8798003 0.5600999
[31,] 0.43000139 0.8600028 0.5699986
[32,] 0.36292976 0.7258595 0.6370702
[33,] 0.43985074 0.8797015 0.5601493
[34,] 0.48221895 0.9644379 0.5177811
[35,] 0.46563026 0.9312605 0.5343697
[36,] 0.42305643 0.8461129 0.5769436
[37,] 0.40558180 0.8111636 0.5944182
[38,] 0.33763515 0.6752703 0.6623649
[39,] 0.31206895 0.6241379 0.6879311
[40,] 0.31710051 0.6342010 0.6828995
[41,] 0.28271429 0.5654286 0.7172857
[42,] 0.28874243 0.5774849 0.7112576
[43,] 0.23335660 0.4667132 0.7666434
[44,] 0.37693154 0.7538631 0.6230685
[45,] 0.34506054 0.6901211 0.6549395
[46,] 0.35280559 0.7056112 0.6471944
[47,] 0.27883774 0.5576755 0.7211623
[48,] 0.23856298 0.4771260 0.7614370
[49,] 0.82199104 0.3560179 0.1780090
[50,] 0.75189455 0.4962109 0.2481054
[51,] 0.68399036 0.6320193 0.3160096
[52,] 0.60886064 0.7822787 0.3911394
[53,] 0.53393346 0.9321331 0.4660665
[54,] 0.42664123 0.8532825 0.5733588
[55,] 0.33211365 0.6642273 0.6678864
[56,] 0.27561806 0.5512361 0.7243819
[57,] 0.44920949 0.8984190 0.5507905
> postscript(file="/var/www/rcomp/tmp/1t9481293631225.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/www/rcomp/tmp/2t9481293631225.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/www/rcomp/tmp/3mj3t1293631225.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/www/rcomp/tmp/4mj3t1293631225.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/www/rcomp/tmp/5mj3t1293631225.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 = 72
Frequency = 1
1 2 3 4 5 6
834.194286 479.767264 1381.582143 175.281478 1081.841159 1130.678371
7 8 9 10 11 12
-101.033200 -245.166619 735.182162 1281.456571 947.517420 6.508883
13 14 15 16 17 18
-161.670999 -266.145998 50.588265 -39.373095 794.572341 1162.985096
19 20 21 22 23 24
-661.528647 -193.507439 766.394304 652.461609 478.223456 -409.042162
25 26 27 28 29 30
-343.330513 -620.217178 382.712694 -1585.282777 976.153751 938.335720
31 32 33 34 35 36
-252.249680 -566.997196 442.268291 383.981836 183.639938 -1092.789610
37 38 39 40 41 42
-14.235321 -766.308414 -340.763522 -1107.801832 897.325037 222.570440
43 44 45 46 47 48
-854.345037 -1237.207789 -147.982572 554.385297 119.677194 -1061.543507
49 50 51 52 53 54
22.247973 -745.788892 -1684.734986 -135.601051 72.340573 674.453411
55 56 57 58 59 60
-56.288249 -1819.871879 890.454485 1035.144055 -352.808909 -721.207750
61 62 63 64 65 66
-275.104120 -619.225967 770.933946 -610.467952 504.029434 737.094890
67 68 69 70 71 72
-805.461126 -616.949907 636.097433 212.403444 -611.195192 -1492.255562
> postscript(file="/var/www/rcomp/tmp/6wskd1293631225.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 834.194286 NA
1 479.767264 834.194286
2 1381.582143 479.767264
3 175.281478 1381.582143
4 1081.841159 175.281478
5 1130.678371 1081.841159
6 -101.033200 1130.678371
7 -245.166619 -101.033200
8 735.182162 -245.166619
9 1281.456571 735.182162
10 947.517420 1281.456571
11 6.508883 947.517420
12 -161.670999 6.508883
13 -266.145998 -161.670999
14 50.588265 -266.145998
15 -39.373095 50.588265
16 794.572341 -39.373095
17 1162.985096 794.572341
18 -661.528647 1162.985096
19 -193.507439 -661.528647
20 766.394304 -193.507439
21 652.461609 766.394304
22 478.223456 652.461609
23 -409.042162 478.223456
24 -343.330513 -409.042162
25 -620.217178 -343.330513
26 382.712694 -620.217178
27 -1585.282777 382.712694
28 976.153751 -1585.282777
29 938.335720 976.153751
30 -252.249680 938.335720
31 -566.997196 -252.249680
32 442.268291 -566.997196
33 383.981836 442.268291
34 183.639938 383.981836
35 -1092.789610 183.639938
36 -14.235321 -1092.789610
37 -766.308414 -14.235321
38 -340.763522 -766.308414
39 -1107.801832 -340.763522
40 897.325037 -1107.801832
41 222.570440 897.325037
42 -854.345037 222.570440
43 -1237.207789 -854.345037
44 -147.982572 -1237.207789
45 554.385297 -147.982572
46 119.677194 554.385297
47 -1061.543507 119.677194
48 22.247973 -1061.543507
49 -745.788892 22.247973
50 -1684.734986 -745.788892
51 -135.601051 -1684.734986
52 72.340573 -135.601051
53 674.453411 72.340573
54 -56.288249 674.453411
55 -1819.871879 -56.288249
56 890.454485 -1819.871879
57 1035.144055 890.454485
58 -352.808909 1035.144055
59 -721.207750 -352.808909
60 -275.104120 -721.207750
61 -619.225967 -275.104120
62 770.933946 -619.225967
63 -610.467952 770.933946
64 504.029434 -610.467952
65 737.094890 504.029434
66 -805.461126 737.094890
67 -616.949907 -805.461126
68 636.097433 -616.949907
69 212.403444 636.097433
70 -611.195192 212.403444
71 -1492.255562 -611.195192
72 NA -1492.255562
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 479.767264 834.194286
[2,] 1381.582143 479.767264
[3,] 175.281478 1381.582143
[4,] 1081.841159 175.281478
[5,] 1130.678371 1081.841159
[6,] -101.033200 1130.678371
[7,] -245.166619 -101.033200
[8,] 735.182162 -245.166619
[9,] 1281.456571 735.182162
[10,] 947.517420 1281.456571
[11,] 6.508883 947.517420
[12,] -161.670999 6.508883
[13,] -266.145998 -161.670999
[14,] 50.588265 -266.145998
[15,] -39.373095 50.588265
[16,] 794.572341 -39.373095
[17,] 1162.985096 794.572341
[18,] -661.528647 1162.985096
[19,] -193.507439 -661.528647
[20,] 766.394304 -193.507439
[21,] 652.461609 766.394304
[22,] 478.223456 652.461609
[23,] -409.042162 478.223456
[24,] -343.330513 -409.042162
[25,] -620.217178 -343.330513
[26,] 382.712694 -620.217178
[27,] -1585.282777 382.712694
[28,] 976.153751 -1585.282777
[29,] 938.335720 976.153751
[30,] -252.249680 938.335720
[31,] -566.997196 -252.249680
[32,] 442.268291 -566.997196
[33,] 383.981836 442.268291
[34,] 183.639938 383.981836
[35,] -1092.789610 183.639938
[36,] -14.235321 -1092.789610
[37,] -766.308414 -14.235321
[38,] -340.763522 -766.308414
[39,] -1107.801832 -340.763522
[40,] 897.325037 -1107.801832
[41,] 222.570440 897.325037
[42,] -854.345037 222.570440
[43,] -1237.207789 -854.345037
[44,] -147.982572 -1237.207789
[45,] 554.385297 -147.982572
[46,] 119.677194 554.385297
[47,] -1061.543507 119.677194
[48,] 22.247973 -1061.543507
[49,] -745.788892 22.247973
[50,] -1684.734986 -745.788892
[51,] -135.601051 -1684.734986
[52,] 72.340573 -135.601051
[53,] 674.453411 72.340573
[54,] -56.288249 674.453411
[55,] -1819.871879 -56.288249
[56,] 890.454485 -1819.871879
[57,] 1035.144055 890.454485
[58,] -352.808909 1035.144055
[59,] -721.207750 -352.808909
[60,] -275.104120 -721.207750
[61,] -619.225967 -275.104120
[62,] 770.933946 -619.225967
[63,] -610.467952 770.933946
[64,] 504.029434 -610.467952
[65,] 737.094890 504.029434
[66,] -805.461126 737.094890
[67,] -616.949907 -805.461126
[68,] 636.097433 -616.949907
[69,] 212.403444 636.097433
[70,] -611.195192 212.403444
[71,] -1492.255562 -611.195192
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 479.767264 834.194286
2 1381.582143 479.767264
3 175.281478 1381.582143
4 1081.841159 175.281478
5 1130.678371 1081.841159
6 -101.033200 1130.678371
7 -245.166619 -101.033200
8 735.182162 -245.166619
9 1281.456571 735.182162
10 947.517420 1281.456571
11 6.508883 947.517420
12 -161.670999 6.508883
13 -266.145998 -161.670999
14 50.588265 -266.145998
15 -39.373095 50.588265
16 794.572341 -39.373095
17 1162.985096 794.572341
18 -661.528647 1162.985096
19 -193.507439 -661.528647
20 766.394304 -193.507439
21 652.461609 766.394304
22 478.223456 652.461609
23 -409.042162 478.223456
24 -343.330513 -409.042162
25 -620.217178 -343.330513
26 382.712694 -620.217178
27 -1585.282777 382.712694
28 976.153751 -1585.282777
29 938.335720 976.153751
30 -252.249680 938.335720
31 -566.997196 -252.249680
32 442.268291 -566.997196
33 383.981836 442.268291
34 183.639938 383.981836
35 -1092.789610 183.639938
36 -14.235321 -1092.789610
37 -766.308414 -14.235321
38 -340.763522 -766.308414
39 -1107.801832 -340.763522
40 897.325037 -1107.801832
41 222.570440 897.325037
42 -854.345037 222.570440
43 -1237.207789 -854.345037
44 -147.982572 -1237.207789
45 554.385297 -147.982572
46 119.677194 554.385297
47 -1061.543507 119.677194
48 22.247973 -1061.543507
49 -745.788892 22.247973
50 -1684.734986 -745.788892
51 -135.601051 -1684.734986
52 72.340573 -135.601051
53 674.453411 72.340573
54 -56.288249 674.453411
55 -1819.871879 -56.288249
56 890.454485 -1819.871879
57 1035.144055 890.454485
58 -352.808909 1035.144055
59 -721.207750 -352.808909
60 -275.104120 -721.207750
61 -619.225967 -275.104120
62 770.933946 -619.225967
63 -610.467952 770.933946
64 504.029434 -610.467952
65 737.094890 504.029434
66 -805.461126 737.094890
67 -616.949907 -805.461126
68 636.097433 -616.949907
69 212.403444 636.097433
70 -611.195192 212.403444
71 -1492.255562 -611.195192
> 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/www/rcomp/tmp/7wskd1293631225.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/www/rcomp/tmp/8pjjy1293631225.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/www/rcomp/tmp/9pjjy1293631225.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/www/rcomp/tmp/10ia1j1293631225.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/www/rcomp/tmp/11lbh71293631225.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/www/rcomp/tmp/12pcgd1293631225.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/www/rcomp/tmp/13klv41293631225.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/www/rcomp/tmp/146mus1293631225.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/www/rcomp/tmp/15gdbu1293631225.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/www/rcomp/tmp/16v5931293631225.tab")
+ }
> try(system("convert tmp/1t9481293631225.ps tmp/1t9481293631225.png",intern=TRUE))
character(0)
> try(system("convert tmp/2t9481293631225.ps tmp/2t9481293631225.png",intern=TRUE))
character(0)
> try(system("convert tmp/3mj3t1293631225.ps tmp/3mj3t1293631225.png",intern=TRUE))
character(0)
> try(system("convert tmp/4mj3t1293631225.ps tmp/4mj3t1293631225.png",intern=TRUE))
character(0)
> try(system("convert tmp/5mj3t1293631225.ps tmp/5mj3t1293631225.png",intern=TRUE))
character(0)
> try(system("convert tmp/6wskd1293631225.ps tmp/6wskd1293631225.png",intern=TRUE))
character(0)
> try(system("convert tmp/7wskd1293631225.ps tmp/7wskd1293631225.png",intern=TRUE))
character(0)
> try(system("convert tmp/8pjjy1293631225.ps tmp/8pjjy1293631225.png",intern=TRUE))
character(0)
> try(system("convert tmp/9pjjy1293631225.ps tmp/9pjjy1293631225.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ia1j1293631225.ps tmp/10ia1j1293631225.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.200 1.780 4.979