R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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.
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(363,14.3,364,14.2,363,15.9,358,15.3,357,15.5,357,15.1,380,15,378,12.1,376,15.8,380,16.9,379,15.1,384,13.7,392,14.8,394,14.7,392,16,396,15.4,392,15,396,15.5,419,15.1,421,11.7,420,16.3,418,16.7,410,15,418,14.9,426,14.6,428,15.3,430,17.9,424,16.4,423,15.4,427,17.9,441,15.9,449,13.9,452,17.8,462,17.9,455,17.4,461,16.7,461,16,463,16.6,462,19.1,456,17.8,455,17.2,456,18.6,472,16.3,472,15.1,471,19.2,465,17.7,459,19.1,465,18,468,17.5,467,17.8,463,21.1,460,17.2,462,19.4,461,19.8,476,17.6,476,16.2,471,19.5,453,19.9,443,20,442,17.3),dim=c(2,60),dimnames=list(c('WK>25j','ExpBE'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('WK>25j','ExpBE'),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 = 'Include Monthly 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
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
WK>25j ExpBE M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 363 14.3 1 0 0 0 0 0 0 0 0 0 0
2 364 14.2 0 1 0 0 0 0 0 0 0 0 0
3 363 15.9 0 0 1 0 0 0 0 0 0 0 0
4 358 15.3 0 0 0 1 0 0 0 0 0 0 0
5 357 15.5 0 0 0 0 1 0 0 0 0 0 0
6 357 15.1 0 0 0 0 0 1 0 0 0 0 0
7 380 15.0 0 0 0 0 0 0 1 0 0 0 0
8 378 12.1 0 0 0 0 0 0 0 1 0 0 0
9 376 15.8 0 0 0 0 0 0 0 0 1 0 0
10 380 16.9 0 0 0 0 0 0 0 0 0 1 0
11 379 15.1 0 0 0 0 0 0 0 0 0 0 1
12 384 13.7 0 0 0 0 0 0 0 0 0 0 0
13 392 14.8 1 0 0 0 0 0 0 0 0 0 0
14 394 14.7 0 1 0 0 0 0 0 0 0 0 0
15 392 16.0 0 0 1 0 0 0 0 0 0 0 0
16 396 15.4 0 0 0 1 0 0 0 0 0 0 0
17 392 15.0 0 0 0 0 1 0 0 0 0 0 0
18 396 15.5 0 0 0 0 0 1 0 0 0 0 0
19 419 15.1 0 0 0 0 0 0 1 0 0 0 0
20 421 11.7 0 0 0 0 0 0 0 1 0 0 0
21 420 16.3 0 0 0 0 0 0 0 0 1 0 0
22 418 16.7 0 0 0 0 0 0 0 0 0 1 0
23 410 15.0 0 0 0 0 0 0 0 0 0 0 1
24 418 14.9 0 0 0 0 0 0 0 0 0 0 0
25 426 14.6 1 0 0 0 0 0 0 0 0 0 0
26 428 15.3 0 1 0 0 0 0 0 0 0 0 0
27 430 17.9 0 0 1 0 0 0 0 0 0 0 0
28 424 16.4 0 0 0 1 0 0 0 0 0 0 0
29 423 15.4 0 0 0 0 1 0 0 0 0 0 0
30 427 17.9 0 0 0 0 0 1 0 0 0 0 0
31 441 15.9 0 0 0 0 0 0 1 0 0 0 0
32 449 13.9 0 0 0 0 0 0 0 1 0 0 0
33 452 17.8 0 0 0 0 0 0 0 0 1 0 0
34 462 17.9 0 0 0 0 0 0 0 0 0 1 0
35 455 17.4 0 0 0 0 0 0 0 0 0 0 1
36 461 16.7 0 0 0 0 0 0 0 0 0 0 0
37 461 16.0 1 0 0 0 0 0 0 0 0 0 0
38 463 16.6 0 1 0 0 0 0 0 0 0 0 0
39 462 19.1 0 0 1 0 0 0 0 0 0 0 0
40 456 17.8 0 0 0 1 0 0 0 0 0 0 0
41 455 17.2 0 0 0 0 1 0 0 0 0 0 0
42 456 18.6 0 0 0 0 0 1 0 0 0 0 0
43 472 16.3 0 0 0 0 0 0 1 0 0 0 0
44 472 15.1 0 0 0 0 0 0 0 1 0 0 0
45 471 19.2 0 0 0 0 0 0 0 0 1 0 0
46 465 17.7 0 0 0 0 0 0 0 0 0 1 0
47 459 19.1 0 0 0 0 0 0 0 0 0 0 1
48 465 18.0 0 0 0 0 0 0 0 0 0 0 0
49 468 17.5 1 0 0 0 0 0 0 0 0 0 0
50 467 17.8 0 1 0 0 0 0 0 0 0 0 0
51 463 21.1 0 0 1 0 0 0 0 0 0 0 0
52 460 17.2 0 0 0 1 0 0 0 0 0 0 0
53 462 19.4 0 0 0 0 1 0 0 0 0 0 0
54 461 19.8 0 0 0 0 0 1 0 0 0 0 0
55 476 17.6 0 0 0 0 0 0 1 0 0 0 0
56 476 16.2 0 0 0 0 0 0 0 1 0 0 0
57 471 19.5 0 0 0 0 0 0 0 0 1 0 0
58 453 19.9 0 0 0 0 0 0 0 0 0 1 0
59 443 20.0 0 0 0 0 0 0 0 0 0 0 1
60 442 17.3 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) ExpBE M1 M2 M3 M4
109.356 20.139 1.695 -2.744 -49.862 -21.242
M5 M6 M7 M8 M9 M10
-23.853 -39.975 6.419 51.923 -28.223 -32.637
M11
-28.967
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-40.661 -14.594 4.461 13.559 31.817
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 109.356 32.210 3.395 0.00140 **
ExpBE 20.139 1.899 10.606 4.64e-14 ***
M1 1.695 14.237 0.119 0.90576
M2 -2.744 14.199 -0.193 0.84757
M3 -49.862 14.621 -3.410 0.00134 **
M4 -21.242 14.190 -1.497 0.14109
M5 -23.853 14.197 -1.680 0.09956 .
M6 -39.975 14.379 -2.780 0.00779 **
M7 6.419 14.181 0.453 0.65286
M8 51.923 14.847 3.497 0.00104 **
M9 -28.223 14.500 -1.946 0.05760 .
M10 -32.637 14.541 -2.244 0.02955 *
M11 -28.967 14.360 -2.017 0.04941 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 22.42 on 47 degrees of freedom
Multiple R-squared: 0.7188, Adjusted R-squared: 0.647
F-statistic: 10.01 on 12 and 47 DF, p-value: 2.535e-09
> 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.7359729 5.280541e-01 2.640271e-01
[2,] 0.9817544 3.649127e-02 1.824563e-02
[3,] 0.9793523 4.129536e-02 2.064768e-02
[4,] 0.9895442 2.091152e-02 1.045576e-02
[5,] 0.9969972 6.005685e-03 3.002842e-03
[6,] 0.9970597 5.880627e-03 2.940313e-03
[7,] 0.9987066 2.586762e-03 1.293381e-03
[8,] 0.9982922 3.415660e-03 1.707830e-03
[9,] 0.9970866 5.826741e-03 2.913370e-03
[10,] 0.9987745 2.450902e-03 1.225451e-03
[11,] 0.9988149 2.370257e-03 1.185128e-03
[12,] 0.9982290 3.542065e-03 1.771033e-03
[13,] 0.9987713 2.457489e-03 1.228744e-03
[14,] 0.9993688 1.262452e-03 6.312258e-04
[15,] 0.9997220 5.560457e-04 2.780228e-04
[16,] 0.9999531 9.384137e-05 4.692069e-05
[17,] 0.9999819 3.627589e-05 1.813795e-05
[18,] 0.9999909 1.810427e-05 9.052134e-06
[19,] 0.9999777 4.458095e-05 2.229047e-05
[20,] 0.9999304 1.392649e-04 6.963243e-05
[21,] 0.9998397 3.205156e-04 1.602578e-04
[22,] 0.9996190 7.620827e-04 3.810413e-04
[23,] 0.9988665 2.266975e-03 1.133487e-03
[24,] 0.9965751 6.849773e-03 3.424886e-03
[25,] 0.9907355 1.852904e-02 9.264519e-03
[26,] 0.9787315 4.253704e-02 2.126852e-02
[27,] 0.9517641 9.647176e-02 4.823588e-02
[28,] 0.9046001 1.907998e-01 9.539991e-02
[29,] 0.8099349 3.801302e-01 1.900651e-01
> postscript(file="/var/www/html/rcomp/tmp/1jpuq1258734605.ps",horizontal=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/rcomp/tmp/2u1ey1258734605.ps",horizontal=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/rcomp/tmp/3sz5m1258734605.ps",horizontal=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/rcomp/tmp/45k4l1258734605.ps",horizontal=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/rcomp/tmp/5e05n1258734605.ps",horizontal=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
-36.0413213 -28.5884284 -16.7076972 -38.2441052 -40.6608082 -16.4826426
7 8 9 10 11 12
-37.8635920 -26.9633739 -23.3327517 -37.0719435 -5.4909941 -1.2631558
13 14 15 16 17 18
-17.1109172 -8.6580243 10.2783837 -2.2580243 4.4087877 14.4616806
19 20 21 22 23 24
-0.8775112 24.0923028 10.5976524 4.9558948 27.5229250 8.5698140
25 26 27 28 29 30
20.9169211 13.2584606 10.0139192 5.6027838 27.3531110 -2.8723797
31 32 33 34 35 36
5.0111353 7.7860808 12.3888647 24.7888647 24.1888647 15.3192687
37 38 39 40 41 42
27.7220526 22.0775112 17.8468890 9.4079153 23.1025657 12.0301860
43 44 45 46 47 48
27.9554586 6.6190506 3.1939961 31.8167030 -6.0477615 -6.8616806
49 50 51 52 53 54
4.5132648 1.9104810 -21.4314947 25.4914304 -14.2036563 -7.1368442
55 56 57 58 59 60
5.7745092 -11.5340604 -2.8477615 -24.4895190 -40.1730341 -15.7642464
> postscript(file="/var/www/html/rcomp/tmp/6p5ky1258734605.ps",horizontal=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 -36.0413213 NA
1 -28.5884284 -36.0413213
2 -16.7076972 -28.5884284
3 -38.2441052 -16.7076972
4 -40.6608082 -38.2441052
5 -16.4826426 -40.6608082
6 -37.8635920 -16.4826426
7 -26.9633739 -37.8635920
8 -23.3327517 -26.9633739
9 -37.0719435 -23.3327517
10 -5.4909941 -37.0719435
11 -1.2631558 -5.4909941
12 -17.1109172 -1.2631558
13 -8.6580243 -17.1109172
14 10.2783837 -8.6580243
15 -2.2580243 10.2783837
16 4.4087877 -2.2580243
17 14.4616806 4.4087877
18 -0.8775112 14.4616806
19 24.0923028 -0.8775112
20 10.5976524 24.0923028
21 4.9558948 10.5976524
22 27.5229250 4.9558948
23 8.5698140 27.5229250
24 20.9169211 8.5698140
25 13.2584606 20.9169211
26 10.0139192 13.2584606
27 5.6027838 10.0139192
28 27.3531110 5.6027838
29 -2.8723797 27.3531110
30 5.0111353 -2.8723797
31 7.7860808 5.0111353
32 12.3888647 7.7860808
33 24.7888647 12.3888647
34 24.1888647 24.7888647
35 15.3192687 24.1888647
36 27.7220526 15.3192687
37 22.0775112 27.7220526
38 17.8468890 22.0775112
39 9.4079153 17.8468890
40 23.1025657 9.4079153
41 12.0301860 23.1025657
42 27.9554586 12.0301860
43 6.6190506 27.9554586
44 3.1939961 6.6190506
45 31.8167030 3.1939961
46 -6.0477615 31.8167030
47 -6.8616806 -6.0477615
48 4.5132648 -6.8616806
49 1.9104810 4.5132648
50 -21.4314947 1.9104810
51 25.4914304 -21.4314947
52 -14.2036563 25.4914304
53 -7.1368442 -14.2036563
54 5.7745092 -7.1368442
55 -11.5340604 5.7745092
56 -2.8477615 -11.5340604
57 -24.4895190 -2.8477615
58 -40.1730341 -24.4895190
59 -15.7642464 -40.1730341
60 NA -15.7642464
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -28.5884284 -36.0413213
[2,] -16.7076972 -28.5884284
[3,] -38.2441052 -16.7076972
[4,] -40.6608082 -38.2441052
[5,] -16.4826426 -40.6608082
[6,] -37.8635920 -16.4826426
[7,] -26.9633739 -37.8635920
[8,] -23.3327517 -26.9633739
[9,] -37.0719435 -23.3327517
[10,] -5.4909941 -37.0719435
[11,] -1.2631558 -5.4909941
[12,] -17.1109172 -1.2631558
[13,] -8.6580243 -17.1109172
[14,] 10.2783837 -8.6580243
[15,] -2.2580243 10.2783837
[16,] 4.4087877 -2.2580243
[17,] 14.4616806 4.4087877
[18,] -0.8775112 14.4616806
[19,] 24.0923028 -0.8775112
[20,] 10.5976524 24.0923028
[21,] 4.9558948 10.5976524
[22,] 27.5229250 4.9558948
[23,] 8.5698140 27.5229250
[24,] 20.9169211 8.5698140
[25,] 13.2584606 20.9169211
[26,] 10.0139192 13.2584606
[27,] 5.6027838 10.0139192
[28,] 27.3531110 5.6027838
[29,] -2.8723797 27.3531110
[30,] 5.0111353 -2.8723797
[31,] 7.7860808 5.0111353
[32,] 12.3888647 7.7860808
[33,] 24.7888647 12.3888647
[34,] 24.1888647 24.7888647
[35,] 15.3192687 24.1888647
[36,] 27.7220526 15.3192687
[37,] 22.0775112 27.7220526
[38,] 17.8468890 22.0775112
[39,] 9.4079153 17.8468890
[40,] 23.1025657 9.4079153
[41,] 12.0301860 23.1025657
[42,] 27.9554586 12.0301860
[43,] 6.6190506 27.9554586
[44,] 3.1939961 6.6190506
[45,] 31.8167030 3.1939961
[46,] -6.0477615 31.8167030
[47,] -6.8616806 -6.0477615
[48,] 4.5132648 -6.8616806
[49,] 1.9104810 4.5132648
[50,] -21.4314947 1.9104810
[51,] 25.4914304 -21.4314947
[52,] -14.2036563 25.4914304
[53,] -7.1368442 -14.2036563
[54,] 5.7745092 -7.1368442
[55,] -11.5340604 5.7745092
[56,] -2.8477615 -11.5340604
[57,] -24.4895190 -2.8477615
[58,] -40.1730341 -24.4895190
[59,] -15.7642464 -40.1730341
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -28.5884284 -36.0413213
2 -16.7076972 -28.5884284
3 -38.2441052 -16.7076972
4 -40.6608082 -38.2441052
5 -16.4826426 -40.6608082
6 -37.8635920 -16.4826426
7 -26.9633739 -37.8635920
8 -23.3327517 -26.9633739
9 -37.0719435 -23.3327517
10 -5.4909941 -37.0719435
11 -1.2631558 -5.4909941
12 -17.1109172 -1.2631558
13 -8.6580243 -17.1109172
14 10.2783837 -8.6580243
15 -2.2580243 10.2783837
16 4.4087877 -2.2580243
17 14.4616806 4.4087877
18 -0.8775112 14.4616806
19 24.0923028 -0.8775112
20 10.5976524 24.0923028
21 4.9558948 10.5976524
22 27.5229250 4.9558948
23 8.5698140 27.5229250
24 20.9169211 8.5698140
25 13.2584606 20.9169211
26 10.0139192 13.2584606
27 5.6027838 10.0139192
28 27.3531110 5.6027838
29 -2.8723797 27.3531110
30 5.0111353 -2.8723797
31 7.7860808 5.0111353
32 12.3888647 7.7860808
33 24.7888647 12.3888647
34 24.1888647 24.7888647
35 15.3192687 24.1888647
36 27.7220526 15.3192687
37 22.0775112 27.7220526
38 17.8468890 22.0775112
39 9.4079153 17.8468890
40 23.1025657 9.4079153
41 12.0301860 23.1025657
42 27.9554586 12.0301860
43 6.6190506 27.9554586
44 3.1939961 6.6190506
45 31.8167030 3.1939961
46 -6.0477615 31.8167030
47 -6.8616806 -6.0477615
48 4.5132648 -6.8616806
49 1.9104810 4.5132648
50 -21.4314947 1.9104810
51 25.4914304 -21.4314947
52 -14.2036563 25.4914304
53 -7.1368442 -14.2036563
54 5.7745092 -7.1368442
55 -11.5340604 5.7745092
56 -2.8477615 -11.5340604
57 -24.4895190 -2.8477615
58 -40.1730341 -24.4895190
59 -15.7642464 -40.1730341
> 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/rcomp/tmp/7sxcj1258734605.ps",horizontal=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/rcomp/tmp/833m31258734605.ps",horizontal=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/rcomp/tmp/90rma1258734605.ps",horizontal=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/rcomp/tmp/10vzbp1258734605.ps",horizontal=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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/rcomp/tmp/118ffp1258734605.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/rcomp/tmp/126cce1258734605.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/rcomp/tmp/13lh2w1258734605.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/rcomp/tmp/14ds9f1258734605.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/rcomp/tmp/15vvnb1258734605.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/rcomp/tmp/16ldzr1258734605.tab")
+ }
>
> system("convert tmp/1jpuq1258734605.ps tmp/1jpuq1258734605.png")
> system("convert tmp/2u1ey1258734605.ps tmp/2u1ey1258734605.png")
> system("convert tmp/3sz5m1258734605.ps tmp/3sz5m1258734605.png")
> system("convert tmp/45k4l1258734605.ps tmp/45k4l1258734605.png")
> system("convert tmp/5e05n1258734605.ps tmp/5e05n1258734605.png")
> system("convert tmp/6p5ky1258734605.ps tmp/6p5ky1258734605.png")
> system("convert tmp/7sxcj1258734605.ps tmp/7sxcj1258734605.png")
> system("convert tmp/833m31258734605.ps tmp/833m31258734605.png")
> system("convert tmp/90rma1258734605.ps tmp/90rma1258734605.png")
> system("convert tmp/10vzbp1258734605.ps tmp/10vzbp1258734605.png")
>
>
> proc.time()
user system elapsed
2.384 1.541 2.787