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(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,1,534,1,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,564,1,564,1),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),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 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 595 0 1 0 0 0 0 0 0 0 0 0 0 1
2 591 0 0 1 0 0 0 0 0 0 0 0 0 2
3 589 0 0 0 1 0 0 0 0 0 0 0 0 3
4 584 0 0 0 0 1 0 0 0 0 0 0 0 4
5 573 0 0 0 0 0 1 0 0 0 0 0 0 5
6 567 0 0 0 0 0 0 1 0 0 0 0 0 6
7 569 0 0 0 0 0 0 0 1 0 0 0 0 7
8 621 0 0 0 0 0 0 0 0 1 0 0 0 8
9 629 0 0 0 0 0 0 0 0 0 1 0 0 9
10 628 0 0 0 0 0 0 0 0 0 0 1 0 10
11 612 0 0 0 0 0 0 0 0 0 0 0 1 11
12 595 0 0 0 0 0 0 0 0 0 0 0 0 12
13 597 0 1 0 0 0 0 0 0 0 0 0 0 13
14 593 0 0 1 0 0 0 0 0 0 0 0 0 14
15 590 0 0 0 1 0 0 0 0 0 0 0 0 15
16 580 0 0 0 0 1 0 0 0 0 0 0 0 16
17 574 0 0 0 0 0 1 0 0 0 0 0 0 17
18 573 0 0 0 0 0 0 1 0 0 0 0 0 18
19 573 0 0 0 0 0 0 0 1 0 0 0 0 19
20 620 0 0 0 0 0 0 0 0 1 0 0 0 20
21 626 0 0 0 0 0 0 0 0 0 1 0 0 21
22 620 0 0 0 0 0 0 0 0 0 0 1 0 22
23 588 0 0 0 0 0 0 0 0 0 0 0 1 23
24 566 0 0 0 0 0 0 0 0 0 0 0 0 24
25 557 0 1 0 0 0 0 0 0 0 0 0 0 25
26 561 0 0 1 0 0 0 0 0 0 0 0 0 26
27 549 0 0 0 1 0 0 0 0 0 0 0 0 27
28 532 0 0 0 0 1 0 0 0 0 0 0 0 28
29 526 0 0 0 0 0 1 0 0 0 0 0 0 29
30 511 0 0 0 0 0 0 1 0 0 0 0 0 30
31 499 0 0 0 0 0 0 0 1 0 0 0 0 31
32 555 0 0 0 0 0 0 0 0 1 0 0 0 32
33 565 0 0 0 0 0 0 0 0 0 1 0 0 33
34 542 0 0 0 0 0 0 0 0 0 0 1 0 34
35 527 0 0 0 0 0 0 0 0 0 0 0 1 35
36 510 0 0 0 0 0 0 0 0 0 0 0 0 36
37 514 0 1 0 0 0 0 0 0 0 0 0 0 37
38 517 0 0 1 0 0 0 0 0 0 0 0 0 38
39 508 0 0 0 1 0 0 0 0 0 0 0 0 39
40 493 0 0 0 0 1 0 0 0 0 0 0 0 40
41 490 0 0 0 0 0 1 0 0 0 0 0 0 41
42 469 0 0 0 0 0 0 1 0 0 0 0 0 42
43 478 0 0 0 0 0 0 0 1 0 0 0 0 43
44 528 1 0 0 0 0 0 0 0 1 0 0 0 44
45 534 1 0 0 0 0 0 0 0 0 1 0 0 45
46 518 1 0 0 0 0 0 0 0 0 0 1 0 46
47 506 1 0 0 0 0 0 0 0 0 0 0 1 47
48 502 1 0 0 0 0 0 0 0 0 0 0 0 48
49 516 1 1 0 0 0 0 0 0 0 0 0 0 49
50 528 1 0 1 0 0 0 0 0 0 0 0 0 50
51 533 1 0 0 1 0 0 0 0 0 0 0 0 51
52 536 1 0 0 0 1 0 0 0 0 0 0 0 52
53 537 1 0 0 0 0 1 0 0 0 0 0 0 53
54 524 1 0 0 0 0 0 1 0 0 0 0 0 54
55 536 1 0 0 0 0 0 0 1 0 0 0 0 55
56 587 1 0 0 0 0 0 0 0 1 0 0 0 56
57 597 1 0 0 0 0 0 0 0 0 1 0 0 57
58 581 1 0 0 0 0 0 0 0 0 0 1 0 58
59 564 1 0 0 0 0 0 0 0 0 0 0 1 59
60 564 1 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
613.892 39.744 -8.826 -4.337 -6.249 -12.760
M5 M6 M7 M8 M9 M10
-15.471 -24.383 -19.894 25.646 35.934 25.823
M11 t
9.711 -2.289
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-56.18 -13.98 -2.64 16.77 47.68
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 613.8918 14.4834 42.386 < 2e-16 ***
X 39.7441 12.3815 3.210 0.00242 **
M1 -8.8257 16.9282 -0.521 0.60461
M2 -4.3371 16.9008 -0.257 0.79861
M3 -6.2485 16.8795 -0.370 0.71294
M4 -12.7600 16.8642 -0.757 0.45313
M5 -15.4714 16.8551 -0.918 0.36346
M6 -24.3828 16.8520 -1.447 0.15471
M7 -19.8942 16.8551 -1.180 0.24395
M8 25.6456 16.8289 1.524 0.13438
M9 35.9342 16.8074 2.138 0.03786 *
M10 25.8228 16.7921 1.538 0.13095
M11 9.7114 16.7829 0.579 0.56565
t -2.2886 0.3208 -7.134 5.74e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 26.53 on 46 degrees of freedom
Multiple R-squared: 0.6721, Adjusted R-squared: 0.5794
F-statistic: 7.252 on 13 and 46 DF, p-value: 2.022e-07
> 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,] 3.749832e-04 7.499664e-04 9.996250e-01
[2,] 8.339102e-05 1.667820e-04 9.999166e-01
[3,] 7.829314e-06 1.565863e-05 9.999922e-01
[4,] 8.096800e-07 1.619360e-06 9.999992e-01
[5,] 1.456647e-07 2.913295e-07 9.999999e-01
[6,] 2.099077e-07 4.198154e-07 9.999998e-01
[7,] 2.225565e-05 4.451131e-05 9.999777e-01
[8,] 1.577426e-04 3.154853e-04 9.998423e-01
[9,] 1.304684e-03 2.609368e-03 9.986953e-01
[10,] 1.706078e-03 3.412155e-03 9.982939e-01
[11,] 3.584413e-03 7.168826e-03 9.964156e-01
[12,] 9.759324e-03 1.951865e-02 9.902407e-01
[13,] 1.666804e-02 3.333608e-02 9.833320e-01
[14,] 5.940350e-02 1.188070e-01 9.405965e-01
[15,] 1.627554e-01 3.255108e-01 8.372446e-01
[16,] 2.284625e-01 4.569251e-01 7.715375e-01
[17,] 3.261752e-01 6.523504e-01 6.738248e-01
[18,] 5.430668e-01 9.138665e-01 4.569332e-01
[19,] 7.863808e-01 4.272383e-01 2.136192e-01
[20,] 9.246890e-01 1.506220e-01 7.531099e-02
[21,] 9.745404e-01 5.091915e-02 2.545957e-02
[22,] 9.968982e-01 6.203560e-03 3.101780e-03
[23,] 9.999110e-01 1.780711e-04 8.903553e-05
[24,] 9.999303e-01 1.394603e-04 6.973016e-05
[25,] 9.999743e-01 5.135606e-05 2.567803e-05
[26,] 9.998124e-01 3.752811e-04 1.876405e-04
[27,] 9.979278e-01 4.144317e-03 2.072158e-03
> postscript(file="/var/www/html/rcomp/tmp/1k6e81259785433.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/2ghzs1259785433.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/3swfg1259785433.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/40eo81259785433.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/59x4h1259785433.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
-7.7774955 -13.9774955 -11.7774955 -7.9774955 -13.9774955 -8.7774955
7 8 9 10 11 12
-8.9774955 -0.2286751 -0.2286751 11.1713249 13.5713249 8.5713249
13 14 15 16 17 18
21.6856624 15.4856624 16.6856624 15.4856624 14.4856624 24.6856624
19 20 21 22 23 24
22.4856624 26.2344828 24.2344828 30.6344828 17.0344828 7.0344828
25 26 27 28 29 30
9.1488203 10.9488203 3.1488203 -5.0511797 -6.0511797 -9.8511797
31 32 33 34 35 36
-24.0511797 -11.3023593 -9.3023593 -19.9023593 -16.5023593 -21.5023593
37 38 39 40 41 42
-6.3880218 -5.5880218 -10.3880218 -16.5880218 -14.5880218 -24.3880218
43 44 45 46 47 48
-17.5880218 -50.5833031 -52.5833031 -56.1833031 -49.7833031 -41.7833031
49 50 51 52 53 54
-16.6689655 -6.8689655 2.3310345 14.1310345 20.1310345 18.3310345
55 56 57 58 59 60
28.1310345 35.8798548 37.8798548 34.2798548 35.6798548 47.6798548
> postscript(file="/var/www/html/rcomp/tmp/6rrt11259785433.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 -7.7774955 NA
1 -13.9774955 -7.7774955
2 -11.7774955 -13.9774955
3 -7.9774955 -11.7774955
4 -13.9774955 -7.9774955
5 -8.7774955 -13.9774955
6 -8.9774955 -8.7774955
7 -0.2286751 -8.9774955
8 -0.2286751 -0.2286751
9 11.1713249 -0.2286751
10 13.5713249 11.1713249
11 8.5713249 13.5713249
12 21.6856624 8.5713249
13 15.4856624 21.6856624
14 16.6856624 15.4856624
15 15.4856624 16.6856624
16 14.4856624 15.4856624
17 24.6856624 14.4856624
18 22.4856624 24.6856624
19 26.2344828 22.4856624
20 24.2344828 26.2344828
21 30.6344828 24.2344828
22 17.0344828 30.6344828
23 7.0344828 17.0344828
24 9.1488203 7.0344828
25 10.9488203 9.1488203
26 3.1488203 10.9488203
27 -5.0511797 3.1488203
28 -6.0511797 -5.0511797
29 -9.8511797 -6.0511797
30 -24.0511797 -9.8511797
31 -11.3023593 -24.0511797
32 -9.3023593 -11.3023593
33 -19.9023593 -9.3023593
34 -16.5023593 -19.9023593
35 -21.5023593 -16.5023593
36 -6.3880218 -21.5023593
37 -5.5880218 -6.3880218
38 -10.3880218 -5.5880218
39 -16.5880218 -10.3880218
40 -14.5880218 -16.5880218
41 -24.3880218 -14.5880218
42 -17.5880218 -24.3880218
43 -50.5833031 -17.5880218
44 -52.5833031 -50.5833031
45 -56.1833031 -52.5833031
46 -49.7833031 -56.1833031
47 -41.7833031 -49.7833031
48 -16.6689655 -41.7833031
49 -6.8689655 -16.6689655
50 2.3310345 -6.8689655
51 14.1310345 2.3310345
52 20.1310345 14.1310345
53 18.3310345 20.1310345
54 28.1310345 18.3310345
55 35.8798548 28.1310345
56 37.8798548 35.8798548
57 34.2798548 37.8798548
58 35.6798548 34.2798548
59 47.6798548 35.6798548
60 NA 47.6798548
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -13.9774955 -7.7774955
[2,] -11.7774955 -13.9774955
[3,] -7.9774955 -11.7774955
[4,] -13.9774955 -7.9774955
[5,] -8.7774955 -13.9774955
[6,] -8.9774955 -8.7774955
[7,] -0.2286751 -8.9774955
[8,] -0.2286751 -0.2286751
[9,] 11.1713249 -0.2286751
[10,] 13.5713249 11.1713249
[11,] 8.5713249 13.5713249
[12,] 21.6856624 8.5713249
[13,] 15.4856624 21.6856624
[14,] 16.6856624 15.4856624
[15,] 15.4856624 16.6856624
[16,] 14.4856624 15.4856624
[17,] 24.6856624 14.4856624
[18,] 22.4856624 24.6856624
[19,] 26.2344828 22.4856624
[20,] 24.2344828 26.2344828
[21,] 30.6344828 24.2344828
[22,] 17.0344828 30.6344828
[23,] 7.0344828 17.0344828
[24,] 9.1488203 7.0344828
[25,] 10.9488203 9.1488203
[26,] 3.1488203 10.9488203
[27,] -5.0511797 3.1488203
[28,] -6.0511797 -5.0511797
[29,] -9.8511797 -6.0511797
[30,] -24.0511797 -9.8511797
[31,] -11.3023593 -24.0511797
[32,] -9.3023593 -11.3023593
[33,] -19.9023593 -9.3023593
[34,] -16.5023593 -19.9023593
[35,] -21.5023593 -16.5023593
[36,] -6.3880218 -21.5023593
[37,] -5.5880218 -6.3880218
[38,] -10.3880218 -5.5880218
[39,] -16.5880218 -10.3880218
[40,] -14.5880218 -16.5880218
[41,] -24.3880218 -14.5880218
[42,] -17.5880218 -24.3880218
[43,] -50.5833031 -17.5880218
[44,] -52.5833031 -50.5833031
[45,] -56.1833031 -52.5833031
[46,] -49.7833031 -56.1833031
[47,] -41.7833031 -49.7833031
[48,] -16.6689655 -41.7833031
[49,] -6.8689655 -16.6689655
[50,] 2.3310345 -6.8689655
[51,] 14.1310345 2.3310345
[52,] 20.1310345 14.1310345
[53,] 18.3310345 20.1310345
[54,] 28.1310345 18.3310345
[55,] 35.8798548 28.1310345
[56,] 37.8798548 35.8798548
[57,] 34.2798548 37.8798548
[58,] 35.6798548 34.2798548
[59,] 47.6798548 35.6798548
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -13.9774955 -7.7774955
2 -11.7774955 -13.9774955
3 -7.9774955 -11.7774955
4 -13.9774955 -7.9774955
5 -8.7774955 -13.9774955
6 -8.9774955 -8.7774955
7 -0.2286751 -8.9774955
8 -0.2286751 -0.2286751
9 11.1713249 -0.2286751
10 13.5713249 11.1713249
11 8.5713249 13.5713249
12 21.6856624 8.5713249
13 15.4856624 21.6856624
14 16.6856624 15.4856624
15 15.4856624 16.6856624
16 14.4856624 15.4856624
17 24.6856624 14.4856624
18 22.4856624 24.6856624
19 26.2344828 22.4856624
20 24.2344828 26.2344828
21 30.6344828 24.2344828
22 17.0344828 30.6344828
23 7.0344828 17.0344828
24 9.1488203 7.0344828
25 10.9488203 9.1488203
26 3.1488203 10.9488203
27 -5.0511797 3.1488203
28 -6.0511797 -5.0511797
29 -9.8511797 -6.0511797
30 -24.0511797 -9.8511797
31 -11.3023593 -24.0511797
32 -9.3023593 -11.3023593
33 -19.9023593 -9.3023593
34 -16.5023593 -19.9023593
35 -21.5023593 -16.5023593
36 -6.3880218 -21.5023593
37 -5.5880218 -6.3880218
38 -10.3880218 -5.5880218
39 -16.5880218 -10.3880218
40 -14.5880218 -16.5880218
41 -24.3880218 -14.5880218
42 -17.5880218 -24.3880218
43 -50.5833031 -17.5880218
44 -52.5833031 -50.5833031
45 -56.1833031 -52.5833031
46 -49.7833031 -56.1833031
47 -41.7833031 -49.7833031
48 -16.6689655 -41.7833031
49 -6.8689655 -16.6689655
50 2.3310345 -6.8689655
51 14.1310345 2.3310345
52 20.1310345 14.1310345
53 18.3310345 20.1310345
54 28.1310345 18.3310345
55 35.8798548 28.1310345
56 37.8798548 35.8798548
57 34.2798548 37.8798548
58 35.6798548 34.2798548
59 47.6798548 35.6798548
> 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/73jm31259785433.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/8lppy1259785433.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/9g64l1259785433.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/10a2st1259785433.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/11gwp51259785433.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/12whtm1259785433.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/13q2w01259785433.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/145foh1259785433.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/15c1251259785433.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/168u3l1259785433.tab")
+ }
> system("convert tmp/1k6e81259785433.ps tmp/1k6e81259785433.png")
> system("convert tmp/2ghzs1259785433.ps tmp/2ghzs1259785433.png")
> system("convert tmp/3swfg1259785433.ps tmp/3swfg1259785433.png")
> system("convert tmp/40eo81259785433.ps tmp/40eo81259785433.png")
> system("convert tmp/59x4h1259785433.ps tmp/59x4h1259785433.png")
> system("convert tmp/6rrt11259785433.ps tmp/6rrt11259785433.png")
> system("convert tmp/73jm31259785433.ps tmp/73jm31259785433.png")
> system("convert tmp/8lppy1259785433.ps tmp/8lppy1259785433.png")
> system("convert tmp/9g64l1259785433.ps tmp/9g64l1259785433.png")
> system("convert tmp/10a2st1259785433.ps tmp/10a2st1259785433.png")
>
>
> proc.time()
user system elapsed
2.379 1.547 3.248