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(9.3,145.0,8.7,137.7,8.2,148.3,8.3,152.2,8.5,169.4,8.6,168.6,8.5,161.1,8.2,174.1,8.1,179.0,7.9,190.6,8.6,190.0,8.7,181.6,8.7,174.8,8.5,180.5,8.4,196.8,8.5,193.8,8.7,197.0,8.7,216.3,8.6,221.4,8.5,217.9,8.3,229.7,8,227.4,8.2,204.2,8.1,196.6,8.1,198.8,8,207.5,7.9,190.7,7.9,201.6,8,210.5,8,223.5,7.9,223.8,8,231.2,7.7,244.0,7.2,234.7,7.5,250.2,7.3,265.7,7,287.6,7,283.3,7,295.4,7.2,312.3,7.3,333.8,7.1,347.7,6.8,383.2,6.4,407.1,6.1,413.6,6.5,362.7,7.7,321.9,7.9,239.4,7.5,191.0,6.9,159.7,6.6,163.4,6.9,157.6,7.7,166.2,8,176.7,8,198.3,7.7,226.2,7.3,216.2,7.4,235.9,8.1,226.9,8.3,242.3,8.2,253.1),dim=c(2,61),dimnames=list(c('TW','GP'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('TW','GP'),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 = '2'
> #'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
GP TW M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 145.0 9.3 1 0 0 0 0 0 0 0 0 0 0
2 137.7 8.7 0 1 0 0 0 0 0 0 0 0 0
3 148.3 8.2 0 0 1 0 0 0 0 0 0 0 0
4 152.2 8.3 0 0 0 1 0 0 0 0 0 0 0
5 169.4 8.5 0 0 0 0 1 0 0 0 0 0 0
6 168.6 8.6 0 0 0 0 0 1 0 0 0 0 0
7 161.1 8.5 0 0 0 0 0 0 1 0 0 0 0
8 174.1 8.2 0 0 0 0 0 0 0 1 0 0 0
9 179.0 8.1 0 0 0 0 0 0 0 0 1 0 0
10 190.6 7.9 0 0 0 0 0 0 0 0 0 1 0
11 190.0 8.6 0 0 0 0 0 0 0 0 0 0 1
12 181.6 8.7 0 0 0 0 0 0 0 0 0 0 0
13 174.8 8.7 1 0 0 0 0 0 0 0 0 0 0
14 180.5 8.5 0 1 0 0 0 0 0 0 0 0 0
15 196.8 8.4 0 0 1 0 0 0 0 0 0 0 0
16 193.8 8.5 0 0 0 1 0 0 0 0 0 0 0
17 197.0 8.7 0 0 0 0 1 0 0 0 0 0 0
18 216.3 8.7 0 0 0 0 0 1 0 0 0 0 0
19 221.4 8.6 0 0 0 0 0 0 1 0 0 0 0
20 217.9 8.5 0 0 0 0 0 0 0 1 0 0 0
21 229.7 8.3 0 0 0 0 0 0 0 0 1 0 0
22 227.4 8.0 0 0 0 0 0 0 0 0 0 1 0
23 204.2 8.2 0 0 0 0 0 0 0 0 0 0 1
24 196.6 8.1 0 0 0 0 0 0 0 0 0 0 0
25 198.8 8.1 1 0 0 0 0 0 0 0 0 0 0
26 207.5 8.0 0 1 0 0 0 0 0 0 0 0 0
27 190.7 7.9 0 0 1 0 0 0 0 0 0 0 0
28 201.6 7.9 0 0 0 1 0 0 0 0 0 0 0
29 210.5 8.0 0 0 0 0 1 0 0 0 0 0 0
30 223.5 8.0 0 0 0 0 0 1 0 0 0 0 0
31 223.8 7.9 0 0 0 0 0 0 1 0 0 0 0
32 231.2 8.0 0 0 0 0 0 0 0 1 0 0 0
33 244.0 7.7 0 0 0 0 0 0 0 0 1 0 0
34 234.7 7.2 0 0 0 0 0 0 0 0 0 1 0
35 250.2 7.5 0 0 0 0 0 0 0 0 0 0 1
36 265.7 7.3 0 0 0 0 0 0 0 0 0 0 0
37 287.6 7.0 1 0 0 0 0 0 0 0 0 0 0
38 283.3 7.0 0 1 0 0 0 0 0 0 0 0 0
39 295.4 7.0 0 0 1 0 0 0 0 0 0 0 0
40 312.3 7.2 0 0 0 1 0 0 0 0 0 0 0
41 333.8 7.3 0 0 0 0 1 0 0 0 0 0 0
42 347.7 7.1 0 0 0 0 0 1 0 0 0 0 0
43 383.2 6.8 0 0 0 0 0 0 1 0 0 0 0
44 407.1 6.4 0 0 0 0 0 0 0 1 0 0 0
45 413.6 6.1 0 0 0 0 0 0 0 0 1 0 0
46 362.7 6.5 0 0 0 0 0 0 0 0 0 1 0
47 321.9 7.7 0 0 0 0 0 0 0 0 0 0 1
48 239.4 7.9 0 0 0 0 0 0 0 0 0 0 0
49 191.0 7.5 1 0 0 0 0 0 0 0 0 0 0
50 159.7 6.9 0 1 0 0 0 0 0 0 0 0 0
51 163.4 6.6 0 0 1 0 0 0 0 0 0 0 0
52 157.6 6.9 0 0 0 1 0 0 0 0 0 0 0
53 166.2 7.7 0 0 0 0 1 0 0 0 0 0 0
54 176.7 8.0 0 0 0 0 0 1 0 0 0 0 0
55 198.3 8.0 0 0 0 0 0 0 1 0 0 0 0
56 226.2 7.7 0 0 0 0 0 0 0 1 0 0 0
57 216.2 7.3 0 0 0 0 0 0 0 0 1 0 0
58 235.9 7.4 0 0 0 0 0 0 0 0 0 1 0
59 226.9 8.1 0 0 0 0 0 0 0 0 0 0 1
60 242.3 8.3 0 0 0 0 0 0 0 0 0 0 0
61 253.1 8.2 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) TW M1 M2 M3 M4
778.153 -68.615 -11.705 -47.847 -56.390 -42.204
M5 M6 M7 M8 M9 M10
-11.112 2.812 5.579 5.596 -7.044 -20.146
M11
10.775
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-105.5069 -24.1196 0.3933 32.2810 70.3758
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 778.153 80.142 9.710 6.59e-13 ***
TW -68.615 9.588 -7.156 4.24e-09 ***
M1 -11.705 28.738 -0.407 0.6856
M2 -47.847 30.095 -1.590 0.1184
M3 -56.390 30.302 -1.861 0.0689 .
M4 -42.204 30.145 -1.400 0.1679
M5 -11.112 30.008 -0.370 0.7128
M6 2.812 30.008 0.094 0.9257
M7 5.579 30.023 0.186 0.8534
M8 5.596 30.145 0.186 0.8535
M9 -7.044 30.484 -0.231 0.8182
M10 -20.146 30.667 -0.657 0.5144
M11 10.775 30.010 0.359 0.7211
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 47.45 on 48 degrees of freedom
Multiple R-squared: 0.5667, Adjusted R-squared: 0.4584
F-statistic: 5.232 on 12 and 48 DF, p-value: 1.560e-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.350753e-01 4.701506e-01 0.7649247
[2,] 1.375901e-01 2.751801e-01 0.8624099
[3,] 1.075400e-01 2.150801e-01 0.8924600
[4,] 1.027983e-01 2.055966e-01 0.8972017
[5,] 6.635669e-02 1.327134e-01 0.9336433
[6,] 4.757089e-02 9.514178e-02 0.9524291
[7,] 2.888827e-02 5.777654e-02 0.9711117
[8,] 1.691417e-02 3.382835e-02 0.9830858
[9,] 9.639287e-03 1.927857e-02 0.9903607
[10,] 6.628171e-03 1.325634e-02 0.9933718
[11,] 5.415486e-03 1.083097e-02 0.9945845
[12,] 2.887743e-03 5.775486e-03 0.9971123
[13,] 1.480823e-03 2.961645e-03 0.9985192
[14,] 6.251843e-04 1.250369e-03 0.9993748
[15,] 2.517377e-04 5.034753e-04 0.9997483
[16,] 9.778322e-05 1.955664e-04 0.9999022
[17,] 4.212735e-05 8.425470e-05 0.9999579
[18,] 1.806300e-05 3.612599e-05 0.9999819
[19,] 6.445914e-06 1.289183e-05 0.9999936
[20,] 3.150893e-06 6.301786e-06 0.9999968
[21,] 1.883176e-06 3.766353e-06 0.9999981
[22,] 9.189252e-07 1.837850e-06 0.9999991
[23,] 1.554607e-06 3.109213e-06 0.9999984
[24,] 2.446124e-05 4.892247e-05 0.9999755
[25,] 2.808599e-03 5.617199e-03 0.9971914
[26,] 1.687504e-02 3.375008e-02 0.9831250
[27,] 1.657393e-02 3.314787e-02 0.9834261
[28,] 1.493434e-02 2.986869e-02 0.9850657
[29,] 9.371503e-03 1.874301e-02 0.9906285
[30,] 1.450720e-02 2.901441e-02 0.9854928
> postscript(file="/var/www/html/rcomp/tmp/118xb1260981371.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/21bv11260981371.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/3fy6n1260981371.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/4snfl1260981371.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/5lrmz1260981371.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
16.6670007 4.3408234 -10.8235482 -14.2481311 -14.4172969 -22.2804226
7 8 9 10 11 12
-39.4081311 -47.0095883 -36.3312568 -25.3527140 -8.8435482 0.3933261
13 14 15 16 17 18
5.2982575 33.4179090 51.3993662 41.0747833 26.9056175 32.2810346
19 20 21 22 23 24
27.7533261 17.3747833 28.0916576 18.3087432 -22.0893770 -25.7754171
25 26 27 28 29 30
-11.8704857 26.1106230 10.9920802 7.7060401 -7.6245829 -8.5491658
31 32 33 34 35 36
-17.8768743 -3.6325027 1.2229144 -29.2829144 -24.1195774 -11.5670747
37 38 39 40 41 42
1.4534851 33.2960510 53.9389654 70.3758397 67.6452167 53.8977195
43 44 45 46 47 48
66.0470965 62.4841821 61.0395992 50.6868852 61.3033370 3.3016685
49 50 51 52 53 54
-60.8392289 -97.1654062 -105.5068634 -104.9085319 -72.5089545 -55.3491658
55 56 57 58 59 60
-36.5154171 -29.2168743 -54.0229144 -14.3600000 -6.2508342 33.6474973
61
49.2909715
> postscript(file="/var/www/html/rcomp/tmp/6juoe1260981371.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 16.6670007 NA
1 4.3408234 16.6670007
2 -10.8235482 4.3408234
3 -14.2481311 -10.8235482
4 -14.4172969 -14.2481311
5 -22.2804226 -14.4172969
6 -39.4081311 -22.2804226
7 -47.0095883 -39.4081311
8 -36.3312568 -47.0095883
9 -25.3527140 -36.3312568
10 -8.8435482 -25.3527140
11 0.3933261 -8.8435482
12 5.2982575 0.3933261
13 33.4179090 5.2982575
14 51.3993662 33.4179090
15 41.0747833 51.3993662
16 26.9056175 41.0747833
17 32.2810346 26.9056175
18 27.7533261 32.2810346
19 17.3747833 27.7533261
20 28.0916576 17.3747833
21 18.3087432 28.0916576
22 -22.0893770 18.3087432
23 -25.7754171 -22.0893770
24 -11.8704857 -25.7754171
25 26.1106230 -11.8704857
26 10.9920802 26.1106230
27 7.7060401 10.9920802
28 -7.6245829 7.7060401
29 -8.5491658 -7.6245829
30 -17.8768743 -8.5491658
31 -3.6325027 -17.8768743
32 1.2229144 -3.6325027
33 -29.2829144 1.2229144
34 -24.1195774 -29.2829144
35 -11.5670747 -24.1195774
36 1.4534851 -11.5670747
37 33.2960510 1.4534851
38 53.9389654 33.2960510
39 70.3758397 53.9389654
40 67.6452167 70.3758397
41 53.8977195 67.6452167
42 66.0470965 53.8977195
43 62.4841821 66.0470965
44 61.0395992 62.4841821
45 50.6868852 61.0395992
46 61.3033370 50.6868852
47 3.3016685 61.3033370
48 -60.8392289 3.3016685
49 -97.1654062 -60.8392289
50 -105.5068634 -97.1654062
51 -104.9085319 -105.5068634
52 -72.5089545 -104.9085319
53 -55.3491658 -72.5089545
54 -36.5154171 -55.3491658
55 -29.2168743 -36.5154171
56 -54.0229144 -29.2168743
57 -14.3600000 -54.0229144
58 -6.2508342 -14.3600000
59 33.6474973 -6.2508342
60 49.2909715 33.6474973
61 NA 49.2909715
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4.3408234 16.6670007
[2,] -10.8235482 4.3408234
[3,] -14.2481311 -10.8235482
[4,] -14.4172969 -14.2481311
[5,] -22.2804226 -14.4172969
[6,] -39.4081311 -22.2804226
[7,] -47.0095883 -39.4081311
[8,] -36.3312568 -47.0095883
[9,] -25.3527140 -36.3312568
[10,] -8.8435482 -25.3527140
[11,] 0.3933261 -8.8435482
[12,] 5.2982575 0.3933261
[13,] 33.4179090 5.2982575
[14,] 51.3993662 33.4179090
[15,] 41.0747833 51.3993662
[16,] 26.9056175 41.0747833
[17,] 32.2810346 26.9056175
[18,] 27.7533261 32.2810346
[19,] 17.3747833 27.7533261
[20,] 28.0916576 17.3747833
[21,] 18.3087432 28.0916576
[22,] -22.0893770 18.3087432
[23,] -25.7754171 -22.0893770
[24,] -11.8704857 -25.7754171
[25,] 26.1106230 -11.8704857
[26,] 10.9920802 26.1106230
[27,] 7.7060401 10.9920802
[28,] -7.6245829 7.7060401
[29,] -8.5491658 -7.6245829
[30,] -17.8768743 -8.5491658
[31,] -3.6325027 -17.8768743
[32,] 1.2229144 -3.6325027
[33,] -29.2829144 1.2229144
[34,] -24.1195774 -29.2829144
[35,] -11.5670747 -24.1195774
[36,] 1.4534851 -11.5670747
[37,] 33.2960510 1.4534851
[38,] 53.9389654 33.2960510
[39,] 70.3758397 53.9389654
[40,] 67.6452167 70.3758397
[41,] 53.8977195 67.6452167
[42,] 66.0470965 53.8977195
[43,] 62.4841821 66.0470965
[44,] 61.0395992 62.4841821
[45,] 50.6868852 61.0395992
[46,] 61.3033370 50.6868852
[47,] 3.3016685 61.3033370
[48,] -60.8392289 3.3016685
[49,] -97.1654062 -60.8392289
[50,] -105.5068634 -97.1654062
[51,] -104.9085319 -105.5068634
[52,] -72.5089545 -104.9085319
[53,] -55.3491658 -72.5089545
[54,] -36.5154171 -55.3491658
[55,] -29.2168743 -36.5154171
[56,] -54.0229144 -29.2168743
[57,] -14.3600000 -54.0229144
[58,] -6.2508342 -14.3600000
[59,] 33.6474973 -6.2508342
[60,] 49.2909715 33.6474973
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4.3408234 16.6670007
2 -10.8235482 4.3408234
3 -14.2481311 -10.8235482
4 -14.4172969 -14.2481311
5 -22.2804226 -14.4172969
6 -39.4081311 -22.2804226
7 -47.0095883 -39.4081311
8 -36.3312568 -47.0095883
9 -25.3527140 -36.3312568
10 -8.8435482 -25.3527140
11 0.3933261 -8.8435482
12 5.2982575 0.3933261
13 33.4179090 5.2982575
14 51.3993662 33.4179090
15 41.0747833 51.3993662
16 26.9056175 41.0747833
17 32.2810346 26.9056175
18 27.7533261 32.2810346
19 17.3747833 27.7533261
20 28.0916576 17.3747833
21 18.3087432 28.0916576
22 -22.0893770 18.3087432
23 -25.7754171 -22.0893770
24 -11.8704857 -25.7754171
25 26.1106230 -11.8704857
26 10.9920802 26.1106230
27 7.7060401 10.9920802
28 -7.6245829 7.7060401
29 -8.5491658 -7.6245829
30 -17.8768743 -8.5491658
31 -3.6325027 -17.8768743
32 1.2229144 -3.6325027
33 -29.2829144 1.2229144
34 -24.1195774 -29.2829144
35 -11.5670747 -24.1195774
36 1.4534851 -11.5670747
37 33.2960510 1.4534851
38 53.9389654 33.2960510
39 70.3758397 53.9389654
40 67.6452167 70.3758397
41 53.8977195 67.6452167
42 66.0470965 53.8977195
43 62.4841821 66.0470965
44 61.0395992 62.4841821
45 50.6868852 61.0395992
46 61.3033370 50.6868852
47 3.3016685 61.3033370
48 -60.8392289 3.3016685
49 -97.1654062 -60.8392289
50 -105.5068634 -97.1654062
51 -104.9085319 -105.5068634
52 -72.5089545 -104.9085319
53 -55.3491658 -72.5089545
54 -36.5154171 -55.3491658
55 -29.2168743 -36.5154171
56 -54.0229144 -29.2168743
57 -14.3600000 -54.0229144
58 -6.2508342 -14.3600000
59 33.6474973 -6.2508342
60 49.2909715 33.6474973
> 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/7snur1260981371.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/8nc471260981371.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/9rix71260981371.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/1006n01260981371.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/116t6a1260981371.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/12zx4t1260981371.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/13u2bm1260981371.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/1424vd1260981371.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/15a2nk1260981371.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/166m1x1260981371.tab")
+ }
>
> try(system("convert tmp/118xb1260981371.ps tmp/118xb1260981371.png",intern=TRUE))
character(0)
> try(system("convert tmp/21bv11260981371.ps tmp/21bv11260981371.png",intern=TRUE))
character(0)
> try(system("convert tmp/3fy6n1260981371.ps tmp/3fy6n1260981371.png",intern=TRUE))
character(0)
> try(system("convert tmp/4snfl1260981371.ps tmp/4snfl1260981371.png",intern=TRUE))
character(0)
> try(system("convert tmp/5lrmz1260981371.ps tmp/5lrmz1260981371.png",intern=TRUE))
character(0)
> try(system("convert tmp/6juoe1260981371.ps tmp/6juoe1260981371.png",intern=TRUE))
character(0)
> try(system("convert tmp/7snur1260981371.ps tmp/7snur1260981371.png",intern=TRUE))
character(0)
> try(system("convert tmp/8nc471260981371.ps tmp/8nc471260981371.png",intern=TRUE))
character(0)
> try(system("convert tmp/9rix71260981371.ps tmp/9rix71260981371.png",intern=TRUE))
character(0)
> try(system("convert tmp/1006n01260981371.ps tmp/1006n01260981371.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.353 1.498 3.991