R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(493,797,514,840,522,988,490,819,484,831,506,904,501,814,462,798,465,828,454,789,464,930,427,744,460,832,473,826,465,907,422,776,415,835,413,715,420,729,363,733,376,736,380,712,384,711,346,667,389,799,407,661,393,692,346,649,348,729,353,622,364,671,305,635,307,648,312,745,312,624,286,477,324,710,336,515,327,461,302,590,299,415,311,554,315,585,264,513,278,591,278,561,287,684,279,668,324,795,354,776,354,1,043,360,964,363,762,385,1,030,412,939,370,779,389,918,395,839,417,874,404,840),dim=c(2,60),dimnames=list(c('WLH','Faill'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('WLH','Faill'),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
WLH Faill M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 493 797 1 0 0 0 0 0 0 0 0 0 0
2 514 840 0 1 0 0 0 0 0 0 0 0 0
3 522 988 0 0 1 0 0 0 0 0 0 0 0
4 490 819 0 0 0 1 0 0 0 0 0 0 0
5 484 831 0 0 0 0 1 0 0 0 0 0 0
6 506 904 0 0 0 0 0 1 0 0 0 0 0
7 501 814 0 0 0 0 0 0 1 0 0 0 0
8 462 798 0 0 0 0 0 0 0 1 0 0 0
9 465 828 0 0 0 0 0 0 0 0 1 0 0
10 454 789 0 0 0 0 0 0 0 0 0 1 0
11 464 930 0 0 0 0 0 0 0 0 0 0 1
12 427 744 0 0 0 0 0 0 0 0 0 0 0
13 460 832 1 0 0 0 0 0 0 0 0 0 0
14 473 826 0 1 0 0 0 0 0 0 0 0 0
15 465 907 0 0 1 0 0 0 0 0 0 0 0
16 422 776 0 0 0 1 0 0 0 0 0 0 0
17 415 835 0 0 0 0 1 0 0 0 0 0 0
18 413 715 0 0 0 0 0 1 0 0 0 0 0
19 420 729 0 0 0 0 0 0 1 0 0 0 0
20 363 733 0 0 0 0 0 0 0 1 0 0 0
21 376 736 0 0 0 0 0 0 0 0 1 0 0
22 380 712 0 0 0 0 0 0 0 0 0 1 0
23 384 711 0 0 0 0 0 0 0 0 0 0 1
24 346 667 0 0 0 0 0 0 0 0 0 0 0
25 389 799 1 0 0 0 0 0 0 0 0 0 0
26 407 661 0 1 0 0 0 0 0 0 0 0 0
27 393 692 0 0 1 0 0 0 0 0 0 0 0
28 346 649 0 0 0 1 0 0 0 0 0 0 0
29 348 729 0 0 0 0 1 0 0 0 0 0 0
30 353 622 0 0 0 0 0 1 0 0 0 0 0
31 364 671 0 0 0 0 0 0 1 0 0 0 0
32 305 635 0 0 0 0 0 0 0 1 0 0 0
33 307 648 0 0 0 0 0 0 0 0 1 0 0
34 312 745 0 0 0 0 0 0 0 0 0 1 0
35 312 624 0 0 0 0 0 0 0 0 0 0 1
36 286 477 0 0 0 0 0 0 0 0 0 0 0
37 324 710 1 0 0 0 0 0 0 0 0 0 0
38 336 515 0 1 0 0 0 0 0 0 0 0 0
39 327 461 0 0 1 0 0 0 0 0 0 0 0
40 302 590 0 0 0 1 0 0 0 0 0 0 0
41 299 415 0 0 0 0 1 0 0 0 0 0 0
42 311 554 0 0 0 0 0 1 0 0 0 0 0
43 315 585 0 0 0 0 0 0 1 0 0 0 0
44 264 513 0 0 0 0 0 0 0 1 0 0 0
45 278 591 0 0 0 0 0 0 0 0 1 0 0
46 278 561 0 0 0 0 0 0 0 0 0 1 0
47 287 684 0 0 0 0 0 0 0 0 0 0 1
48 279 668 0 0 0 0 0 0 0 0 0 0 0
49 324 795 1 0 0 0 0 0 0 0 0 0 0
50 354 776 0 1 0 0 0 0 0 0 0 0 0
51 354 1 0 0 1 0 0 0 0 0 0 0 0
52 43 360 0 0 0 1 0 0 0 0 0 0 0
53 964 363 0 0 0 0 1 0 0 0 0 0 0
54 762 385 0 0 0 0 0 1 0 0 0 0 0
55 1 30 0 0 0 0 0 0 1 0 0 0 0
56 412 939 0 0 0 0 0 0 0 1 0 0 0
57 370 779 0 0 0 0 0 0 0 0 1 0 0
58 389 918 0 0 0 0 0 0 0 0 0 1 0
59 395 839 0 0 0 0 0 0 0 0 0 0 1
60 417 874 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) Faill M1 M2 M3 M4
177.6398 0.2527 21.5772 56.2980 80.4566 -18.4720
M5 M6 M7 M8 M9 M10
163.9894 130.6356 -0.4241 0.6980 0.5176 -3.3100
M11
-0.6942
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-207.14 -49.78 -8.16 35.37 530.64
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 177.63976 83.00979 2.140 0.03758 *
Faill 0.25271 0.08956 2.822 0.00698 **
M1 21.57720 79.45439 0.272 0.78714
M2 56.29804 79.01369 0.713 0.47967
M3 80.45663 79.23634 1.015 0.31511
M4 -18.47201 79.05500 -0.234 0.81626
M5 163.98938 79.07601 2.074 0.04360 *
M6 130.63559 79.06881 1.652 0.10516
M7 -0.42405 79.67253 -0.005 0.99578
M8 0.69804 79.01369 0.009 0.99299
M9 0.51756 78.98883 0.007 0.99480
M10 -3.30999 79.11855 -0.042 0.96681
M11 -0.69416 79.20192 -0.009 0.99304
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 124.8 on 47 degrees of freedom
Multiple R-squared: 0.287, Adjusted R-squared: 0.105
F-statistic: 1.577 on 12 and 47 DF, p-value: 0.1316
> 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,] 1.458573e-02 2.917147e-02 0.9854143
[2,] 9.473738e-03 1.894748e-02 0.9905263
[3,] 2.260547e-03 4.521094e-03 0.9977395
[4,] 6.412733e-04 1.282547e-03 0.9993587
[5,] 2.857484e-04 5.714969e-04 0.9997143
[6,] 6.960476e-05 1.392095e-04 0.9999304
[7,] 1.515221e-05 3.030442e-05 0.9999848
[8,] 7.199677e-06 1.439935e-05 0.9999928
[9,] 1.874062e-06 3.748125e-06 0.9999981
[10,] 2.144981e-06 4.289963e-06 0.9999979
[11,] 4.789318e-07 9.578637e-07 0.9999995
[12,] 1.449290e-07 2.898580e-07 0.9999999
[13,] 4.283996e-08 8.567991e-08 1.0000000
[14,] 4.849387e-08 9.698774e-08 1.0000000
[15,] 1.748617e-08 3.497235e-08 1.0000000
[16,] 6.116741e-09 1.223348e-08 1.0000000
[17,] 1.519659e-09 3.039318e-09 1.0000000
[18,] 3.772644e-10 7.545289e-10 1.0000000
[19,] 1.069614e-09 2.139227e-09 1.0000000
[20,] 2.002474e-10 4.004947e-10 1.0000000
[21,] 4.762594e-11 9.525187e-11 1.0000000
[22,] 2.459286e-11 4.918571e-11 1.0000000
[23,] 4.891527e-12 9.783055e-12 1.0000000
[24,] 7.622631e-12 1.524526e-11 1.0000000
[25,] 2.601087e-12 5.202174e-12 1.0000000
[26,] 9.627130e-08 1.925426e-07 0.9999999
[27,] 5.650318e-01 8.699364e-01 0.4349682
[28,] 6.685220e-01 6.629560e-01 0.3314780
[29,] 5.557998e-01 8.884004e-01 0.4442002
> postscript(file="/var/www/html/freestat/rcomp/tmp/1bi8x1293456238.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/2bi8x1293456238.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3y70b1293456238.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4y70b1293456238.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5y70b1293456238.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
92.3717982 67.7843564 14.2244296 123.8613490 -67.6325807 -30.7267395
7 8 9 10 11 12
118.0769525 81.9982484 77.5973727 80.2806846 52.0325004 61.3427206
13 14 15 16 17 18
50.5268882 30.3223204 -22.3059216 66.7279528 -137.6434276 -75.9642254
19 20 21 22 23 24
58.5574482 -0.5754901 11.8468504 25.7394866 27.3763659 -0.1984774
25 26 27 28 29 30
-12.1336253 6.0197533 -39.9729029 22.8223405 -177.8559859 -112.4620360
31 32 33 34 35 36
17.2147276 -33.8097421 -34.9145187 -50.6000000 -22.6377149 -12.1832517
37 38 39 40 41 42
-54.6422827 -28.0843364 -47.5964969 -6.2676683 -147.5045075 -137.2776394
43 44 45 46 47 48
-10.0520649 -43.9789129 -49.5099510 -38.1010445 -62.8004178 -67.4511891
49 50 51 52 53 54
-76.1227784 -76.0420938 95.6508918 -207.1439740 530.6365017 356.4306403
55 56 57 58 59 60
-183.7970634 -3.6341033 -5.0197533 -17.3191266 6.0292664 18.4901977
> postscript(file="/var/www/html/freestat/rcomp/tmp/6rgze1293456238.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 92.3717982 NA
1 67.7843564 92.3717982
2 14.2244296 67.7843564
3 123.8613490 14.2244296
4 -67.6325807 123.8613490
5 -30.7267395 -67.6325807
6 118.0769525 -30.7267395
7 81.9982484 118.0769525
8 77.5973727 81.9982484
9 80.2806846 77.5973727
10 52.0325004 80.2806846
11 61.3427206 52.0325004
12 50.5268882 61.3427206
13 30.3223204 50.5268882
14 -22.3059216 30.3223204
15 66.7279528 -22.3059216
16 -137.6434276 66.7279528
17 -75.9642254 -137.6434276
18 58.5574482 -75.9642254
19 -0.5754901 58.5574482
20 11.8468504 -0.5754901
21 25.7394866 11.8468504
22 27.3763659 25.7394866
23 -0.1984774 27.3763659
24 -12.1336253 -0.1984774
25 6.0197533 -12.1336253
26 -39.9729029 6.0197533
27 22.8223405 -39.9729029
28 -177.8559859 22.8223405
29 -112.4620360 -177.8559859
30 17.2147276 -112.4620360
31 -33.8097421 17.2147276
32 -34.9145187 -33.8097421
33 -50.6000000 -34.9145187
34 -22.6377149 -50.6000000
35 -12.1832517 -22.6377149
36 -54.6422827 -12.1832517
37 -28.0843364 -54.6422827
38 -47.5964969 -28.0843364
39 -6.2676683 -47.5964969
40 -147.5045075 -6.2676683
41 -137.2776394 -147.5045075
42 -10.0520649 -137.2776394
43 -43.9789129 -10.0520649
44 -49.5099510 -43.9789129
45 -38.1010445 -49.5099510
46 -62.8004178 -38.1010445
47 -67.4511891 -62.8004178
48 -76.1227784 -67.4511891
49 -76.0420938 -76.1227784
50 95.6508918 -76.0420938
51 -207.1439740 95.6508918
52 530.6365017 -207.1439740
53 356.4306403 530.6365017
54 -183.7970634 356.4306403
55 -3.6341033 -183.7970634
56 -5.0197533 -3.6341033
57 -17.3191266 -5.0197533
58 6.0292664 -17.3191266
59 18.4901977 6.0292664
60 NA 18.4901977
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 67.7843564 92.3717982
[2,] 14.2244296 67.7843564
[3,] 123.8613490 14.2244296
[4,] -67.6325807 123.8613490
[5,] -30.7267395 -67.6325807
[6,] 118.0769525 -30.7267395
[7,] 81.9982484 118.0769525
[8,] 77.5973727 81.9982484
[9,] 80.2806846 77.5973727
[10,] 52.0325004 80.2806846
[11,] 61.3427206 52.0325004
[12,] 50.5268882 61.3427206
[13,] 30.3223204 50.5268882
[14,] -22.3059216 30.3223204
[15,] 66.7279528 -22.3059216
[16,] -137.6434276 66.7279528
[17,] -75.9642254 -137.6434276
[18,] 58.5574482 -75.9642254
[19,] -0.5754901 58.5574482
[20,] 11.8468504 -0.5754901
[21,] 25.7394866 11.8468504
[22,] 27.3763659 25.7394866
[23,] -0.1984774 27.3763659
[24,] -12.1336253 -0.1984774
[25,] 6.0197533 -12.1336253
[26,] -39.9729029 6.0197533
[27,] 22.8223405 -39.9729029
[28,] -177.8559859 22.8223405
[29,] -112.4620360 -177.8559859
[30,] 17.2147276 -112.4620360
[31,] -33.8097421 17.2147276
[32,] -34.9145187 -33.8097421
[33,] -50.6000000 -34.9145187
[34,] -22.6377149 -50.6000000
[35,] -12.1832517 -22.6377149
[36,] -54.6422827 -12.1832517
[37,] -28.0843364 -54.6422827
[38,] -47.5964969 -28.0843364
[39,] -6.2676683 -47.5964969
[40,] -147.5045075 -6.2676683
[41,] -137.2776394 -147.5045075
[42,] -10.0520649 -137.2776394
[43,] -43.9789129 -10.0520649
[44,] -49.5099510 -43.9789129
[45,] -38.1010445 -49.5099510
[46,] -62.8004178 -38.1010445
[47,] -67.4511891 -62.8004178
[48,] -76.1227784 -67.4511891
[49,] -76.0420938 -76.1227784
[50,] 95.6508918 -76.0420938
[51,] -207.1439740 95.6508918
[52,] 530.6365017 -207.1439740
[53,] 356.4306403 530.6365017
[54,] -183.7970634 356.4306403
[55,] -3.6341033 -183.7970634
[56,] -5.0197533 -3.6341033
[57,] -17.3191266 -5.0197533
[58,] 6.0292664 -17.3191266
[59,] 18.4901977 6.0292664
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 67.7843564 92.3717982
2 14.2244296 67.7843564
3 123.8613490 14.2244296
4 -67.6325807 123.8613490
5 -30.7267395 -67.6325807
6 118.0769525 -30.7267395
7 81.9982484 118.0769525
8 77.5973727 81.9982484
9 80.2806846 77.5973727
10 52.0325004 80.2806846
11 61.3427206 52.0325004
12 50.5268882 61.3427206
13 30.3223204 50.5268882
14 -22.3059216 30.3223204
15 66.7279528 -22.3059216
16 -137.6434276 66.7279528
17 -75.9642254 -137.6434276
18 58.5574482 -75.9642254
19 -0.5754901 58.5574482
20 11.8468504 -0.5754901
21 25.7394866 11.8468504
22 27.3763659 25.7394866
23 -0.1984774 27.3763659
24 -12.1336253 -0.1984774
25 6.0197533 -12.1336253
26 -39.9729029 6.0197533
27 22.8223405 -39.9729029
28 -177.8559859 22.8223405
29 -112.4620360 -177.8559859
30 17.2147276 -112.4620360
31 -33.8097421 17.2147276
32 -34.9145187 -33.8097421
33 -50.6000000 -34.9145187
34 -22.6377149 -50.6000000
35 -12.1832517 -22.6377149
36 -54.6422827 -12.1832517
37 -28.0843364 -54.6422827
38 -47.5964969 -28.0843364
39 -6.2676683 -47.5964969
40 -147.5045075 -6.2676683
41 -137.2776394 -147.5045075
42 -10.0520649 -137.2776394
43 -43.9789129 -10.0520649
44 -49.5099510 -43.9789129
45 -38.1010445 -49.5099510
46 -62.8004178 -38.1010445
47 -67.4511891 -62.8004178
48 -76.1227784 -67.4511891
49 -76.0420938 -76.1227784
50 95.6508918 -76.0420938
51 -207.1439740 95.6508918
52 530.6365017 -207.1439740
53 356.4306403 530.6365017
54 -183.7970634 356.4306403
55 -3.6341033 -183.7970634
56 -5.0197533 -3.6341033
57 -17.3191266 -5.0197533
58 6.0292664 -17.3191266
59 18.4901977 6.0292664
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/72pyh1293456238.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/82pyh1293456238.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/92pyh1293456238.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10cgfk1293456238.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/11yheq1293456238.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/121hde1293456238.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/13f9a41293456238.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/141s9s1293456238.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/15t1qd1293456238.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/16psom1293456238.tab")
+ }
> try(system("convert tmp/1bi8x1293456238.ps tmp/1bi8x1293456238.png",intern=TRUE))
character(0)
> try(system("convert tmp/2bi8x1293456238.ps tmp/2bi8x1293456238.png",intern=TRUE))
character(0)
> try(system("convert tmp/3y70b1293456238.ps tmp/3y70b1293456238.png",intern=TRUE))
character(0)
> try(system("convert tmp/4y70b1293456238.ps tmp/4y70b1293456238.png",intern=TRUE))
character(0)
> try(system("convert tmp/5y70b1293456238.ps tmp/5y70b1293456238.png",intern=TRUE))
character(0)
> try(system("convert tmp/6rgze1293456238.ps tmp/6rgze1293456238.png",intern=TRUE))
character(0)
> try(system("convert tmp/72pyh1293456238.ps tmp/72pyh1293456238.png",intern=TRUE))
character(0)
> try(system("convert tmp/82pyh1293456238.ps tmp/82pyh1293456238.png",intern=TRUE))
character(0)
> try(system("convert tmp/92pyh1293456238.ps tmp/92pyh1293456238.png",intern=TRUE))
character(0)
> try(system("convert tmp/10cgfk1293456238.ps tmp/10cgfk1293456238.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.754 2.424 4.040