R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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.
Natural language support but running in an English locale
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
Attaching package: 'zoo'
The following object(s) are masked from package:base :
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
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.872 -610.650 -3.863 657.960 1381.582
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/html/freestat/rcomp/tmp/1v90i1293626052.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/html/freestat/rcomp/tmp/2v90i1293626052.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/html/freestat/rcomp/tmp/3n0zl1293626052.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/html/freestat/rcomp/tmp/4n0zl1293626052.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/html/freestat/rcomp/tmp/5n0zl1293626052.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/html/freestat/rcomp/tmp/6yrgo1293626052.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/html/freestat/rcomp/tmp/7rjfr1293626052.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/html/freestat/rcomp/tmp/8rjfr1293626052.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/html/freestat/rcomp/tmp/9rjfr1293626052.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/html/freestat/rcomp/tmp/10jsfc1293626052.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/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/html/freestat/rcomp/tmp/115svi1293626052.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/html/freestat/rcomp/tmp/12qtun1293626052.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/html/freestat/rcomp/tmp/13mlae1293626052.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/html/freestat/rcomp/tmp/14q3q21293626052.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/html/freestat/rcomp/tmp/15t4o81293626052.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/html/freestat/rcomp/tmp/16e4nw1293626052.tab")
+ }
> try(system("convert tmp/1v90i1293626052.ps tmp/1v90i1293626052.png",intern=TRUE))
character(0)
> try(system("convert tmp/2v90i1293626052.ps tmp/2v90i1293626052.png",intern=TRUE))
character(0)
> try(system("convert tmp/3n0zl1293626052.ps tmp/3n0zl1293626052.png",intern=TRUE))
character(0)
> try(system("convert tmp/4n0zl1293626052.ps tmp/4n0zl1293626052.png",intern=TRUE))
character(0)
> try(system("convert tmp/5n0zl1293626052.ps tmp/5n0zl1293626052.png",intern=TRUE))
character(0)
> try(system("convert tmp/6yrgo1293626052.ps tmp/6yrgo1293626052.png",intern=TRUE))
character(0)
> try(system("convert tmp/7rjfr1293626052.ps tmp/7rjfr1293626052.png",intern=TRUE))
character(0)
> try(system("convert tmp/8rjfr1293626052.ps tmp/8rjfr1293626052.png",intern=TRUE))
character(0)
> try(system("convert tmp/9rjfr1293626052.ps tmp/9rjfr1293626052.png",intern=TRUE))
character(0)
> try(system("convert tmp/10jsfc1293626052.ps tmp/10jsfc1293626052.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.072 2.553 4.411