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(5124
+ ,119.880
+ ,98
+ ,4742
+ ,131.468
+ ,107
+ ,5434
+ ,155.089
+ ,101
+ ,5684
+ ,149.581
+ ,114
+ ,6332
+ ,122.788
+ ,118
+ ,6334
+ ,143.900
+ ,123
+ ,5636
+ ,112.115
+ ,137
+ ,5940
+ ,109.600
+ ,102
+ ,6195
+ ,117.446
+ ,136
+ ,6022
+ ,118.456
+ ,116
+ ,4535
+ ,101.901
+ ,108
+ ,4320
+ ,89.940
+ ,95
+ ,4872
+ ,129.143
+ ,97
+ ,4662
+ ,126.102
+ ,73
+ ,4663
+ ,143.048
+ ,78
+ ,5491
+ ,142.258
+ ,90
+ ,6018
+ ,131.011
+ ,97
+ ,6393
+ ,146.471
+ ,122
+ ,5610
+ ,114.073
+ ,101
+ ,5777
+ ,114.642
+ ,76
+ ,6094
+ ,118.226
+ ,98
+ ,6478
+ ,111.338
+ ,98
+ ,5216
+ ,108.701
+ ,79
+ ,5201
+ ,80.512
+ ,80
+ ,4784
+ ,146.865
+ ,70
+ ,4205
+ ,137.179
+ ,87
+ ,4681
+ ,166.536
+ ,85
+ ,4896
+ ,137.070
+ ,83
+ ,5752
+ ,127.090
+ ,83
+ ,6452
+ ,139.966
+ ,86
+ ,5995
+ ,122.243
+ ,96
+ ,5601
+ ,109.097
+ ,78
+ ,6119
+ ,116.591
+ ,119
+ ,6569
+ ,111.964
+ ,98
+ ,5798
+ ,109.754
+ ,88
+ ,5492
+ ,77.609
+ ,102
+ ,5018
+ ,138.445
+ ,75
+ ,4773
+ ,127.901
+ ,75
+ ,5502
+ ,156.615
+ ,89
+ ,5908
+ ,133.264
+ ,95
+ ,5902
+ ,143.521
+ ,86
+ ,6125
+ ,152.139
+ ,89
+ ,5419
+ ,131.523
+ ,101
+ ,5559
+ ,113.925
+ ,88
+ ,5962
+ ,86.495
+ ,82
+ ,6023
+ ,127.877
+ ,95
+ ,5346
+ ,107.017
+ ,100
+ ,5379
+ ,78.716
+ ,96
+ ,4859
+ ,138.278
+ ,71
+ ,5156
+ ,144.238
+ ,74
+ ,5010
+ ,143.679
+ ,76
+ ,5508
+ ,159.932
+ ,69
+ ,6426
+ ,136.781
+ ,83
+ ,6043
+ ,148.173
+ ,84
+ ,5499
+ ,125.673
+ ,100
+ ,5191
+ ,105.573
+ ,90
+ ,5790
+ ,122.405
+ ,70
+ ,5949
+ ,128.045
+ ,82
+ ,5219
+ ,94.467
+ ,79
+ ,4729
+ ,85.573
+ ,66)
+ ,dim=c(3
+ ,60)
+ ,dimnames=list(c('verkeersongevallen'
+ ,'auto-inschrijvingen'
+ ,'verkeersdoden')
+ ,1:60))
> y <- array(NA,dim=c(3,60),dimnames=list(c('verkeersongevallen','auto-inschrijvingen','verkeersdoden'),1:60))
> 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 = '1'
> #'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
verkeersongevallen auto-inschrijvingen verkeersdoden
1 5124 119.880 98
2 4742 131.468 107
3 5434 155.089 101
4 5684 149.581 114
5 6332 122.788 118
6 6334 143.900 123
7 5636 112.115 137
8 5940 109.600 102
9 6195 117.446 136
10 6022 118.456 116
11 4535 101.901 108
12 4320 89.940 95
13 4872 129.143 97
14 4662 126.102 73
15 4663 143.048 78
16 5491 142.258 90
17 6018 131.011 97
18 6393 146.471 122
19 5610 114.073 101
20 5777 114.642 76
21 6094 118.226 98
22 6478 111.338 98
23 5216 108.701 79
24 5201 80.512 80
25 4784 146.865 70
26 4205 137.179 87
27 4681 166.536 85
28 4896 137.070 83
29 5752 127.090 83
30 6452 139.966 86
31 5995 122.243 96
32 5601 109.097 78
33 6119 116.591 119
34 6569 111.964 98
35 5798 109.754 88
36 5492 77.609 102
37 5018 138.445 75
38 4773 127.901 75
39 5502 156.615 89
40 5908 133.264 95
41 5902 143.521 86
42 6125 152.139 89
43 5419 131.523 101
44 5559 113.925 88
45 5962 86.495 82
46 6023 127.877 95
47 5346 107.017 100
48 5379 78.716 96
49 4859 138.278 71
50 5156 144.238 74
51 5010 143.679 76
52 5508 159.932 69
53 6426 136.781 83
54 6043 148.173 84
55 5499 125.673 100
56 5191 105.573 90
57 5790 122.405 70
58 5949 128.045 82
59 5219 94.467 79
60 4729 85.573 66
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `auto-inschrijvingen` verkeersdoden
3889.451 2.173 14.784
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1268.72 -316.95 -40.82 412.56 1012.28
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3889.451 622.897 6.244 5.71e-08 ***
`auto-inschrijvingen` 2.173 3.454 0.629 0.53177
verkeersdoden 14.784 4.453 3.320 0.00158 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 553.6 on 57 degrees of freedom
Multiple R-squared: 0.1636, Adjusted R-squared: 0.1342
F-statistic: 5.573 on 2 and 57 DF, p-value: 0.006155
> 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.4894215 0.97884303 0.51057849
[2,] 0.5923764 0.81524720 0.40762360
[3,] 0.6359461 0.72810784 0.36405392
[4,] 0.5089158 0.98216835 0.49108418
[5,] 0.4116997 0.82339938 0.58830031
[6,] 0.6211486 0.75770284 0.37885142
[7,] 0.6555225 0.68895507 0.34447754
[8,] 0.6280678 0.74386434 0.37193217
[9,] 0.5753937 0.84921262 0.42460631
[10,] 0.5240200 0.95195992 0.47597996
[11,] 0.4712534 0.94250680 0.52874660
[12,] 0.5494274 0.90114529 0.45057265
[13,] 0.4753298 0.95065962 0.52467019
[14,] 0.4452909 0.89058180 0.55470910
[15,] 0.6431609 0.71367820 0.35683910
[16,] 0.6787633 0.64247345 0.32123672
[17,] 0.8225980 0.35480400 0.17740200
[18,] 0.7718029 0.45639420 0.22819710
[19,] 0.7191783 0.56164331 0.28082166
[20,] 0.6794809 0.64103814 0.32051907
[21,] 0.8941449 0.21171012 0.10585506
[22,] 0.9395383 0.12092341 0.06046171
[23,] 0.9450401 0.10991984 0.05495992
[24,] 0.9385363 0.12292732 0.06146366
[25,] 0.9784178 0.04316443 0.02158221
[26,] 0.9727450 0.05451002 0.02725501
[27,] 0.9648422 0.07031556 0.03515778
[28,] 0.9483540 0.10329196 0.05164598
[29,] 0.9789202 0.04215954 0.02107977
[30,] 0.9721734 0.05565327 0.02782663
[31,] 0.9568361 0.08632779 0.04316390
[32,] 0.9465805 0.10683907 0.05341953
[33,] 0.9541009 0.09179825 0.04589912
[34,] 0.9393358 0.12132832 0.06066416
[35,] 0.9151907 0.16961858 0.08480929
[36,] 0.8910011 0.21799789 0.10899895
[37,] 0.8796343 0.24073145 0.12036573
[38,] 0.8555483 0.28890336 0.14445168
[39,] 0.7962115 0.40757695 0.20378848
[40,] 0.8565884 0.28682327 0.14341164
[41,] 0.8195954 0.36080917 0.18040459
[42,] 0.7679778 0.46404446 0.23202223
[43,] 0.6787884 0.64242314 0.32121157
[44,] 0.6849964 0.63000711 0.31500355
[45,] 0.6563885 0.68722292 0.34361146
[46,] 0.7784160 0.44316805 0.22158403
[47,] 0.9567476 0.08650476 0.04325238
[48,] 0.9780572 0.04388555 0.02194277
[49,] 0.9590563 0.08188735 0.04094368
> postscript(file="/var/www/rcomp/tmp/17a7q1324642408.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/2gtnz1324642408.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/3p0f51324642408.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/4k5ft1324642408.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/5qnaw1324642408.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 = 60
Frequency = 1
1 2 3 4 5 6
-474.74888 -1014.98260 -285.61164 -215.82858 431.26060 313.46467
7 8 9 10 11 12
-522.43371 304.45618 39.76514 160.24152 -1172.51461 -1169.33609
13 14 15 16 17 18
-732.09459 -580.68082 -690.42368 -38.10966 409.84609 381.66124
19 20 21 22 23 24
-20.48047 514.87204 498.84540 897.81360 -77.56835 -46.09487
25 26 27 28 29 30
-459.44986 -1268.72189 -826.94998 -518.35079 359.33658 987.00528
31 32 33 34 35 36
420.68324 321.35466 216.94364 987.45325 369.09135 -74.02472
37 38 39 40 41 42
-281.07030 -503.15731 -43.52506 324.51725 429.27997 589.20166
43 44 45 46 47 48
-249.40077 121.02742 672.33645 451.22365 -254.36363 -100.72897
49 50 51 52 53 54
-380.57316 -140.87541 -315.22778 250.93802 1012.27723 589.73791
55 56 57 58 59 60
-141.90467 -258.39010 599.70375 569.04485 -43.63669 -322.12301
> postscript(file="/var/www/rcomp/tmp/6n0zl1324642408.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -474.74888 NA
1 -1014.98260 -474.74888
2 -285.61164 -1014.98260
3 -215.82858 -285.61164
4 431.26060 -215.82858
5 313.46467 431.26060
6 -522.43371 313.46467
7 304.45618 -522.43371
8 39.76514 304.45618
9 160.24152 39.76514
10 -1172.51461 160.24152
11 -1169.33609 -1172.51461
12 -732.09459 -1169.33609
13 -580.68082 -732.09459
14 -690.42368 -580.68082
15 -38.10966 -690.42368
16 409.84609 -38.10966
17 381.66124 409.84609
18 -20.48047 381.66124
19 514.87204 -20.48047
20 498.84540 514.87204
21 897.81360 498.84540
22 -77.56835 897.81360
23 -46.09487 -77.56835
24 -459.44986 -46.09487
25 -1268.72189 -459.44986
26 -826.94998 -1268.72189
27 -518.35079 -826.94998
28 359.33658 -518.35079
29 987.00528 359.33658
30 420.68324 987.00528
31 321.35466 420.68324
32 216.94364 321.35466
33 987.45325 216.94364
34 369.09135 987.45325
35 -74.02472 369.09135
36 -281.07030 -74.02472
37 -503.15731 -281.07030
38 -43.52506 -503.15731
39 324.51725 -43.52506
40 429.27997 324.51725
41 589.20166 429.27997
42 -249.40077 589.20166
43 121.02742 -249.40077
44 672.33645 121.02742
45 451.22365 672.33645
46 -254.36363 451.22365
47 -100.72897 -254.36363
48 -380.57316 -100.72897
49 -140.87541 -380.57316
50 -315.22778 -140.87541
51 250.93802 -315.22778
52 1012.27723 250.93802
53 589.73791 1012.27723
54 -141.90467 589.73791
55 -258.39010 -141.90467
56 599.70375 -258.39010
57 569.04485 599.70375
58 -43.63669 569.04485
59 -322.12301 -43.63669
60 NA -322.12301
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1014.98260 -474.74888
[2,] -285.61164 -1014.98260
[3,] -215.82858 -285.61164
[4,] 431.26060 -215.82858
[5,] 313.46467 431.26060
[6,] -522.43371 313.46467
[7,] 304.45618 -522.43371
[8,] 39.76514 304.45618
[9,] 160.24152 39.76514
[10,] -1172.51461 160.24152
[11,] -1169.33609 -1172.51461
[12,] -732.09459 -1169.33609
[13,] -580.68082 -732.09459
[14,] -690.42368 -580.68082
[15,] -38.10966 -690.42368
[16,] 409.84609 -38.10966
[17,] 381.66124 409.84609
[18,] -20.48047 381.66124
[19,] 514.87204 -20.48047
[20,] 498.84540 514.87204
[21,] 897.81360 498.84540
[22,] -77.56835 897.81360
[23,] -46.09487 -77.56835
[24,] -459.44986 -46.09487
[25,] -1268.72189 -459.44986
[26,] -826.94998 -1268.72189
[27,] -518.35079 -826.94998
[28,] 359.33658 -518.35079
[29,] 987.00528 359.33658
[30,] 420.68324 987.00528
[31,] 321.35466 420.68324
[32,] 216.94364 321.35466
[33,] 987.45325 216.94364
[34,] 369.09135 987.45325
[35,] -74.02472 369.09135
[36,] -281.07030 -74.02472
[37,] -503.15731 -281.07030
[38,] -43.52506 -503.15731
[39,] 324.51725 -43.52506
[40,] 429.27997 324.51725
[41,] 589.20166 429.27997
[42,] -249.40077 589.20166
[43,] 121.02742 -249.40077
[44,] 672.33645 121.02742
[45,] 451.22365 672.33645
[46,] -254.36363 451.22365
[47,] -100.72897 -254.36363
[48,] -380.57316 -100.72897
[49,] -140.87541 -380.57316
[50,] -315.22778 -140.87541
[51,] 250.93802 -315.22778
[52,] 1012.27723 250.93802
[53,] 589.73791 1012.27723
[54,] -141.90467 589.73791
[55,] -258.39010 -141.90467
[56,] 599.70375 -258.39010
[57,] 569.04485 599.70375
[58,] -43.63669 569.04485
[59,] -322.12301 -43.63669
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1014.98260 -474.74888
2 -285.61164 -1014.98260
3 -215.82858 -285.61164
4 431.26060 -215.82858
5 313.46467 431.26060
6 -522.43371 313.46467
7 304.45618 -522.43371
8 39.76514 304.45618
9 160.24152 39.76514
10 -1172.51461 160.24152
11 -1169.33609 -1172.51461
12 -732.09459 -1169.33609
13 -580.68082 -732.09459
14 -690.42368 -580.68082
15 -38.10966 -690.42368
16 409.84609 -38.10966
17 381.66124 409.84609
18 -20.48047 381.66124
19 514.87204 -20.48047
20 498.84540 514.87204
21 897.81360 498.84540
22 -77.56835 897.81360
23 -46.09487 -77.56835
24 -459.44986 -46.09487
25 -1268.72189 -459.44986
26 -826.94998 -1268.72189
27 -518.35079 -826.94998
28 359.33658 -518.35079
29 987.00528 359.33658
30 420.68324 987.00528
31 321.35466 420.68324
32 216.94364 321.35466
33 987.45325 216.94364
34 369.09135 987.45325
35 -74.02472 369.09135
36 -281.07030 -74.02472
37 -503.15731 -281.07030
38 -43.52506 -503.15731
39 324.51725 -43.52506
40 429.27997 324.51725
41 589.20166 429.27997
42 -249.40077 589.20166
43 121.02742 -249.40077
44 672.33645 121.02742
45 451.22365 672.33645
46 -254.36363 451.22365
47 -100.72897 -254.36363
48 -380.57316 -100.72897
49 -140.87541 -380.57316
50 -315.22778 -140.87541
51 250.93802 -315.22778
52 1012.27723 250.93802
53 589.73791 1012.27723
54 -141.90467 589.73791
55 -258.39010 -141.90467
56 599.70375 -258.39010
57 569.04485 599.70375
58 -43.63669 569.04485
59 -322.12301 -43.63669
> 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/7c0tj1324642408.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/8lfps1324642408.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/9qwj61324642408.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/10p66f1324642408.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/117llr1324642408.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/12r0li1324642408.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/13ovec1324642408.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/14wex81324642408.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/15err01324642408.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/166htx1324642408.tab")
+ }
>
> try(system("convert tmp/17a7q1324642408.ps tmp/17a7q1324642408.png",intern=TRUE))
character(0)
> try(system("convert tmp/2gtnz1324642408.ps tmp/2gtnz1324642408.png",intern=TRUE))
character(0)
> try(system("convert tmp/3p0f51324642408.ps tmp/3p0f51324642408.png",intern=TRUE))
character(0)
> try(system("convert tmp/4k5ft1324642408.ps tmp/4k5ft1324642408.png",intern=TRUE))
character(0)
> try(system("convert tmp/5qnaw1324642408.ps tmp/5qnaw1324642408.png",intern=TRUE))
character(0)
> try(system("convert tmp/6n0zl1324642408.ps tmp/6n0zl1324642408.png",intern=TRUE))
character(0)
> try(system("convert tmp/7c0tj1324642408.ps tmp/7c0tj1324642408.png",intern=TRUE))
character(0)
> try(system("convert tmp/8lfps1324642408.ps tmp/8lfps1324642408.png",intern=TRUE))
character(0)
> try(system("convert tmp/9qwj61324642408.ps tmp/9qwj61324642408.png",intern=TRUE))
character(0)
> try(system("convert tmp/10p66f1324642408.ps tmp/10p66f1324642408.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.940 0.360 5.267