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(101.3,0,106.3,0,94,0,102.8,0,102,0,105.1,1,92.4,0,81.4,0,105.8,0,120.3,1,100.7,0,88.8,0,94.3,0,99.9,0,103.4,0,103.3,0,98.8,0,104.2,0,91.2,0,74.7,0,108.5,0,114.5,0,96.9,0,89.6,0,97.1,0,100.3,0,122.6,0,115.4,1,109,0,129.1,1,102.8,1,96.2,0,127.7,1,128.9,1,126.5,1,119.8,1,113.2,1,114.1,1,134.1,1,130,1,121.8,1,132.1,1,105.3,1,103,1,117.1,1,126.3,1,138.1,1,119.5,1,138,1,135.5,1,178.6,1,162.2,1,176.9,1,204.9,1,132.2,1,142.5,1,164.3,1,174.9,1,175.4,1,143,1),dim=c(2,60),dimnames=list(c('Omzet','Uitvoer'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Omzet','Uitvoer'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = ''
> par9 = ''
> par8 = ''
> par7 = ''
> par6 = ''
> par5 = ''
> par4 = ''
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> ylab = ''
> xlab = ''
> main = ''
> #'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
Omzet Uitvoer M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 101.3 0 1 0 0 0 0 0 0 0 0 0 0
2 106.3 0 0 1 0 0 0 0 0 0 0 0 0
3 94.0 0 0 0 1 0 0 0 0 0 0 0 0
4 102.8 0 0 0 0 1 0 0 0 0 0 0 0
5 102.0 0 0 0 0 0 1 0 0 0 0 0 0
6 105.1 1 0 0 0 0 0 1 0 0 0 0 0
7 92.4 0 0 0 0 0 0 0 1 0 0 0 0
8 81.4 0 0 0 0 0 0 0 0 1 0 0 0
9 105.8 0 0 0 0 0 0 0 0 0 1 0 0
10 120.3 1 0 0 0 0 0 0 0 0 0 1 0
11 100.7 0 0 0 0 0 0 0 0 0 0 0 1
12 88.8 0 0 0 0 0 0 0 0 0 0 0 0
13 94.3 0 1 0 0 0 0 0 0 0 0 0 0
14 99.9 0 0 1 0 0 0 0 0 0 0 0 0
15 103.4 0 0 0 1 0 0 0 0 0 0 0 0
16 103.3 0 0 0 0 1 0 0 0 0 0 0 0
17 98.8 0 0 0 0 0 1 0 0 0 0 0 0
18 104.2 0 0 0 0 0 0 1 0 0 0 0 0
19 91.2 0 0 0 0 0 0 0 1 0 0 0 0
20 74.7 0 0 0 0 0 0 0 0 1 0 0 0
21 108.5 0 0 0 0 0 0 0 0 0 1 0 0
22 114.5 0 0 0 0 0 0 0 0 0 0 1 0
23 96.9 0 0 0 0 0 0 0 0 0 0 0 1
24 89.6 0 0 0 0 0 0 0 0 0 0 0 0
25 97.1 0 1 0 0 0 0 0 0 0 0 0 0
26 100.3 0 0 1 0 0 0 0 0 0 0 0 0
27 122.6 0 0 0 1 0 0 0 0 0 0 0 0
28 115.4 1 0 0 0 1 0 0 0 0 0 0 0
29 109.0 0 0 0 0 0 1 0 0 0 0 0 0
30 129.1 1 0 0 0 0 0 1 0 0 0 0 0
31 102.8 1 0 0 0 0 0 0 1 0 0 0 0
32 96.2 0 0 0 0 0 0 0 0 1 0 0 0
33 127.7 1 0 0 0 0 0 0 0 0 1 0 0
34 128.9 1 0 0 0 0 0 0 0 0 0 1 0
35 126.5 1 0 0 0 0 0 0 0 0 0 0 1
36 119.8 1 0 0 0 0 0 0 0 0 0 0 0
37 113.2 1 1 0 0 0 0 0 0 0 0 0 0
38 114.1 1 0 1 0 0 0 0 0 0 0 0 0
39 134.1 1 0 0 1 0 0 0 0 0 0 0 0
40 130.0 1 0 0 0 1 0 0 0 0 0 0 0
41 121.8 1 0 0 0 0 1 0 0 0 0 0 0
42 132.1 1 0 0 0 0 0 1 0 0 0 0 0
43 105.3 1 0 0 0 0 0 0 1 0 0 0 0
44 103.0 1 0 0 0 0 0 0 0 1 0 0 0
45 117.1 1 0 0 0 0 0 0 0 0 1 0 0
46 126.3 1 0 0 0 0 0 0 0 0 0 1 0
47 138.1 1 0 0 0 0 0 0 0 0 0 0 1
48 119.5 1 0 0 0 0 0 0 0 0 0 0 0
49 138.0 1 1 0 0 0 0 0 0 0 0 0 0
50 135.5 1 0 1 0 0 0 0 0 0 0 0 0
51 178.6 1 0 0 1 0 0 0 0 0 0 0 0
52 162.2 1 0 0 0 1 0 0 0 0 0 0 0
53 176.9 1 0 0 0 0 1 0 0 0 0 0 0
54 204.9 1 0 0 0 0 0 1 0 0 0 0 0
55 132.2 1 0 0 0 0 0 0 1 0 0 0 0
56 142.5 1 0 0 0 0 0 0 0 1 0 0 0
57 164.3 1 0 0 0 0 0 0 0 0 1 0 0
58 174.9 1 0 0 0 0 0 0 0 0 0 1 0
59 175.4 1 0 0 0 0 0 0 0 0 0 0 1
60 143.0 1 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) Uitvoer M1 M2 M3 M4
91.176 34.940 3.628 6.068 21.388 10.600
M5 M6 M7 M8 M9 M10
16.548 15.952 -7.360 -5.592 12.540 13.852
M11
15.380
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-36.968 -11.543 -2.652 8.338 62.832
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 91.176 9.291 9.813 5.87e-13 ***
Uitvoer 34.940 5.294 6.600 3.31e-08 ***
M1 3.628 12.393 0.293 0.771
M2 6.068 12.393 0.490 0.627
M3 21.388 12.393 1.726 0.091 .
M4 10.600 12.348 0.858 0.395
M5 16.548 12.393 1.335 0.188
M6 15.952 12.393 1.287 0.204
M7 -7.360 12.348 -0.596 0.554
M8 -5.592 12.393 -0.451 0.654
M9 12.540 12.348 1.016 0.315
M10 13.852 12.393 1.118 0.269
M11 15.380 12.348 1.246 0.219
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.52 on 47 degrees of freedom
Multiple R-squared: 0.57, Adjusted R-squared: 0.4602
F-statistic: 5.192 on 12 and 47 DF, p-value: 1.860e-05
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 2.107221e-02 4.214443e-02 0.9789278
[2,] 4.733095e-03 9.466190e-03 0.9952669
[3,] 8.955368e-04 1.791074e-03 0.9991045
[4,] 1.555568e-04 3.111137e-04 0.9998444
[5,] 5.154968e-05 1.030994e-04 0.9999485
[6,] 9.251391e-06 1.850278e-05 0.9999907
[7,] 1.718440e-06 3.436879e-06 0.9999983
[8,] 3.439650e-07 6.879301e-07 0.9999997
[9,] 4.870446e-08 9.740893e-08 1.0000000
[10,] 6.426356e-09 1.285271e-08 1.0000000
[11,] 1.018870e-09 2.037740e-09 1.0000000
[12,] 7.784842e-07 1.556968e-06 0.9999992
[13,] 3.479978e-07 6.959956e-07 0.9999997
[14,] 1.428702e-07 2.857404e-07 0.9999999
[15,] 9.323925e-07 1.864785e-06 0.9999991
[16,] 2.567499e-07 5.134999e-07 0.9999997
[17,] 3.905163e-07 7.810326e-07 0.9999996
[18,] 1.534115e-07 3.068229e-07 0.9999998
[19,] 5.304257e-08 1.060851e-07 0.9999999
[20,] 4.951446e-08 9.902892e-08 1.0000000
[21,] 3.086753e-08 6.173506e-08 1.0000000
[22,] 9.552180e-09 1.910436e-08 1.0000000
[23,] 2.997622e-09 5.995245e-09 1.0000000
[24,] 3.212341e-09 6.424683e-09 1.0000000
[25,] 2.225932e-09 4.451864e-09 1.0000000
[26,] 3.790440e-09 7.580879e-09 1.0000000
[27,] 2.691203e-07 5.382406e-07 0.9999997
[28,] 1.424542e-07 2.849085e-07 0.9999999
[29,] 2.405609e-07 4.811217e-07 0.9999998
> postscript(file="/var/www/html/rcomp/tmp/1wg471258567496.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/23d0j1258567496.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/35fu61258567496.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/4mrgf1258567496.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/5s33f1258567496.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
6.4958824 9.0558824 -18.5641176 1.0238235 -5.7241176 -36.9679412
7 8 9 10 11 12
8.5838235 -4.1841176 2.0838235 -19.6679412 -5.8561765 -2.3761765
13 14 15 16 17 18
-0.5041176 2.6558824 -9.1641176 1.5238235 -8.9241176 -2.9282353
19 20 21 22 23 24
7.3838235 -10.8841176 4.7838235 9.4717647 -9.6561765 -1.5761765
25 26 27 28 29 30
2.2958824 3.0558824 10.0358824 -21.3158824 1.2758824 -12.9679412
31 32 33 34 35 36
-15.9558824 10.6158824 -10.9558824 -11.0679412 -14.9958824 -6.3158824
37 38 39 40 41 42
-16.5438235 -18.0838235 -13.4038235 -6.7158824 -20.8638235 -9.9679412
43 44 45 46 47 48
-13.4558824 -17.5238235 -21.5558824 -13.6679412 -3.3958824 -6.6158824
49 50 51 52 53 54
8.2561765 3.3161765 31.0961765 25.4841176 34.2361765 62.8320588
55 56 57 58 59 60
13.4441176 21.9761765 25.6441176 34.9320588 33.9041176 16.8841176
> postscript(file="/var/www/html/rcomp/tmp/6bx4k1258567496.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 6.4958824 NA
1 9.0558824 6.4958824
2 -18.5641176 9.0558824
3 1.0238235 -18.5641176
4 -5.7241176 1.0238235
5 -36.9679412 -5.7241176
6 8.5838235 -36.9679412
7 -4.1841176 8.5838235
8 2.0838235 -4.1841176
9 -19.6679412 2.0838235
10 -5.8561765 -19.6679412
11 -2.3761765 -5.8561765
12 -0.5041176 -2.3761765
13 2.6558824 -0.5041176
14 -9.1641176 2.6558824
15 1.5238235 -9.1641176
16 -8.9241176 1.5238235
17 -2.9282353 -8.9241176
18 7.3838235 -2.9282353
19 -10.8841176 7.3838235
20 4.7838235 -10.8841176
21 9.4717647 4.7838235
22 -9.6561765 9.4717647
23 -1.5761765 -9.6561765
24 2.2958824 -1.5761765
25 3.0558824 2.2958824
26 10.0358824 3.0558824
27 -21.3158824 10.0358824
28 1.2758824 -21.3158824
29 -12.9679412 1.2758824
30 -15.9558824 -12.9679412
31 10.6158824 -15.9558824
32 -10.9558824 10.6158824
33 -11.0679412 -10.9558824
34 -14.9958824 -11.0679412
35 -6.3158824 -14.9958824
36 -16.5438235 -6.3158824
37 -18.0838235 -16.5438235
38 -13.4038235 -18.0838235
39 -6.7158824 -13.4038235
40 -20.8638235 -6.7158824
41 -9.9679412 -20.8638235
42 -13.4558824 -9.9679412
43 -17.5238235 -13.4558824
44 -21.5558824 -17.5238235
45 -13.6679412 -21.5558824
46 -3.3958824 -13.6679412
47 -6.6158824 -3.3958824
48 8.2561765 -6.6158824
49 3.3161765 8.2561765
50 31.0961765 3.3161765
51 25.4841176 31.0961765
52 34.2361765 25.4841176
53 62.8320588 34.2361765
54 13.4441176 62.8320588
55 21.9761765 13.4441176
56 25.6441176 21.9761765
57 34.9320588 25.6441176
58 33.9041176 34.9320588
59 16.8841176 33.9041176
60 NA 16.8841176
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 9.0558824 6.4958824
[2,] -18.5641176 9.0558824
[3,] 1.0238235 -18.5641176
[4,] -5.7241176 1.0238235
[5,] -36.9679412 -5.7241176
[6,] 8.5838235 -36.9679412
[7,] -4.1841176 8.5838235
[8,] 2.0838235 -4.1841176
[9,] -19.6679412 2.0838235
[10,] -5.8561765 -19.6679412
[11,] -2.3761765 -5.8561765
[12,] -0.5041176 -2.3761765
[13,] 2.6558824 -0.5041176
[14,] -9.1641176 2.6558824
[15,] 1.5238235 -9.1641176
[16,] -8.9241176 1.5238235
[17,] -2.9282353 -8.9241176
[18,] 7.3838235 -2.9282353
[19,] -10.8841176 7.3838235
[20,] 4.7838235 -10.8841176
[21,] 9.4717647 4.7838235
[22,] -9.6561765 9.4717647
[23,] -1.5761765 -9.6561765
[24,] 2.2958824 -1.5761765
[25,] 3.0558824 2.2958824
[26,] 10.0358824 3.0558824
[27,] -21.3158824 10.0358824
[28,] 1.2758824 -21.3158824
[29,] -12.9679412 1.2758824
[30,] -15.9558824 -12.9679412
[31,] 10.6158824 -15.9558824
[32,] -10.9558824 10.6158824
[33,] -11.0679412 -10.9558824
[34,] -14.9958824 -11.0679412
[35,] -6.3158824 -14.9958824
[36,] -16.5438235 -6.3158824
[37,] -18.0838235 -16.5438235
[38,] -13.4038235 -18.0838235
[39,] -6.7158824 -13.4038235
[40,] -20.8638235 -6.7158824
[41,] -9.9679412 -20.8638235
[42,] -13.4558824 -9.9679412
[43,] -17.5238235 -13.4558824
[44,] -21.5558824 -17.5238235
[45,] -13.6679412 -21.5558824
[46,] -3.3958824 -13.6679412
[47,] -6.6158824 -3.3958824
[48,] 8.2561765 -6.6158824
[49,] 3.3161765 8.2561765
[50,] 31.0961765 3.3161765
[51,] 25.4841176 31.0961765
[52,] 34.2361765 25.4841176
[53,] 62.8320588 34.2361765
[54,] 13.4441176 62.8320588
[55,] 21.9761765 13.4441176
[56,] 25.6441176 21.9761765
[57,] 34.9320588 25.6441176
[58,] 33.9041176 34.9320588
[59,] 16.8841176 33.9041176
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 9.0558824 6.4958824
2 -18.5641176 9.0558824
3 1.0238235 -18.5641176
4 -5.7241176 1.0238235
5 -36.9679412 -5.7241176
6 8.5838235 -36.9679412
7 -4.1841176 8.5838235
8 2.0838235 -4.1841176
9 -19.6679412 2.0838235
10 -5.8561765 -19.6679412
11 -2.3761765 -5.8561765
12 -0.5041176 -2.3761765
13 2.6558824 -0.5041176
14 -9.1641176 2.6558824
15 1.5238235 -9.1641176
16 -8.9241176 1.5238235
17 -2.9282353 -8.9241176
18 7.3838235 -2.9282353
19 -10.8841176 7.3838235
20 4.7838235 -10.8841176
21 9.4717647 4.7838235
22 -9.6561765 9.4717647
23 -1.5761765 -9.6561765
24 2.2958824 -1.5761765
25 3.0558824 2.2958824
26 10.0358824 3.0558824
27 -21.3158824 10.0358824
28 1.2758824 -21.3158824
29 -12.9679412 1.2758824
30 -15.9558824 -12.9679412
31 10.6158824 -15.9558824
32 -10.9558824 10.6158824
33 -11.0679412 -10.9558824
34 -14.9958824 -11.0679412
35 -6.3158824 -14.9958824
36 -16.5438235 -6.3158824
37 -18.0838235 -16.5438235
38 -13.4038235 -18.0838235
39 -6.7158824 -13.4038235
40 -20.8638235 -6.7158824
41 -9.9679412 -20.8638235
42 -13.4558824 -9.9679412
43 -17.5238235 -13.4558824
44 -21.5558824 -17.5238235
45 -13.6679412 -21.5558824
46 -3.3958824 -13.6679412
47 -6.6158824 -3.3958824
48 8.2561765 -6.6158824
49 3.3161765 8.2561765
50 31.0961765 3.3161765
51 25.4841176 31.0961765
52 34.2361765 25.4841176
53 62.8320588 34.2361765
54 13.4441176 62.8320588
55 21.9761765 13.4441176
56 25.6441176 21.9761765
57 34.9320588 25.6441176
58 33.9041176 34.9320588
59 16.8841176 33.9041176
> 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/7rlfr1258567496.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/856291258567496.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/9itzs1258567496.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/10z2z71258567496.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/11v3es1258567496.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/12dp2c1258567496.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/13c2f61258567496.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/1435i31258567496.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/15adr11258567497.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/16ueoa1258567497.tab")
+ }
>
> system("convert tmp/1wg471258567496.ps tmp/1wg471258567496.png")
> system("convert tmp/23d0j1258567496.ps tmp/23d0j1258567496.png")
> system("convert tmp/35fu61258567496.ps tmp/35fu61258567496.png")
> system("convert tmp/4mrgf1258567496.ps tmp/4mrgf1258567496.png")
> system("convert tmp/5s33f1258567496.ps tmp/5s33f1258567496.png")
> system("convert tmp/6bx4k1258567496.ps tmp/6bx4k1258567496.png")
> system("convert tmp/7rlfr1258567496.ps tmp/7rlfr1258567496.png")
> system("convert tmp/856291258567496.ps tmp/856291258567496.png")
> system("convert tmp/9itzs1258567496.ps tmp/9itzs1258567496.png")
> system("convert tmp/10z2z71258567496.ps tmp/10z2z71258567496.png")
>
>
> proc.time()
user system elapsed
2.413 1.556 3.591