R version 2.12.1 (2010-12-16)
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(6.8
+ ,225
+ ,0.442
+ ,0.672
+ ,9.2
+ ,6.3
+ ,180
+ ,0.435
+ ,0.797
+ ,11.7
+ ,6.4
+ ,190
+ ,0.456
+ ,0.761
+ ,15.8
+ ,6.2
+ ,180
+ ,0.416
+ ,0.651
+ ,8.6
+ ,6.9
+ ,205
+ ,0.449
+ ,0.9
+ ,23.2
+ ,6.4
+ ,225
+ ,0.431
+ ,0.78
+ ,27.4
+ ,6.3
+ ,185
+ ,0.487
+ ,0.771
+ ,9.3
+ ,6.8
+ ,235
+ ,0.469
+ ,0.75
+ ,16
+ ,6.9
+ ,235
+ ,0.435
+ ,0.818
+ ,4.7
+ ,6.7
+ ,210
+ ,0.48
+ ,0.825
+ ,12.5
+ ,6.9
+ ,245
+ ,0.516
+ ,0.632
+ ,20.1
+ ,6.9
+ ,245
+ ,0.493
+ ,0.757
+ ,9.1
+ ,6.3
+ ,185
+ ,0.374
+ ,0.709
+ ,8.1
+ ,6.1
+ ,185
+ ,0.424
+ ,0.782
+ ,8.6
+ ,6.2
+ ,180
+ ,0.441
+ ,0.775
+ ,20.3
+ ,6.8
+ ,220
+ ,0.503
+ ,0.88
+ ,25
+ ,6.5
+ ,194
+ ,0.503
+ ,0.833
+ ,19.2
+ ,7.6
+ ,225
+ ,0.425
+ ,0.571
+ ,3.3
+ ,6.3
+ ,210
+ ,0.371
+ ,0.816
+ ,11.2
+ ,7.1
+ ,240
+ ,0.504
+ ,0.714
+ ,10.5
+ ,6.8
+ ,225
+ ,0.4
+ ,0.765
+ ,10.1
+ ,7.3
+ ,263
+ ,0.482
+ ,0.655
+ ,7.2
+ ,6.4
+ ,210
+ ,0.475
+ ,0.244
+ ,13.6
+ ,6.8
+ ,235
+ ,0.428
+ ,0.728
+ ,9
+ ,7.2
+ ,230
+ ,0.559
+ ,0.721
+ ,24.6
+ ,6.4
+ ,190
+ ,0.441
+ ,0.757
+ ,12.6
+ ,6.6
+ ,220
+ ,0.492
+ ,0.747
+ ,5.6
+ ,6.8
+ ,210
+ ,0.402
+ ,0.739
+ ,8.7
+ ,6.1
+ ,180
+ ,0.415
+ ,0.713
+ ,7.7
+ ,6.5
+ ,235
+ ,0.492
+ ,0.742
+ ,24.1
+ ,6.4
+ ,185
+ ,0.484
+ ,0.861
+ ,11.7
+ ,6
+ ,175
+ ,0.387
+ ,0.721
+ ,7.7
+ ,6
+ ,192
+ ,0.436
+ ,0.785
+ ,9.6
+ ,7.3
+ ,263
+ ,0.482
+ ,0.655
+ ,7.2
+ ,6.1
+ ,180
+ ,0.34
+ ,0.821
+ ,12.3
+ ,6.7
+ ,240
+ ,0.516
+ ,0.728
+ ,8.9
+ ,6.4
+ ,210
+ ,0.475
+ ,0.846
+ ,13.6
+ ,5.8
+ ,160
+ ,0.412
+ ,0.813
+ ,11.2
+ ,6.9
+ ,230
+ ,0.411
+ ,0.595
+ ,2.8
+ ,7
+ ,245
+ ,0.407
+ ,0.573
+ ,3.2
+ ,7.3
+ ,228
+ ,0.445
+ ,0.726
+ ,9.4
+ ,5.9
+ ,155
+ ,0.291
+ ,0.707
+ ,11.9
+ ,6.2
+ ,200
+ ,0.449
+ ,0.804
+ ,15.4
+ ,6.8
+ ,235
+ ,0.546
+ ,0.784
+ ,7.4
+ ,7
+ ,235
+ ,0.48
+ ,0.744
+ ,18.9
+ ,5.9
+ ,105
+ ,0.359
+ ,0.839
+ ,7.9
+ ,6.1
+ ,180
+ ,0.528
+ ,0.79
+ ,12.2
+ ,5.7
+ ,185
+ ,0.352
+ ,0.701
+ ,11
+ ,7.1
+ ,245
+ ,0.414
+ ,0.778
+ ,2.8
+ ,5.8
+ ,180
+ ,0.425
+ ,0.872
+ ,11.8
+ ,7.4
+ ,240
+ ,0.599
+ ,0.713
+ ,17.1
+ ,6.8
+ ,225
+ ,0.482
+ ,0.701
+ ,11.6
+ ,6.8
+ ,215
+ ,0.457
+ ,0.734
+ ,5.8
+ ,7
+ ,230
+ ,0.435
+ ,0.764
+ ,8.3)
+ ,dim=c(5
+ ,54)
+ ,dimnames=list(c('hoogte'
+ ,'gewicht'
+ ,'veldgoal'
+ ,'vrijeworp'
+ ,'puntpergame')
+ ,1:54))
> y <- array(NA,dim=c(5,54),dimnames=list(c('hoogte','gewicht','veldgoal','vrijeworp','puntpergame'),1:54))
> 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'
> 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
hoogte gewicht veldgoal vrijeworp puntpergame
1 6.8 225 0.442 0.672 9.2
2 6.3 180 0.435 0.797 11.7
3 6.4 190 0.456 0.761 15.8
4 6.2 180 0.416 0.651 8.6
5 6.9 205 0.449 0.900 23.2
6 6.4 225 0.431 0.780 27.4
7 6.3 185 0.487 0.771 9.3
8 6.8 235 0.469 0.750 16.0
9 6.9 235 0.435 0.818 4.7
10 6.7 210 0.480 0.825 12.5
11 6.9 245 0.516 0.632 20.1
12 6.9 245 0.493 0.757 9.1
13 6.3 185 0.374 0.709 8.1
14 6.1 185 0.424 0.782 8.6
15 6.2 180 0.441 0.775 20.3
16 6.8 220 0.503 0.880 25.0
17 6.5 194 0.503 0.833 19.2
18 7.6 225 0.425 0.571 3.3
19 6.3 210 0.371 0.816 11.2
20 7.1 240 0.504 0.714 10.5
21 6.8 225 0.400 0.765 10.1
22 7.3 263 0.482 0.655 7.2
23 6.4 210 0.475 0.244 13.6
24 6.8 235 0.428 0.728 9.0
25 7.2 230 0.559 0.721 24.6
26 6.4 190 0.441 0.757 12.6
27 6.6 220 0.492 0.747 5.6
28 6.8 210 0.402 0.739 8.7
29 6.1 180 0.415 0.713 7.7
30 6.5 235 0.492 0.742 24.1
31 6.4 185 0.484 0.861 11.7
32 6.0 175 0.387 0.721 7.7
33 6.0 192 0.436 0.785 9.6
34 7.3 263 0.482 0.655 7.2
35 6.1 180 0.340 0.821 12.3
36 6.7 240 0.516 0.728 8.9
37 6.4 210 0.475 0.846 13.6
38 5.8 160 0.412 0.813 11.2
39 6.9 230 0.411 0.595 2.8
40 7.0 245 0.407 0.573 3.2
41 7.3 228 0.445 0.726 9.4
42 5.9 155 0.291 0.707 11.9
43 6.2 200 0.449 0.804 15.4
44 6.8 235 0.546 0.784 7.4
45 7.0 235 0.480 0.744 18.9
46 5.9 105 0.359 0.839 7.9
47 6.1 180 0.528 0.790 12.2
48 5.7 185 0.352 0.701 11.0
49 7.1 245 0.414 0.778 2.8
50 5.8 180 0.425 0.872 11.8
51 7.4 240 0.599 0.713 17.1
52 6.8 225 0.482 0.701 11.6
53 6.8 215 0.457 0.734 5.8
54 7.0 230 0.435 0.764 8.3
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) gewicht veldgoal vrijeworp puntpergame
3.798050 0.011489 1.138890 -0.049316 -0.008273
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.49883 -0.19071 0.01065 0.09117 0.78835
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.798050 0.448437 8.470 3.69e-11 ***
gewicht 0.011489 0.001454 7.899 2.72e-10 ***
veldgoal 1.138890 0.794921 1.433 0.158
vrijeworp -0.049316 0.380348 -0.130 0.897
puntpergame -0.008273 0.006660 -1.242 0.220
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2562 on 49 degrees of freedom
Multiple R-squared: 0.7119, Adjusted R-squared: 0.6883
F-statistic: 30.26 on 4 and 49 DF, p-value: 1.062e-12
> 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.286048075 0.57209615 0.7139519
[2,] 0.305235523 0.61047105 0.6947645
[3,] 0.179743225 0.35948645 0.8202568
[4,] 0.118628635 0.23725727 0.8813714
[5,] 0.074752121 0.14950424 0.9252479
[6,] 0.039224860 0.07844972 0.9607751
[7,] 0.052366150 0.10473230 0.9476339
[8,] 0.028210779 0.05642156 0.9717892
[9,] 0.015033437 0.03006687 0.9849666
[10,] 0.007465951 0.01493190 0.9925340
[11,] 0.561038589 0.87792282 0.4389614
[12,] 0.532156136 0.93568773 0.4678439
[13,] 0.445716043 0.89143209 0.5542840
[14,] 0.363880277 0.72776055 0.6361197
[15,] 0.286166773 0.57233355 0.7138332
[16,] 0.293974204 0.58794841 0.7060258
[17,] 0.231360895 0.46272179 0.7686391
[18,] 0.281519088 0.56303818 0.7184809
[19,] 0.216307554 0.43261511 0.7836924
[20,] 0.207439420 0.41487884 0.7925606
[21,] 0.198719324 0.39743865 0.8012807
[22,] 0.162738845 0.32547769 0.8372612
[23,] 0.174595457 0.34919091 0.8254045
[24,] 0.130458061 0.26091612 0.8695419
[25,] 0.103833898 0.20766780 0.8961661
[26,] 0.144853428 0.28970686 0.8551466
[27,] 0.100128547 0.20025709 0.8998715
[28,] 0.066808043 0.13361609 0.9331920
[29,] 0.082764056 0.16552811 0.9172359
[30,] 0.061132506 0.12226501 0.9388675
[31,] 0.045438522 0.09087704 0.9545615
[32,] 0.027592423 0.05518485 0.9724076
[33,] 0.019095656 0.03819131 0.9809043
[34,] 0.047469306 0.09493861 0.9525307
[35,] 0.028492805 0.05698561 0.9715072
[36,] 0.018616727 0.03723345 0.9813833
[37,] 0.016422955 0.03284591 0.9835770
[38,] 0.024939865 0.04987973 0.9750601
[39,] 0.884221274 0.23155745 0.1157787
> postscript(file="/var/www/rcomp/tmp/1jvg91322162488.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/2euf41322162488.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/3lana1322162488.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/48e391322162488.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/5kuqb1322162488.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 = 54
Frequency = 1
1 2 3 4 5 6
0.022783571 0.074610058 0.067948418 -0.036598486 0.471662791 -0.208787649
7 8 9 10 11 12
-0.063195597 -0.062751397 -0.014164515 0.086688575 -0.103068040 -0.161715953
13 14 15 16 17 18
0.052513344 -0.196694426 0.037842563 0.151732985 0.100144579 0.788351076
19 20 21 22 23 24
-0.200371653 0.092663540 0.082649292 0.023259707 -0.227168630 -0.075055295
25 26 27 28 29 30
0.361914337 0.058359789 -0.202801129 0.239842105 -0.139848044 -0.322326270
31 32 33 34 35 36
0.064515504 -0.150119451 -0.382363033 0.023259707 -0.011047851 -0.333550073
37 38 39 40 41 42
-0.197480669 -0.172762462 0.043897269 -0.021658259 0.489217515 0.123052231
43 44 45 46 47 48
-0.240158447 -0.219919948 0.148417631 0.593475615 -0.227515207 -0.498832904
49 50 51 52 53 54
0.077169863 -0.409475045 0.339023805 -0.001485828 0.095518734 0.170401656
> postscript(file="/var/www/rcomp/tmp/651sw1322162489.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 = 54
Frequency = 1
lag(myerror, k = 1) myerror
0 0.022783571 NA
1 0.074610058 0.022783571
2 0.067948418 0.074610058
3 -0.036598486 0.067948418
4 0.471662791 -0.036598486
5 -0.208787649 0.471662791
6 -0.063195597 -0.208787649
7 -0.062751397 -0.063195597
8 -0.014164515 -0.062751397
9 0.086688575 -0.014164515
10 -0.103068040 0.086688575
11 -0.161715953 -0.103068040
12 0.052513344 -0.161715953
13 -0.196694426 0.052513344
14 0.037842563 -0.196694426
15 0.151732985 0.037842563
16 0.100144579 0.151732985
17 0.788351076 0.100144579
18 -0.200371653 0.788351076
19 0.092663540 -0.200371653
20 0.082649292 0.092663540
21 0.023259707 0.082649292
22 -0.227168630 0.023259707
23 -0.075055295 -0.227168630
24 0.361914337 -0.075055295
25 0.058359789 0.361914337
26 -0.202801129 0.058359789
27 0.239842105 -0.202801129
28 -0.139848044 0.239842105
29 -0.322326270 -0.139848044
30 0.064515504 -0.322326270
31 -0.150119451 0.064515504
32 -0.382363033 -0.150119451
33 0.023259707 -0.382363033
34 -0.011047851 0.023259707
35 -0.333550073 -0.011047851
36 -0.197480669 -0.333550073
37 -0.172762462 -0.197480669
38 0.043897269 -0.172762462
39 -0.021658259 0.043897269
40 0.489217515 -0.021658259
41 0.123052231 0.489217515
42 -0.240158447 0.123052231
43 -0.219919948 -0.240158447
44 0.148417631 -0.219919948
45 0.593475615 0.148417631
46 -0.227515207 0.593475615
47 -0.498832904 -0.227515207
48 0.077169863 -0.498832904
49 -0.409475045 0.077169863
50 0.339023805 -0.409475045
51 -0.001485828 0.339023805
52 0.095518734 -0.001485828
53 0.170401656 0.095518734
54 NA 0.170401656
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.074610058 0.022783571
[2,] 0.067948418 0.074610058
[3,] -0.036598486 0.067948418
[4,] 0.471662791 -0.036598486
[5,] -0.208787649 0.471662791
[6,] -0.063195597 -0.208787649
[7,] -0.062751397 -0.063195597
[8,] -0.014164515 -0.062751397
[9,] 0.086688575 -0.014164515
[10,] -0.103068040 0.086688575
[11,] -0.161715953 -0.103068040
[12,] 0.052513344 -0.161715953
[13,] -0.196694426 0.052513344
[14,] 0.037842563 -0.196694426
[15,] 0.151732985 0.037842563
[16,] 0.100144579 0.151732985
[17,] 0.788351076 0.100144579
[18,] -0.200371653 0.788351076
[19,] 0.092663540 -0.200371653
[20,] 0.082649292 0.092663540
[21,] 0.023259707 0.082649292
[22,] -0.227168630 0.023259707
[23,] -0.075055295 -0.227168630
[24,] 0.361914337 -0.075055295
[25,] 0.058359789 0.361914337
[26,] -0.202801129 0.058359789
[27,] 0.239842105 -0.202801129
[28,] -0.139848044 0.239842105
[29,] -0.322326270 -0.139848044
[30,] 0.064515504 -0.322326270
[31,] -0.150119451 0.064515504
[32,] -0.382363033 -0.150119451
[33,] 0.023259707 -0.382363033
[34,] -0.011047851 0.023259707
[35,] -0.333550073 -0.011047851
[36,] -0.197480669 -0.333550073
[37,] -0.172762462 -0.197480669
[38,] 0.043897269 -0.172762462
[39,] -0.021658259 0.043897269
[40,] 0.489217515 -0.021658259
[41,] 0.123052231 0.489217515
[42,] -0.240158447 0.123052231
[43,] -0.219919948 -0.240158447
[44,] 0.148417631 -0.219919948
[45,] 0.593475615 0.148417631
[46,] -0.227515207 0.593475615
[47,] -0.498832904 -0.227515207
[48,] 0.077169863 -0.498832904
[49,] -0.409475045 0.077169863
[50,] 0.339023805 -0.409475045
[51,] -0.001485828 0.339023805
[52,] 0.095518734 -0.001485828
[53,] 0.170401656 0.095518734
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.074610058 0.022783571
2 0.067948418 0.074610058
3 -0.036598486 0.067948418
4 0.471662791 -0.036598486
5 -0.208787649 0.471662791
6 -0.063195597 -0.208787649
7 -0.062751397 -0.063195597
8 -0.014164515 -0.062751397
9 0.086688575 -0.014164515
10 -0.103068040 0.086688575
11 -0.161715953 -0.103068040
12 0.052513344 -0.161715953
13 -0.196694426 0.052513344
14 0.037842563 -0.196694426
15 0.151732985 0.037842563
16 0.100144579 0.151732985
17 0.788351076 0.100144579
18 -0.200371653 0.788351076
19 0.092663540 -0.200371653
20 0.082649292 0.092663540
21 0.023259707 0.082649292
22 -0.227168630 0.023259707
23 -0.075055295 -0.227168630
24 0.361914337 -0.075055295
25 0.058359789 0.361914337
26 -0.202801129 0.058359789
27 0.239842105 -0.202801129
28 -0.139848044 0.239842105
29 -0.322326270 -0.139848044
30 0.064515504 -0.322326270
31 -0.150119451 0.064515504
32 -0.382363033 -0.150119451
33 0.023259707 -0.382363033
34 -0.011047851 0.023259707
35 -0.333550073 -0.011047851
36 -0.197480669 -0.333550073
37 -0.172762462 -0.197480669
38 0.043897269 -0.172762462
39 -0.021658259 0.043897269
40 0.489217515 -0.021658259
41 0.123052231 0.489217515
42 -0.240158447 0.123052231
43 -0.219919948 -0.240158447
44 0.148417631 -0.219919948
45 0.593475615 0.148417631
46 -0.227515207 0.593475615
47 -0.498832904 -0.227515207
48 0.077169863 -0.498832904
49 -0.409475045 0.077169863
50 0.339023805 -0.409475045
51 -0.001485828 0.339023805
52 0.095518734 -0.001485828
53 0.170401656 0.095518734
> 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/7gyzz1322162489.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/80qgj1322162489.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/9lwt01322162489.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/10r84d1322162489.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/11isor1322162489.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/126r0r1322162489.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/13mr3s1322162489.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/14svp81322162489.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/15to8t1322162489.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/16eod01322162489.tab")
+ }
>
> try(system("convert tmp/1jvg91322162488.ps tmp/1jvg91322162488.png",intern=TRUE))
character(0)
> try(system("convert tmp/2euf41322162488.ps tmp/2euf41322162488.png",intern=TRUE))
character(0)
> try(system("convert tmp/3lana1322162488.ps tmp/3lana1322162488.png",intern=TRUE))
character(0)
> try(system("convert tmp/48e391322162488.ps tmp/48e391322162488.png",intern=TRUE))
character(0)
> try(system("convert tmp/5kuqb1322162488.ps tmp/5kuqb1322162488.png",intern=TRUE))
character(0)
> try(system("convert tmp/651sw1322162489.ps tmp/651sw1322162489.png",intern=TRUE))
character(0)
> try(system("convert tmp/7gyzz1322162489.ps tmp/7gyzz1322162489.png",intern=TRUE))
character(0)
> try(system("convert tmp/80qgj1322162489.ps tmp/80qgj1322162489.png",intern=TRUE))
character(0)
> try(system("convert tmp/9lwt01322162489.ps tmp/9lwt01322162489.png",intern=TRUE))
character(0)
> try(system("convert tmp/10r84d1322162489.ps tmp/10r84d1322162489.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.132 0.708 4.960