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 = '3'
> 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
veldgoal hoogte gewicht vrijeworp puntpergame
1 0.442 6.8 225 0.672 9.2
2 0.435 6.3 180 0.797 11.7
3 0.456 6.4 190 0.761 15.8
4 0.416 6.2 180 0.651 8.6
5 0.449 6.9 205 0.900 23.2
6 0.431 6.4 225 0.780 27.4
7 0.487 6.3 185 0.771 9.3
8 0.469 6.8 235 0.750 16.0
9 0.435 6.9 235 0.818 4.7
10 0.480 6.7 210 0.825 12.5
11 0.516 6.9 245 0.632 20.1
12 0.493 6.9 245 0.757 9.1
13 0.374 6.3 185 0.709 8.1
14 0.424 6.1 185 0.782 8.6
15 0.441 6.2 180 0.775 20.3
16 0.503 6.8 220 0.880 25.0
17 0.503 6.5 194 0.833 19.2
18 0.425 7.6 225 0.571 3.3
19 0.371 6.3 210 0.816 11.2
20 0.504 7.1 240 0.714 10.5
21 0.400 6.8 225 0.765 10.1
22 0.482 7.3 263 0.655 7.2
23 0.475 6.4 210 0.244 13.6
24 0.428 6.8 235 0.728 9.0
25 0.559 7.2 230 0.721 24.6
26 0.441 6.4 190 0.757 12.6
27 0.492 6.6 220 0.747 5.6
28 0.402 6.8 210 0.739 8.7
29 0.415 6.1 180 0.713 7.7
30 0.492 6.5 235 0.742 24.1
31 0.484 6.4 185 0.861 11.7
32 0.387 6.0 175 0.721 7.7
33 0.436 6.0 192 0.785 9.6
34 0.482 7.3 263 0.655 7.2
35 0.340 6.1 180 0.821 12.3
36 0.516 6.7 240 0.728 8.9
37 0.475 6.4 210 0.846 13.6
38 0.412 5.8 160 0.813 11.2
39 0.411 6.9 230 0.595 2.8
40 0.407 7.0 245 0.573 3.2
41 0.445 7.3 228 0.726 9.4
42 0.291 5.9 155 0.707 11.9
43 0.449 6.2 200 0.804 15.4
44 0.546 6.8 235 0.784 7.4
45 0.480 7.0 235 0.744 18.9
46 0.359 5.9 105 0.839 7.9
47 0.528 6.1 180 0.790 12.2
48 0.352 5.7 185 0.701 11.0
49 0.414 7.1 245 0.778 2.8
50 0.425 5.8 180 0.872 11.8
51 0.599 7.4 240 0.713 17.1
52 0.482 6.8 225 0.701 11.6
53 0.457 6.8 215 0.734 5.8
54 0.435 7.0 230 0.764 8.3
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) hoogte gewicht vrijeworp puntpergame
0.0370425 0.0353033 0.0005548 0.0320573 0.0033314
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.102639 -0.030365 0.003615 0.021841 0.109770
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0370425 0.1238188 0.299 0.76608
hoogte 0.0353033 0.0246410 1.433 0.15829
gewicht 0.0005548 0.0003779 1.468 0.14842
vrijeworp 0.0320573 0.0668196 0.480 0.63353
puntpergame 0.0033314 0.0010916 3.052 0.00367 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0451 on 49 degrees of freedom
Multiple R-squared: 0.4119, Adjusted R-squared: 0.3639
F-statistic: 8.579 on 4 and 49 DF, p-value: 2.495e-05
> 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.15219414 0.30438829 0.8478059
[2,] 0.20673531 0.41347062 0.7932647
[3,] 0.14743216 0.29486432 0.8525678
[4,] 0.21406378 0.42812756 0.7859362
[5,] 0.16130571 0.32261142 0.8386943
[6,] 0.25746871 0.51493743 0.7425313
[7,] 0.17554943 0.35109885 0.8244506
[8,] 0.11505960 0.23011920 0.8849404
[9,] 0.08965991 0.17931982 0.9103401
[10,] 0.09504339 0.19008678 0.9049566
[11,] 0.07546769 0.15093537 0.9245323
[12,] 0.16989683 0.33979366 0.8301032
[13,] 0.14610183 0.29220366 0.8538982
[14,] 0.18218492 0.36436985 0.8178151
[15,] 0.13110219 0.26220438 0.8688978
[16,] 0.15318760 0.30637521 0.8468124
[17,] 0.12747446 0.25494893 0.8725255
[18,] 0.12082837 0.24165674 0.8791716
[19,] 0.08267335 0.16534670 0.9173266
[20,] 0.11406076 0.22812152 0.8859392
[21,] 0.11618724 0.23237448 0.8838128
[22,] 0.08783282 0.17566565 0.9121672
[23,] 0.06229681 0.12459362 0.9377032
[24,] 0.06122246 0.12244491 0.9387775
[25,] 0.04266352 0.08532705 0.9573365
[26,] 0.02959394 0.05918789 0.9704061
[27,] 0.01769275 0.03538550 0.9823073
[28,] 0.06389741 0.12779483 0.9361026
[29,] 0.07646842 0.15293683 0.9235316
[30,] 0.05381776 0.10763552 0.9461822
[31,] 0.03417763 0.06835527 0.9658224
[32,] 0.02454501 0.04909001 0.9754550
[33,] 0.01906483 0.03812965 0.9809352
[34,] 0.01372953 0.02745907 0.9862705
[35,] 0.04927704 0.09855408 0.9507230
[36,] 0.02987177 0.05974355 0.9701282
[37,] 0.08497445 0.16994891 0.9150255
[38,] 0.21637504 0.43275008 0.7836250
[39,] 0.61108971 0.77782058 0.3889103
> postscript(file="/var/www/rcomp/tmp/1g2391322161941.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/2soci1322161941.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/36ej01322161941.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/4npcf1322161941.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/55u0i1322161941.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.012133126 0.011150172 0.010566786 0.010688262 -0.051515769 -0.073105830
7 8 9 10 11 12
0.069204914 -0.015835549 -0.017900740 0.021821210 0.012209764 0.021848207
13 14 15 16 17 18
-0.037809831 0.015244929 -0.007264442 -0.007663308 0.038182187 -0.034482580
19 20 21 22 23 24
-0.068438125 0.025276185 -0.060112732 -0.003660509 0.042372920 -0.032810357
25 26 27 28 29 30
0.035096759 0.006355555 0.057290489 -0.044292798 0.014229313 -0.008972603
31 32 33 34 35 36
0.051794025 -0.007722666 0.023463851 -0.003660509 -0.079557402 0.056278958
37 38 39 40 41 42
0.023074425 0.018051208 -0.025648114 -0.042128205 -0.030846628 -0.102638885
43 44 45 46 47 48
0.005033236 0.088724704 -0.021364972 0.002196754 0.109769516 -0.048032528
49 50 51 52 53 54
-0.043897711 0.016064365 0.087729899 0.018941807 0.017754450 -0.028918927
> postscript(file="/var/www/rcomp/tmp/6ydz01322161941.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.012133126 NA
1 0.011150172 -0.012133126
2 0.010566786 0.011150172
3 0.010688262 0.010566786
4 -0.051515769 0.010688262
5 -0.073105830 -0.051515769
6 0.069204914 -0.073105830
7 -0.015835549 0.069204914
8 -0.017900740 -0.015835549
9 0.021821210 -0.017900740
10 0.012209764 0.021821210
11 0.021848207 0.012209764
12 -0.037809831 0.021848207
13 0.015244929 -0.037809831
14 -0.007264442 0.015244929
15 -0.007663308 -0.007264442
16 0.038182187 -0.007663308
17 -0.034482580 0.038182187
18 -0.068438125 -0.034482580
19 0.025276185 -0.068438125
20 -0.060112732 0.025276185
21 -0.003660509 -0.060112732
22 0.042372920 -0.003660509
23 -0.032810357 0.042372920
24 0.035096759 -0.032810357
25 0.006355555 0.035096759
26 0.057290489 0.006355555
27 -0.044292798 0.057290489
28 0.014229313 -0.044292798
29 -0.008972603 0.014229313
30 0.051794025 -0.008972603
31 -0.007722666 0.051794025
32 0.023463851 -0.007722666
33 -0.003660509 0.023463851
34 -0.079557402 -0.003660509
35 0.056278958 -0.079557402
36 0.023074425 0.056278958
37 0.018051208 0.023074425
38 -0.025648114 0.018051208
39 -0.042128205 -0.025648114
40 -0.030846628 -0.042128205
41 -0.102638885 -0.030846628
42 0.005033236 -0.102638885
43 0.088724704 0.005033236
44 -0.021364972 0.088724704
45 0.002196754 -0.021364972
46 0.109769516 0.002196754
47 -0.048032528 0.109769516
48 -0.043897711 -0.048032528
49 0.016064365 -0.043897711
50 0.087729899 0.016064365
51 0.018941807 0.087729899
52 0.017754450 0.018941807
53 -0.028918927 0.017754450
54 NA -0.028918927
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.011150172 -0.012133126
[2,] 0.010566786 0.011150172
[3,] 0.010688262 0.010566786
[4,] -0.051515769 0.010688262
[5,] -0.073105830 -0.051515769
[6,] 0.069204914 -0.073105830
[7,] -0.015835549 0.069204914
[8,] -0.017900740 -0.015835549
[9,] 0.021821210 -0.017900740
[10,] 0.012209764 0.021821210
[11,] 0.021848207 0.012209764
[12,] -0.037809831 0.021848207
[13,] 0.015244929 -0.037809831
[14,] -0.007264442 0.015244929
[15,] -0.007663308 -0.007264442
[16,] 0.038182187 -0.007663308
[17,] -0.034482580 0.038182187
[18,] -0.068438125 -0.034482580
[19,] 0.025276185 -0.068438125
[20,] -0.060112732 0.025276185
[21,] -0.003660509 -0.060112732
[22,] 0.042372920 -0.003660509
[23,] -0.032810357 0.042372920
[24,] 0.035096759 -0.032810357
[25,] 0.006355555 0.035096759
[26,] 0.057290489 0.006355555
[27,] -0.044292798 0.057290489
[28,] 0.014229313 -0.044292798
[29,] -0.008972603 0.014229313
[30,] 0.051794025 -0.008972603
[31,] -0.007722666 0.051794025
[32,] 0.023463851 -0.007722666
[33,] -0.003660509 0.023463851
[34,] -0.079557402 -0.003660509
[35,] 0.056278958 -0.079557402
[36,] 0.023074425 0.056278958
[37,] 0.018051208 0.023074425
[38,] -0.025648114 0.018051208
[39,] -0.042128205 -0.025648114
[40,] -0.030846628 -0.042128205
[41,] -0.102638885 -0.030846628
[42,] 0.005033236 -0.102638885
[43,] 0.088724704 0.005033236
[44,] -0.021364972 0.088724704
[45,] 0.002196754 -0.021364972
[46,] 0.109769516 0.002196754
[47,] -0.048032528 0.109769516
[48,] -0.043897711 -0.048032528
[49,] 0.016064365 -0.043897711
[50,] 0.087729899 0.016064365
[51,] 0.018941807 0.087729899
[52,] 0.017754450 0.018941807
[53,] -0.028918927 0.017754450
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.011150172 -0.012133126
2 0.010566786 0.011150172
3 0.010688262 0.010566786
4 -0.051515769 0.010688262
5 -0.073105830 -0.051515769
6 0.069204914 -0.073105830
7 -0.015835549 0.069204914
8 -0.017900740 -0.015835549
9 0.021821210 -0.017900740
10 0.012209764 0.021821210
11 0.021848207 0.012209764
12 -0.037809831 0.021848207
13 0.015244929 -0.037809831
14 -0.007264442 0.015244929
15 -0.007663308 -0.007264442
16 0.038182187 -0.007663308
17 -0.034482580 0.038182187
18 -0.068438125 -0.034482580
19 0.025276185 -0.068438125
20 -0.060112732 0.025276185
21 -0.003660509 -0.060112732
22 0.042372920 -0.003660509
23 -0.032810357 0.042372920
24 0.035096759 -0.032810357
25 0.006355555 0.035096759
26 0.057290489 0.006355555
27 -0.044292798 0.057290489
28 0.014229313 -0.044292798
29 -0.008972603 0.014229313
30 0.051794025 -0.008972603
31 -0.007722666 0.051794025
32 0.023463851 -0.007722666
33 -0.003660509 0.023463851
34 -0.079557402 -0.003660509
35 0.056278958 -0.079557402
36 0.023074425 0.056278958
37 0.018051208 0.023074425
38 -0.025648114 0.018051208
39 -0.042128205 -0.025648114
40 -0.030846628 -0.042128205
41 -0.102638885 -0.030846628
42 0.005033236 -0.102638885
43 0.088724704 0.005033236
44 -0.021364972 0.088724704
45 0.002196754 -0.021364972
46 0.109769516 0.002196754
47 -0.048032528 0.109769516
48 -0.043897711 -0.048032528
49 0.016064365 -0.043897711
50 0.087729899 0.016064365
51 0.018941807 0.087729899
52 0.017754450 0.018941807
53 -0.028918927 0.017754450
> 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/7red51322161941.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/808sg1322161941.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/9bwq41322161941.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/10k2mq1322161941.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/11g5r41322161941.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/12v4j01322161941.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/13ykjw1322161941.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/14xixo1322161941.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/15xbux1322161941.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/16kmev1322161941.tab")
+ }
>
> try(system("convert tmp/1g2391322161941.ps tmp/1g2391322161941.png",intern=TRUE))
character(0)
> try(system("convert tmp/2soci1322161941.ps tmp/2soci1322161941.png",intern=TRUE))
character(0)
> try(system("convert tmp/36ej01322161941.ps tmp/36ej01322161941.png",intern=TRUE))
character(0)
> try(system("convert tmp/4npcf1322161941.ps tmp/4npcf1322161941.png",intern=TRUE))
character(0)
> try(system("convert tmp/55u0i1322161941.ps tmp/55u0i1322161941.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ydz01322161941.ps tmp/6ydz01322161941.png",intern=TRUE))
character(0)
> try(system("convert tmp/7red51322161941.ps tmp/7red51322161941.png",intern=TRUE))
character(0)
> try(system("convert tmp/808sg1322161941.ps tmp/808sg1322161941.png",intern=TRUE))
character(0)
> try(system("convert tmp/9bwq41322161941.ps tmp/9bwq41322161941.png",intern=TRUE))
character(0)
> try(system("convert tmp/10k2mq1322161941.ps tmp/10k2mq1322161941.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.032 0.668 4.720