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(613,0,611,0,594,0,595,0,591,0,589,0,584,0,573,0,567,0,569,0,621,0,629,0,628,0,612,0,595,0,597,0,593,0,590,0,580,0,574,0,573,0,573,0,620,0,626,0,620,0,588,0,566,0,557,0,561,0,549,0,532,0,526,0,511,0,499,0,555,0,565,0,542,0,527,0,510,0,514,0,517,0,508,0,493,0,490,0,469,0,478,0,528,0,534,0,518,1,506,1,502,1,516,1,528,1,533,1,536,1,537,1,524,1,536,1,587,1,597,1,581,1),dim=c(2,61),dimnames=list(c('WklBe','X'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('WklBe','X'),1:61))
> 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
WklBe X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 613 0 1 0 0 0 0 0 0 0 0 0 0
2 611 0 0 1 0 0 0 0 0 0 0 0 0
3 594 0 0 0 1 0 0 0 0 0 0 0 0
4 595 0 0 0 0 1 0 0 0 0 0 0 0
5 591 0 0 0 0 0 1 0 0 0 0 0 0
6 589 0 0 0 0 0 0 1 0 0 0 0 0
7 584 0 0 0 0 0 0 0 1 0 0 0 0
8 573 0 0 0 0 0 0 0 0 1 0 0 0
9 567 0 0 0 0 0 0 0 0 0 1 0 0
10 569 0 0 0 0 0 0 0 0 0 0 1 0
11 621 0 0 0 0 0 0 0 0 0 0 0 1
12 629 0 0 0 0 0 0 0 0 0 0 0 0
13 628 0 1 0 0 0 0 0 0 0 0 0 0
14 612 0 0 1 0 0 0 0 0 0 0 0 0
15 595 0 0 0 1 0 0 0 0 0 0 0 0
16 597 0 0 0 0 1 0 0 0 0 0 0 0
17 593 0 0 0 0 0 1 0 0 0 0 0 0
18 590 0 0 0 0 0 0 1 0 0 0 0 0
19 580 0 0 0 0 0 0 0 1 0 0 0 0
20 574 0 0 0 0 0 0 0 0 1 0 0 0
21 573 0 0 0 0 0 0 0 0 0 1 0 0
22 573 0 0 0 0 0 0 0 0 0 0 1 0
23 620 0 0 0 0 0 0 0 0 0 0 0 1
24 626 0 0 0 0 0 0 0 0 0 0 0 0
25 620 0 1 0 0 0 0 0 0 0 0 0 0
26 588 0 0 1 0 0 0 0 0 0 0 0 0
27 566 0 0 0 1 0 0 0 0 0 0 0 0
28 557 0 0 0 0 1 0 0 0 0 0 0 0
29 561 0 0 0 0 0 1 0 0 0 0 0 0
30 549 0 0 0 0 0 0 1 0 0 0 0 0
31 532 0 0 0 0 0 0 0 1 0 0 0 0
32 526 0 0 0 0 0 0 0 0 1 0 0 0
33 511 0 0 0 0 0 0 0 0 0 1 0 0
34 499 0 0 0 0 0 0 0 0 0 0 1 0
35 555 0 0 0 0 0 0 0 0 0 0 0 1
36 565 0 0 0 0 0 0 0 0 0 0 0 0
37 542 0 1 0 0 0 0 0 0 0 0 0 0
38 527 0 0 1 0 0 0 0 0 0 0 0 0
39 510 0 0 0 1 0 0 0 0 0 0 0 0
40 514 0 0 0 0 1 0 0 0 0 0 0 0
41 517 0 0 0 0 0 1 0 0 0 0 0 0
42 508 0 0 0 0 0 0 1 0 0 0 0 0
43 493 0 0 0 0 0 0 0 1 0 0 0 0
44 490 0 0 0 0 0 0 0 0 1 0 0 0
45 469 0 0 0 0 0 0 0 0 0 1 0 0
46 478 0 0 0 0 0 0 0 0 0 0 1 0
47 528 0 0 0 0 0 0 0 0 0 0 0 1
48 534 0 0 0 0 0 0 0 0 0 0 0 0
49 518 1 1 0 0 0 0 0 0 0 0 0 0
50 506 1 0 1 0 0 0 0 0 0 0 0 0
51 502 1 0 0 1 0 0 0 0 0 0 0 0
52 516 1 0 0 0 1 0 0 0 0 0 0 0
53 528 1 0 0 0 0 1 0 0 0 0 0 0
54 533 1 0 0 0 0 0 1 0 0 0 0 0
55 536 1 0 0 0 0 0 0 1 0 0 0 0
56 537 1 0 0 0 0 0 0 0 1 0 0 0
57 524 1 0 0 0 0 0 0 0 0 1 0 0
58 536 1 0 0 0 0 0 0 0 0 0 1 0
59 587 1 0 0 0 0 0 0 0 0 0 0 1
60 597 1 0 0 0 0 0 0 0 0 0 0 0
61 581 1 1 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) X M1 M2 M3 M4
595.595 -26.974 -2.937 -21.400 -36.800 -34.400
M5 M6 M7 M8 M9 M10
-32.200 -36.400 -45.200 -50.200 -61.400 -59.200
M11
-8.000
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-65.19 -32.59 15.32 30.81 38.81
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 595.595 17.908 33.258 <2e-16 ***
X -26.974 12.457 -2.165 0.0354 *
M1 -2.937 24.069 -0.122 0.9034
M2 -21.400 25.080 -0.853 0.3977
M3 -36.800 25.080 -1.467 0.1488
M4 -34.400 25.080 -1.372 0.1766
M5 -32.200 25.080 -1.284 0.2053
M6 -36.400 25.080 -1.451 0.1532
M7 -45.200 25.080 -1.802 0.0778 .
M8 -50.200 25.080 -2.002 0.0510 .
M9 -61.400 25.080 -2.448 0.0181 *
M10 -59.200 25.080 -2.360 0.0224 *
M11 -8.000 25.080 -0.319 0.7511
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 39.65 on 48 degrees of freedom
Multiple R-squared: 0.2885, Adjusted R-squared: 0.1107
F-statistic: 1.622 on 12 and 48 DF, p-value: 0.1170
> 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,] 5.796367e-03 1.159273e-02 0.99420363
[2,] 8.448656e-04 1.689731e-03 0.99915513
[3,] 1.158555e-04 2.317111e-04 0.99988414
[4,] 1.854712e-05 3.709424e-05 0.99998145
[5,] 2.470792e-06 4.941583e-06 0.99999753
[6,] 6.482995e-07 1.296599e-06 0.99999935
[7,] 1.436044e-07 2.872088e-07 0.99999986
[8,] 2.674861e-08 5.349721e-08 0.99999997
[9,] 6.140668e-09 1.228134e-08 0.99999999
[10,] 2.266691e-09 4.533382e-09 1.00000000
[11,] 1.901421e-06 3.802842e-06 0.99999810
[12,] 1.254156e-04 2.508311e-04 0.99987458
[13,] 4.366502e-03 8.733004e-03 0.99563350
[14,] 1.686519e-02 3.373037e-02 0.98313481
[15,] 6.644272e-02 1.328854e-01 0.93355728
[16,] 1.895690e-01 3.791381e-01 0.81043095
[17,] 3.023486e-01 6.046973e-01 0.69765136
[18,] 4.865717e-01 9.731434e-01 0.51342828
[19,] 6.402811e-01 7.194377e-01 0.35971885
[20,] 6.999020e-01 6.001961e-01 0.30009805
[21,] 7.224876e-01 5.550248e-01 0.27751240
[22,] 7.938850e-01 4.122300e-01 0.20611500
[23,] 8.984568e-01 2.030865e-01 0.10154325
[24,] 9.409544e-01 1.180912e-01 0.05904562
[25,] 9.577526e-01 8.449481e-02 0.04224740
[26,] 9.637622e-01 7.247556e-02 0.03623778
[27,] 9.564565e-01 8.708701e-02 0.04354350
[28,] 9.255026e-01 1.489949e-01 0.07449745
[29,] 8.647219e-01 2.705562e-01 0.13527810
[30,] 7.583393e-01 4.833214e-01 0.24166068
> postscript(file="/var/www/html/rcomp/tmp/1tm6p1258726007.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/2hd4j1258726007.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/3d8bz1258726007.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/45gcw1258726007.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/59vuw1258726007.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 = 61
Frequency = 1
1 2 3 4 5 6
20.3421053 36.8052632 35.2052632 33.8052632 27.6052632 29.8052632
7 8 9 10 11 12
33.6052632 27.6052632 32.8052632 32.6052632 33.4052632 33.4052632
13 14 15 16 17 18
35.3421053 37.8052632 36.2052632 35.8052632 29.6052632 30.8052632
19 20 21 22 23 24
29.6052632 28.6052632 38.8052632 36.6052632 32.4052632 30.4052632
25 26 27 28 29 30
27.3421053 13.8052632 7.2052632 -4.1947368 -2.3947368 -10.1947368
31 32 33 34 35 36
-18.3947368 -19.3947368 -23.1947368 -37.3947368 -32.5947368 -30.5947368
37 38 39 40 41 42
-50.6578947 -47.1947368 -48.7947368 -47.1947368 -46.3947368 -51.1947368
43 44 45 46 47 48
-57.3947368 -55.3947368 -65.1947368 -58.3947368 -59.5947368 -61.5947368
49 50 51 52 53 54
-47.6842105 -41.2210526 -29.8210526 -18.2210526 -8.4210526 0.7789474
55 56 57 58 59 60
12.5789474 18.5789474 16.7789474 26.5789474 26.3789474 28.3789474
61
15.3157895
> postscript(file="/var/www/html/rcomp/tmp/61ky61258726007.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 20.3421053 NA
1 36.8052632 20.3421053
2 35.2052632 36.8052632
3 33.8052632 35.2052632
4 27.6052632 33.8052632
5 29.8052632 27.6052632
6 33.6052632 29.8052632
7 27.6052632 33.6052632
8 32.8052632 27.6052632
9 32.6052632 32.8052632
10 33.4052632 32.6052632
11 33.4052632 33.4052632
12 35.3421053 33.4052632
13 37.8052632 35.3421053
14 36.2052632 37.8052632
15 35.8052632 36.2052632
16 29.6052632 35.8052632
17 30.8052632 29.6052632
18 29.6052632 30.8052632
19 28.6052632 29.6052632
20 38.8052632 28.6052632
21 36.6052632 38.8052632
22 32.4052632 36.6052632
23 30.4052632 32.4052632
24 27.3421053 30.4052632
25 13.8052632 27.3421053
26 7.2052632 13.8052632
27 -4.1947368 7.2052632
28 -2.3947368 -4.1947368
29 -10.1947368 -2.3947368
30 -18.3947368 -10.1947368
31 -19.3947368 -18.3947368
32 -23.1947368 -19.3947368
33 -37.3947368 -23.1947368
34 -32.5947368 -37.3947368
35 -30.5947368 -32.5947368
36 -50.6578947 -30.5947368
37 -47.1947368 -50.6578947
38 -48.7947368 -47.1947368
39 -47.1947368 -48.7947368
40 -46.3947368 -47.1947368
41 -51.1947368 -46.3947368
42 -57.3947368 -51.1947368
43 -55.3947368 -57.3947368
44 -65.1947368 -55.3947368
45 -58.3947368 -65.1947368
46 -59.5947368 -58.3947368
47 -61.5947368 -59.5947368
48 -47.6842105 -61.5947368
49 -41.2210526 -47.6842105
50 -29.8210526 -41.2210526
51 -18.2210526 -29.8210526
52 -8.4210526 -18.2210526
53 0.7789474 -8.4210526
54 12.5789474 0.7789474
55 18.5789474 12.5789474
56 16.7789474 18.5789474
57 26.5789474 16.7789474
58 26.3789474 26.5789474
59 28.3789474 26.3789474
60 15.3157895 28.3789474
61 NA 15.3157895
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 36.8052632 20.3421053
[2,] 35.2052632 36.8052632
[3,] 33.8052632 35.2052632
[4,] 27.6052632 33.8052632
[5,] 29.8052632 27.6052632
[6,] 33.6052632 29.8052632
[7,] 27.6052632 33.6052632
[8,] 32.8052632 27.6052632
[9,] 32.6052632 32.8052632
[10,] 33.4052632 32.6052632
[11,] 33.4052632 33.4052632
[12,] 35.3421053 33.4052632
[13,] 37.8052632 35.3421053
[14,] 36.2052632 37.8052632
[15,] 35.8052632 36.2052632
[16,] 29.6052632 35.8052632
[17,] 30.8052632 29.6052632
[18,] 29.6052632 30.8052632
[19,] 28.6052632 29.6052632
[20,] 38.8052632 28.6052632
[21,] 36.6052632 38.8052632
[22,] 32.4052632 36.6052632
[23,] 30.4052632 32.4052632
[24,] 27.3421053 30.4052632
[25,] 13.8052632 27.3421053
[26,] 7.2052632 13.8052632
[27,] -4.1947368 7.2052632
[28,] -2.3947368 -4.1947368
[29,] -10.1947368 -2.3947368
[30,] -18.3947368 -10.1947368
[31,] -19.3947368 -18.3947368
[32,] -23.1947368 -19.3947368
[33,] -37.3947368 -23.1947368
[34,] -32.5947368 -37.3947368
[35,] -30.5947368 -32.5947368
[36,] -50.6578947 -30.5947368
[37,] -47.1947368 -50.6578947
[38,] -48.7947368 -47.1947368
[39,] -47.1947368 -48.7947368
[40,] -46.3947368 -47.1947368
[41,] -51.1947368 -46.3947368
[42,] -57.3947368 -51.1947368
[43,] -55.3947368 -57.3947368
[44,] -65.1947368 -55.3947368
[45,] -58.3947368 -65.1947368
[46,] -59.5947368 -58.3947368
[47,] -61.5947368 -59.5947368
[48,] -47.6842105 -61.5947368
[49,] -41.2210526 -47.6842105
[50,] -29.8210526 -41.2210526
[51,] -18.2210526 -29.8210526
[52,] -8.4210526 -18.2210526
[53,] 0.7789474 -8.4210526
[54,] 12.5789474 0.7789474
[55,] 18.5789474 12.5789474
[56,] 16.7789474 18.5789474
[57,] 26.5789474 16.7789474
[58,] 26.3789474 26.5789474
[59,] 28.3789474 26.3789474
[60,] 15.3157895 28.3789474
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 36.8052632 20.3421053
2 35.2052632 36.8052632
3 33.8052632 35.2052632
4 27.6052632 33.8052632
5 29.8052632 27.6052632
6 33.6052632 29.8052632
7 27.6052632 33.6052632
8 32.8052632 27.6052632
9 32.6052632 32.8052632
10 33.4052632 32.6052632
11 33.4052632 33.4052632
12 35.3421053 33.4052632
13 37.8052632 35.3421053
14 36.2052632 37.8052632
15 35.8052632 36.2052632
16 29.6052632 35.8052632
17 30.8052632 29.6052632
18 29.6052632 30.8052632
19 28.6052632 29.6052632
20 38.8052632 28.6052632
21 36.6052632 38.8052632
22 32.4052632 36.6052632
23 30.4052632 32.4052632
24 27.3421053 30.4052632
25 13.8052632 27.3421053
26 7.2052632 13.8052632
27 -4.1947368 7.2052632
28 -2.3947368 -4.1947368
29 -10.1947368 -2.3947368
30 -18.3947368 -10.1947368
31 -19.3947368 -18.3947368
32 -23.1947368 -19.3947368
33 -37.3947368 -23.1947368
34 -32.5947368 -37.3947368
35 -30.5947368 -32.5947368
36 -50.6578947 -30.5947368
37 -47.1947368 -50.6578947
38 -48.7947368 -47.1947368
39 -47.1947368 -48.7947368
40 -46.3947368 -47.1947368
41 -51.1947368 -46.3947368
42 -57.3947368 -51.1947368
43 -55.3947368 -57.3947368
44 -65.1947368 -55.3947368
45 -58.3947368 -65.1947368
46 -59.5947368 -58.3947368
47 -61.5947368 -59.5947368
48 -47.6842105 -61.5947368
49 -41.2210526 -47.6842105
50 -29.8210526 -41.2210526
51 -18.2210526 -29.8210526
52 -8.4210526 -18.2210526
53 0.7789474 -8.4210526
54 12.5789474 0.7789474
55 18.5789474 12.5789474
56 16.7789474 18.5789474
57 26.5789474 16.7789474
58 26.3789474 26.5789474
59 28.3789474 26.3789474
60 15.3157895 28.3789474
> 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/7lyfa1258726007.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/8mjp61258726007.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/994791258726007.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/10pr4v1258726007.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/1156on1258726007.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/12dsog1258726007.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/13pz6a1258726007.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/14mrra1258726007.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/15r0ox1258726007.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/162lb31258726007.tab")
+ }
>
> system("convert tmp/1tm6p1258726007.ps tmp/1tm6p1258726007.png")
> system("convert tmp/2hd4j1258726007.ps tmp/2hd4j1258726007.png")
> system("convert tmp/3d8bz1258726007.ps tmp/3d8bz1258726007.png")
> system("convert tmp/45gcw1258726007.ps tmp/45gcw1258726007.png")
> system("convert tmp/59vuw1258726007.ps tmp/59vuw1258726007.png")
> system("convert tmp/61ky61258726007.ps tmp/61ky61258726007.png")
> system("convert tmp/7lyfa1258726007.ps tmp/7lyfa1258726007.png")
> system("convert tmp/8mjp61258726007.ps tmp/8mjp61258726007.png")
> system("convert tmp/994791258726007.ps tmp/994791258726007.png")
> system("convert tmp/10pr4v1258726007.ps tmp/10pr4v1258726007.png")
>
>
> proc.time()
user system elapsed
2.364 1.619 2.869