R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(262,302,218,262,175,218,100,175,77,100,43,77,47,43,49,47,69,49,152,69,205,152,246,205,294,246,242,294,181,242,107,181,56,107,49,56,47,49,47,47,71,47,151,71,244,151,280,244,230,280,185,230,148,185,98,148,61,98,46,61,45,46,55,45,48,55,115,48,185,115,276,185,220,276,181,220,151,181,83,151,55,83,49,55,42,49,46,42,74,46,103,74,200,103,237,200,247,237,215,247,182,215,80,182,46,80,65,46,40,65,44,40,63,44,85,63,185,85,247,185,231,247,167,231,117,167,79,117,45,79,40,45,38,40,41,38,69,41,152,69,232,152,282,232,255,282,161,255,107,161,53,107,40,53,39,40,34,39,35,34,56,35,97,56,210,97,260,210,257,260,210,257,125,210,80,125,42,80,35,42,31,35,32,31,50,32,92,50,189,92,256,189,250,256,198,250,136,198,73,136,39,73,32,39,30,32,31,30,45,31),dim=c(2,105),dimnames=list(c('Gasverbruik','verbruik-1'),1:105))
> y <- array(NA,dim=c(2,105),dimnames=list(c('Gasverbruik','verbruik-1'),1:105))
> 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'
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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
Gasverbruik verbruik-1 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 262 302 1 0 0 0 0 0 0 0 0 0 0 1
2 218 262 0 1 0 0 0 0 0 0 0 0 0 2
3 175 218 0 0 1 0 0 0 0 0 0 0 0 3
4 100 175 0 0 0 1 0 0 0 0 0 0 0 4
5 77 100 0 0 0 0 1 0 0 0 0 0 0 5
6 43 77 0 0 0 0 0 1 0 0 0 0 0 6
7 47 43 0 0 0 0 0 0 1 0 0 0 0 7
8 49 47 0 0 0 0 0 0 0 1 0 0 0 8
9 69 49 0 0 0 0 0 0 0 0 1 0 0 9
10 152 69 0 0 0 0 0 0 0 0 0 1 0 10
11 205 152 0 0 0 0 0 0 0 0 0 0 1 11
12 246 205 0 0 0 0 0 0 0 0 0 0 0 12
13 294 246 1 0 0 0 0 0 0 0 0 0 0 13
14 242 294 0 1 0 0 0 0 0 0 0 0 0 14
15 181 242 0 0 1 0 0 0 0 0 0 0 0 15
16 107 181 0 0 0 1 0 0 0 0 0 0 0 16
17 56 107 0 0 0 0 1 0 0 0 0 0 0 17
18 49 56 0 0 0 0 0 1 0 0 0 0 0 18
19 47 49 0 0 0 0 0 0 1 0 0 0 0 19
20 47 47 0 0 0 0 0 0 0 1 0 0 0 20
21 71 47 0 0 0 0 0 0 0 0 1 0 0 21
22 151 71 0 0 0 0 0 0 0 0 0 1 0 22
23 244 151 0 0 0 0 0 0 0 0 0 0 1 23
24 280 244 0 0 0 0 0 0 0 0 0 0 0 24
25 230 280 1 0 0 0 0 0 0 0 0 0 0 25
26 185 230 0 1 0 0 0 0 0 0 0 0 0 26
27 148 185 0 0 1 0 0 0 0 0 0 0 0 27
28 98 148 0 0 0 1 0 0 0 0 0 0 0 28
29 61 98 0 0 0 0 1 0 0 0 0 0 0 29
30 46 61 0 0 0 0 0 1 0 0 0 0 0 30
31 45 46 0 0 0 0 0 0 1 0 0 0 0 31
32 55 45 0 0 0 0 0 0 0 1 0 0 0 32
33 48 55 0 0 0 0 0 0 0 0 1 0 0 33
34 115 48 0 0 0 0 0 0 0 0 0 1 0 34
35 185 115 0 0 0 0 0 0 0 0 0 0 1 35
36 276 185 0 0 0 0 0 0 0 0 0 0 0 36
37 220 276 1 0 0 0 0 0 0 0 0 0 0 37
38 181 220 0 1 0 0 0 0 0 0 0 0 0 38
39 151 181 0 0 1 0 0 0 0 0 0 0 0 39
40 83 151 0 0 0 1 0 0 0 0 0 0 0 40
41 55 83 0 0 0 0 1 0 0 0 0 0 0 41
42 49 55 0 0 0 0 0 1 0 0 0 0 0 42
43 42 49 0 0 0 0 0 0 1 0 0 0 0 43
44 46 42 0 0 0 0 0 0 0 1 0 0 0 44
45 74 46 0 0 0 0 0 0 0 0 1 0 0 45
46 103 74 0 0 0 0 0 0 0 0 0 1 0 46
47 200 103 0 0 0 0 0 0 0 0 0 0 1 47
48 237 200 0 0 0 0 0 0 0 0 0 0 0 48
49 247 237 1 0 0 0 0 0 0 0 0 0 0 49
50 215 247 0 1 0 0 0 0 0 0 0 0 0 50
51 182 215 0 0 1 0 0 0 0 0 0 0 0 51
52 80 182 0 0 0 1 0 0 0 0 0 0 0 52
53 46 80 0 0 0 0 1 0 0 0 0 0 0 53
54 65 46 0 0 0 0 0 1 0 0 0 0 0 54
55 40 65 0 0 0 0 0 0 1 0 0 0 0 55
56 44 40 0 0 0 0 0 0 0 1 0 0 0 56
57 63 44 0 0 0 0 0 0 0 0 1 0 0 57
58 85 63 0 0 0 0 0 0 0 0 0 1 0 58
59 185 85 0 0 0 0 0 0 0 0 0 0 1 59
60 247 185 0 0 0 0 0 0 0 0 0 0 0 60
61 231 247 1 0 0 0 0 0 0 0 0 0 0 61
62 167 231 0 1 0 0 0 0 0 0 0 0 0 62
63 117 167 0 0 1 0 0 0 0 0 0 0 0 63
64 79 117 0 0 0 1 0 0 0 0 0 0 0 64
65 45 79 0 0 0 0 1 0 0 0 0 0 0 65
66 40 45 0 0 0 0 0 1 0 0 0 0 0 66
67 38 40 0 0 0 0 0 0 1 0 0 0 0 67
68 41 38 0 0 0 0 0 0 0 1 0 0 0 68
69 69 41 0 0 0 0 0 0 0 0 1 0 0 69
70 152 69 0 0 0 0 0 0 0 0 0 1 0 70
71 232 152 0 0 0 0 0 0 0 0 0 0 1 71
72 282 232 0 0 0 0 0 0 0 0 0 0 0 72
73 255 282 1 0 0 0 0 0 0 0 0 0 0 73
74 161 255 0 1 0 0 0 0 0 0 0 0 0 74
75 107 161 0 0 1 0 0 0 0 0 0 0 0 75
76 53 107 0 0 0 1 0 0 0 0 0 0 0 76
77 40 53 0 0 0 0 1 0 0 0 0 0 0 77
78 39 40 0 0 0 0 0 1 0 0 0 0 0 78
79 34 39 0 0 0 0 0 0 1 0 0 0 0 79
80 35 34 0 0 0 0 0 0 0 1 0 0 0 80
81 56 35 0 0 0 0 0 0 0 0 1 0 0 81
82 97 56 0 0 0 0 0 0 0 0 0 1 0 82
83 210 97 0 0 0 0 0 0 0 0 0 0 1 83
84 260 210 0 0 0 0 0 0 0 0 0 0 0 84
85 257 260 1 0 0 0 0 0 0 0 0 0 0 85
86 210 257 0 1 0 0 0 0 0 0 0 0 0 86
87 125 210 0 0 1 0 0 0 0 0 0 0 0 87
88 80 125 0 0 0 1 0 0 0 0 0 0 0 88
89 42 80 0 0 0 0 1 0 0 0 0 0 0 89
90 35 42 0 0 0 0 0 1 0 0 0 0 0 90
91 31 35 0 0 0 0 0 0 1 0 0 0 0 91
92 32 31 0 0 0 0 0 0 0 1 0 0 0 92
93 50 32 0 0 0 0 0 0 0 0 1 0 0 93
94 92 50 0 0 0 0 0 0 0 0 0 1 0 94
95 189 92 0 0 0 0 0 0 0 0 0 0 1 95
96 256 189 0 0 0 0 0 0 0 0 0 0 0 96
97 250 256 1 0 0 0 0 0 0 0 0 0 0 97
98 198 250 0 1 0 0 0 0 0 0 0 0 0 98
99 136 198 0 0 1 0 0 0 0 0 0 0 0 99
100 73 136 0 0 0 1 0 0 0 0 0 0 0 100
101 39 73 0 0 0 0 1 0 0 0 0 0 0 101
102 32 39 0 0 0 0 0 1 0 0 0 0 0 102
103 30 32 0 0 0 0 0 0 1 0 0 0 0 103
104 31 30 0 0 0 0 0 0 0 1 0 0 0 104
105 45 31 0 0 0 0 0 0 0 0 1 0 0 105
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `verbruik-1` M1 M2 M3
184.1878 0.4119 -35.9915 -81.5347 -110.4643
M4 M5 M6 M7 M8
-152.7013 -158.9428 -152.4180 -154.2632 -149.2004
M9 M10 M11 t
-131.8969 -83.2311 -18.2122 -0.1601
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-35.963 -6.093 0.342 8.963 46.552
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 184.18778 21.75058 8.468 4.08e-13 ***
`verbruik-1` 0.41192 0.09418 4.374 3.25e-05 ***
M1 -35.99148 9.13807 -3.939 0.00016 ***
M2 -81.53474 8.36703 -9.745 8.72e-16 ***
M3 -110.46429 7.42120 -14.885 < 2e-16 ***
M4 -152.70128 9.27941 -16.456 < 2e-16 ***
M5 -158.94278 13.71695 -11.587 < 2e-16 ***
M6 -152.41796 16.35271 -9.321 6.73e-15 ***
M7 -154.26324 16.91922 -9.118 1.79e-14 ***
M8 -149.20036 17.31019 -8.619 1.97e-13 ***
M9 -131.89690 17.03983 -7.741 1.32e-11 ***
M10 -83.23106 15.56331 -5.348 6.57e-07 ***
M11 -18.21223 11.24135 -1.620 0.10867
t -0.16014 0.05617 -2.851 0.00539 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.15 on 91 degrees of freedom
Multiple R-squared: 0.9706, Adjusted R-squared: 0.9664
F-statistic: 231.4 on 13 and 91 DF, p-value: < 2.2e-16
> 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,] 0.6683713 6.632573e-01 3.316287e-01
[2,] 0.5289681 9.420638e-01 4.710319e-01
[3,] 0.3969046 7.938093e-01 6.030954e-01
[4,] 0.2964669 5.929338e-01 7.035331e-01
[5,] 0.2006161 4.012321e-01 7.993839e-01
[6,] 0.1578668 3.157337e-01 8.421332e-01
[7,] 0.3365763 6.731525e-01 6.634237e-01
[8,] 0.4359028 8.718057e-01 5.640972e-01
[9,] 0.9311222 1.377556e-01 6.887780e-02
[10,] 0.9660817 6.783657e-02 3.391829e-02
[11,] 0.9541304 9.173913e-02 4.586957e-02
[12,] 0.9390873 1.218253e-01 6.091267e-02
[13,] 0.9114552 1.770895e-01 8.854475e-02
[14,] 0.8803055 2.393890e-01 1.196945e-01
[15,] 0.8385381 3.229237e-01 1.614619e-01
[16,] 0.8042679 3.914643e-01 1.957321e-01
[17,] 0.8188196 3.623608e-01 1.811804e-01
[18,] 0.8320095 3.359810e-01 1.679905e-01
[19,] 0.8535215 2.929570e-01 1.464785e-01
[20,] 0.9069784 1.860431e-01 9.302156e-02
[21,] 0.9700849 5.983027e-02 2.991513e-02
[22,] 0.9596580 8.068399e-02 4.034200e-02
[23,] 0.9536665 9.266708e-02 4.633354e-02
[24,] 0.9351583 1.296833e-01 6.484165e-02
[25,] 0.9144470 1.711060e-01 8.555298e-02
[26,] 0.9012750 1.974501e-01 9.872503e-02
[27,] 0.8708928 2.582144e-01 1.291072e-01
[28,] 0.8363417 3.273166e-01 1.636583e-01
[29,] 0.8337448 3.325103e-01 1.662552e-01
[30,] 0.8636546 2.726908e-01 1.363454e-01
[31,] 0.8300198 3.399604e-01 1.699802e-01
[32,] 0.8524106 2.951789e-01 1.475894e-01
[33,] 0.8376029 3.247942e-01 1.623971e-01
[34,] 0.8851395 2.297210e-01 1.148605e-01
[35,] 0.9680284 6.394330e-02 3.197165e-02
[36,] 0.9772964 4.540711e-02 2.270356e-02
[37,] 0.9670970 6.580603e-02 3.290302e-02
[38,] 0.9850539 2.989217e-02 1.494608e-02
[39,] 0.9805605 3.887904e-02 1.943952e-02
[40,] 0.9721329 5.573421e-02 2.786711e-02
[41,] 0.9605844 7.883115e-02 3.941558e-02
[42,] 0.9902791 1.944171e-02 9.720853e-03
[43,] 0.9863030 2.739400e-02 1.369700e-02
[44,] 0.9797804 4.043927e-02 2.021963e-02
[45,] 0.9763161 4.736786e-02 2.368393e-02
[46,] 0.9776913 4.461747e-02 2.230874e-02
[47,] 0.9749645 5.007100e-02 2.503550e-02
[48,] 0.9685537 6.289262e-02 3.144631e-02
[49,] 0.9549052 9.018970e-02 4.509485e-02
[50,] 0.9360181 1.279638e-01 6.398188e-02
[51,] 0.9121817 1.756366e-01 8.781829e-02
[52,] 0.8807647 2.384705e-01 1.192353e-01
[53,] 0.8654962 2.690076e-01 1.345038e-01
[54,] 0.9936017 1.279655e-02 6.398273e-03
[55,] 0.9914507 1.709870e-02 8.549349e-03
[56,] 0.9899119 2.017620e-02 1.008810e-02
[57,] 0.9841716 3.165680e-02 1.582840e-02
[58,] 0.9999798 4.037503e-05 2.018752e-05
[59,] 0.9999651 6.984035e-05 3.492018e-05
[60,] 0.9999985 2.983980e-06 1.491990e-06
[61,] 0.9999971 5.855317e-06 2.927659e-06
[62,] 0.9999898 2.048949e-05 1.024475e-05
[63,] 0.9999678 6.430287e-05 3.215143e-05
[64,] 0.9999094 1.812706e-04 9.063532e-05
[65,] 0.9996993 6.013289e-04 3.006645e-04
[66,] 0.9990958 1.808498e-03 9.042492e-04
[67,] 0.9997946 4.107320e-04 2.053660e-04
[68,] 0.9994752 1.049625e-03 5.248125e-04
[69,] 0.9984337 3.132642e-03 1.566321e-03
[70,] 0.9995422 9.156152e-04 4.578076e-04
[71,] 0.9999261 1.478148e-04 7.390742e-05
[72,] 0.9992194 1.561108e-03 7.805542e-04
> postscript(file="/var/wessaorg/rcomp/tmp/1zl0y1324653059.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/2s3kl1324653059.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/3okpw1324653059.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/4l52i1324653059.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/5ygj91324653059.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 105
Frequency = 1
1 2 3 4 5 6
-10.43750589 7.74288247 11.95726876 -2.93283108 11.36317422 -19.52722835
7 8 9 10 11 12
0.48363495 -4.06680930 -2.03398709 24.22180401 -21.82666102 -20.71077741
13 14 15 16 17 18
46.55191679 20.48290706 9.99269318 3.51724295 -10.59867673 -2.95517991
19 20 21 22 23 24
-0.06629102 -4.14518540 2.71148676 24.31957795 19.50688785 -0.85422766
25 26 27 28 29 30
-29.53190858 -8.23227043 2.39404083 10.03239113 0.03027197 -6.09318091
31 32 33 34 35 36
1.09110781 6.60028845 -21.66228918 -0.28452366 -22.74218904 21.37096995
37 38 39 40 41 42
-35.96258477 -6.19139676 8.96336464 -4.28175992 2.13077054 1.29999285
43 44 45 46 47 48
-1.22304323 0.75768728 9.96665952 -21.07294920 -0.87746541 -21.88628083
49 50 51 52 53 54
9.02411328 18.60825272 27.87953927 -18.12981035 -3.71183064 22.92894155
55 56 57 58 59 60
-7.89221899 1.50316113 1.71213337 -32.62015054 -6.54119190 -3.78578226
61 62 63 64 65 66
-9.17351261 -20.87932373 -15.42643788 9.56693713 -2.37828176 0.26249043
67 68 69 70 71 72
2.32752936 1.24863498 10.86953220 33.82992348 14.78145845 13.77536765
73 74 75 76 77 78
2.33073704 -34.84389932 -21.03326412 -10.39218919 5.25339157 3.24373921
79 80 81 82 83 84
0.66107823 -1.18204122 2.26270596 -13.89342791 17.35895615 2.75934107
85 86 87 88 89 90
15.31471046 15.25387462 -21.29596416 11.11478509 -1.94695895 0.34151315
91 92 93 94 95 96
1.23040204 -1.02464239 -0.57989521 -14.50025414 0.34020493 9.33138950
97 98 99 100 101 102
11.88403427 8.05897336 -3.43124052 1.50523423 -0.14186021 0.49891198
103 104 105
3.38780087 0.30890648 -3.24634633
> postscript(file="/var/wessaorg/rcomp/tmp/652f11324653059.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 105
Frequency = 1
lag(myerror, k = 1) myerror
0 -10.43750589 NA
1 7.74288247 -10.43750589
2 11.95726876 7.74288247
3 -2.93283108 11.95726876
4 11.36317422 -2.93283108
5 -19.52722835 11.36317422
6 0.48363495 -19.52722835
7 -4.06680930 0.48363495
8 -2.03398709 -4.06680930
9 24.22180401 -2.03398709
10 -21.82666102 24.22180401
11 -20.71077741 -21.82666102
12 46.55191679 -20.71077741
13 20.48290706 46.55191679
14 9.99269318 20.48290706
15 3.51724295 9.99269318
16 -10.59867673 3.51724295
17 -2.95517991 -10.59867673
18 -0.06629102 -2.95517991
19 -4.14518540 -0.06629102
20 2.71148676 -4.14518540
21 24.31957795 2.71148676
22 19.50688785 24.31957795
23 -0.85422766 19.50688785
24 -29.53190858 -0.85422766
25 -8.23227043 -29.53190858
26 2.39404083 -8.23227043
27 10.03239113 2.39404083
28 0.03027197 10.03239113
29 -6.09318091 0.03027197
30 1.09110781 -6.09318091
31 6.60028845 1.09110781
32 -21.66228918 6.60028845
33 -0.28452366 -21.66228918
34 -22.74218904 -0.28452366
35 21.37096995 -22.74218904
36 -35.96258477 21.37096995
37 -6.19139676 -35.96258477
38 8.96336464 -6.19139676
39 -4.28175992 8.96336464
40 2.13077054 -4.28175992
41 1.29999285 2.13077054
42 -1.22304323 1.29999285
43 0.75768728 -1.22304323
44 9.96665952 0.75768728
45 -21.07294920 9.96665952
46 -0.87746541 -21.07294920
47 -21.88628083 -0.87746541
48 9.02411328 -21.88628083
49 18.60825272 9.02411328
50 27.87953927 18.60825272
51 -18.12981035 27.87953927
52 -3.71183064 -18.12981035
53 22.92894155 -3.71183064
54 -7.89221899 22.92894155
55 1.50316113 -7.89221899
56 1.71213337 1.50316113
57 -32.62015054 1.71213337
58 -6.54119190 -32.62015054
59 -3.78578226 -6.54119190
60 -9.17351261 -3.78578226
61 -20.87932373 -9.17351261
62 -15.42643788 -20.87932373
63 9.56693713 -15.42643788
64 -2.37828176 9.56693713
65 0.26249043 -2.37828176
66 2.32752936 0.26249043
67 1.24863498 2.32752936
68 10.86953220 1.24863498
69 33.82992348 10.86953220
70 14.78145845 33.82992348
71 13.77536765 14.78145845
72 2.33073704 13.77536765
73 -34.84389932 2.33073704
74 -21.03326412 -34.84389932
75 -10.39218919 -21.03326412
76 5.25339157 -10.39218919
77 3.24373921 5.25339157
78 0.66107823 3.24373921
79 -1.18204122 0.66107823
80 2.26270596 -1.18204122
81 -13.89342791 2.26270596
82 17.35895615 -13.89342791
83 2.75934107 17.35895615
84 15.31471046 2.75934107
85 15.25387462 15.31471046
86 -21.29596416 15.25387462
87 11.11478509 -21.29596416
88 -1.94695895 11.11478509
89 0.34151315 -1.94695895
90 1.23040204 0.34151315
91 -1.02464239 1.23040204
92 -0.57989521 -1.02464239
93 -14.50025414 -0.57989521
94 0.34020493 -14.50025414
95 9.33138950 0.34020493
96 11.88403427 9.33138950
97 8.05897336 11.88403427
98 -3.43124052 8.05897336
99 1.50523423 -3.43124052
100 -0.14186021 1.50523423
101 0.49891198 -0.14186021
102 3.38780087 0.49891198
103 0.30890648 3.38780087
104 -3.24634633 0.30890648
105 NA -3.24634633
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 7.74288247 -10.43750589
[2,] 11.95726876 7.74288247
[3,] -2.93283108 11.95726876
[4,] 11.36317422 -2.93283108
[5,] -19.52722835 11.36317422
[6,] 0.48363495 -19.52722835
[7,] -4.06680930 0.48363495
[8,] -2.03398709 -4.06680930
[9,] 24.22180401 -2.03398709
[10,] -21.82666102 24.22180401
[11,] -20.71077741 -21.82666102
[12,] 46.55191679 -20.71077741
[13,] 20.48290706 46.55191679
[14,] 9.99269318 20.48290706
[15,] 3.51724295 9.99269318
[16,] -10.59867673 3.51724295
[17,] -2.95517991 -10.59867673
[18,] -0.06629102 -2.95517991
[19,] -4.14518540 -0.06629102
[20,] 2.71148676 -4.14518540
[21,] 24.31957795 2.71148676
[22,] 19.50688785 24.31957795
[23,] -0.85422766 19.50688785
[24,] -29.53190858 -0.85422766
[25,] -8.23227043 -29.53190858
[26,] 2.39404083 -8.23227043
[27,] 10.03239113 2.39404083
[28,] 0.03027197 10.03239113
[29,] -6.09318091 0.03027197
[30,] 1.09110781 -6.09318091
[31,] 6.60028845 1.09110781
[32,] -21.66228918 6.60028845
[33,] -0.28452366 -21.66228918
[34,] -22.74218904 -0.28452366
[35,] 21.37096995 -22.74218904
[36,] -35.96258477 21.37096995
[37,] -6.19139676 -35.96258477
[38,] 8.96336464 -6.19139676
[39,] -4.28175992 8.96336464
[40,] 2.13077054 -4.28175992
[41,] 1.29999285 2.13077054
[42,] -1.22304323 1.29999285
[43,] 0.75768728 -1.22304323
[44,] 9.96665952 0.75768728
[45,] -21.07294920 9.96665952
[46,] -0.87746541 -21.07294920
[47,] -21.88628083 -0.87746541
[48,] 9.02411328 -21.88628083
[49,] 18.60825272 9.02411328
[50,] 27.87953927 18.60825272
[51,] -18.12981035 27.87953927
[52,] -3.71183064 -18.12981035
[53,] 22.92894155 -3.71183064
[54,] -7.89221899 22.92894155
[55,] 1.50316113 -7.89221899
[56,] 1.71213337 1.50316113
[57,] -32.62015054 1.71213337
[58,] -6.54119190 -32.62015054
[59,] -3.78578226 -6.54119190
[60,] -9.17351261 -3.78578226
[61,] -20.87932373 -9.17351261
[62,] -15.42643788 -20.87932373
[63,] 9.56693713 -15.42643788
[64,] -2.37828176 9.56693713
[65,] 0.26249043 -2.37828176
[66,] 2.32752936 0.26249043
[67,] 1.24863498 2.32752936
[68,] 10.86953220 1.24863498
[69,] 33.82992348 10.86953220
[70,] 14.78145845 33.82992348
[71,] 13.77536765 14.78145845
[72,] 2.33073704 13.77536765
[73,] -34.84389932 2.33073704
[74,] -21.03326412 -34.84389932
[75,] -10.39218919 -21.03326412
[76,] 5.25339157 -10.39218919
[77,] 3.24373921 5.25339157
[78,] 0.66107823 3.24373921
[79,] -1.18204122 0.66107823
[80,] 2.26270596 -1.18204122
[81,] -13.89342791 2.26270596
[82,] 17.35895615 -13.89342791
[83,] 2.75934107 17.35895615
[84,] 15.31471046 2.75934107
[85,] 15.25387462 15.31471046
[86,] -21.29596416 15.25387462
[87,] 11.11478509 -21.29596416
[88,] -1.94695895 11.11478509
[89,] 0.34151315 -1.94695895
[90,] 1.23040204 0.34151315
[91,] -1.02464239 1.23040204
[92,] -0.57989521 -1.02464239
[93,] -14.50025414 -0.57989521
[94,] 0.34020493 -14.50025414
[95,] 9.33138950 0.34020493
[96,] 11.88403427 9.33138950
[97,] 8.05897336 11.88403427
[98,] -3.43124052 8.05897336
[99,] 1.50523423 -3.43124052
[100,] -0.14186021 1.50523423
[101,] 0.49891198 -0.14186021
[102,] 3.38780087 0.49891198
[103,] 0.30890648 3.38780087
[104,] -3.24634633 0.30890648
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 7.74288247 -10.43750589
2 11.95726876 7.74288247
3 -2.93283108 11.95726876
4 11.36317422 -2.93283108
5 -19.52722835 11.36317422
6 0.48363495 -19.52722835
7 -4.06680930 0.48363495
8 -2.03398709 -4.06680930
9 24.22180401 -2.03398709
10 -21.82666102 24.22180401
11 -20.71077741 -21.82666102
12 46.55191679 -20.71077741
13 20.48290706 46.55191679
14 9.99269318 20.48290706
15 3.51724295 9.99269318
16 -10.59867673 3.51724295
17 -2.95517991 -10.59867673
18 -0.06629102 -2.95517991
19 -4.14518540 -0.06629102
20 2.71148676 -4.14518540
21 24.31957795 2.71148676
22 19.50688785 24.31957795
23 -0.85422766 19.50688785
24 -29.53190858 -0.85422766
25 -8.23227043 -29.53190858
26 2.39404083 -8.23227043
27 10.03239113 2.39404083
28 0.03027197 10.03239113
29 -6.09318091 0.03027197
30 1.09110781 -6.09318091
31 6.60028845 1.09110781
32 -21.66228918 6.60028845
33 -0.28452366 -21.66228918
34 -22.74218904 -0.28452366
35 21.37096995 -22.74218904
36 -35.96258477 21.37096995
37 -6.19139676 -35.96258477
38 8.96336464 -6.19139676
39 -4.28175992 8.96336464
40 2.13077054 -4.28175992
41 1.29999285 2.13077054
42 -1.22304323 1.29999285
43 0.75768728 -1.22304323
44 9.96665952 0.75768728
45 -21.07294920 9.96665952
46 -0.87746541 -21.07294920
47 -21.88628083 -0.87746541
48 9.02411328 -21.88628083
49 18.60825272 9.02411328
50 27.87953927 18.60825272
51 -18.12981035 27.87953927
52 -3.71183064 -18.12981035
53 22.92894155 -3.71183064
54 -7.89221899 22.92894155
55 1.50316113 -7.89221899
56 1.71213337 1.50316113
57 -32.62015054 1.71213337
58 -6.54119190 -32.62015054
59 -3.78578226 -6.54119190
60 -9.17351261 -3.78578226
61 -20.87932373 -9.17351261
62 -15.42643788 -20.87932373
63 9.56693713 -15.42643788
64 -2.37828176 9.56693713
65 0.26249043 -2.37828176
66 2.32752936 0.26249043
67 1.24863498 2.32752936
68 10.86953220 1.24863498
69 33.82992348 10.86953220
70 14.78145845 33.82992348
71 13.77536765 14.78145845
72 2.33073704 13.77536765
73 -34.84389932 2.33073704
74 -21.03326412 -34.84389932
75 -10.39218919 -21.03326412
76 5.25339157 -10.39218919
77 3.24373921 5.25339157
78 0.66107823 3.24373921
79 -1.18204122 0.66107823
80 2.26270596 -1.18204122
81 -13.89342791 2.26270596
82 17.35895615 -13.89342791
83 2.75934107 17.35895615
84 15.31471046 2.75934107
85 15.25387462 15.31471046
86 -21.29596416 15.25387462
87 11.11478509 -21.29596416
88 -1.94695895 11.11478509
89 0.34151315 -1.94695895
90 1.23040204 0.34151315
91 -1.02464239 1.23040204
92 -0.57989521 -1.02464239
93 -14.50025414 -0.57989521
94 0.34020493 -14.50025414
95 9.33138950 0.34020493
96 11.88403427 9.33138950
97 8.05897336 11.88403427
98 -3.43124052 8.05897336
99 1.50523423 -3.43124052
100 -0.14186021 1.50523423
101 0.49891198 -0.14186021
102 3.38780087 0.49891198
103 0.30890648 3.38780087
104 -3.24634633 0.30890648
> 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/wessaorg/rcomp/tmp/798e51324653059.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/8l4rl1324653059.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/9pski1324653059.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/wessaorg/rcomp/tmp/10v44k1324653059.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11pthb1324653059.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/wessaorg/rcomp/tmp/12yjux1324653059.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/wessaorg/rcomp/tmp/130d651324653059.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/wessaorg/rcomp/tmp/14fl7p1324653059.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/wessaorg/rcomp/tmp/15umvg1324653059.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/wessaorg/rcomp/tmp/16qm1f1324653059.tab")
+ }
>
> try(system("convert tmp/1zl0y1324653059.ps tmp/1zl0y1324653059.png",intern=TRUE))
character(0)
> try(system("convert tmp/2s3kl1324653059.ps tmp/2s3kl1324653059.png",intern=TRUE))
character(0)
> try(system("convert tmp/3okpw1324653059.ps tmp/3okpw1324653059.png",intern=TRUE))
character(0)
> try(system("convert tmp/4l52i1324653059.ps tmp/4l52i1324653059.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ygj91324653059.ps tmp/5ygj91324653059.png",intern=TRUE))
character(0)
> try(system("convert tmp/652f11324653059.ps tmp/652f11324653059.png",intern=TRUE))
character(0)
> try(system("convert tmp/798e51324653059.ps tmp/798e51324653059.png",intern=TRUE))
character(0)
> try(system("convert tmp/8l4rl1324653059.ps tmp/8l4rl1324653059.png",intern=TRUE))
character(0)
> try(system("convert tmp/9pski1324653059.ps tmp/9pski1324653059.png",intern=TRUE))
character(0)
> try(system("convert tmp/10v44k1324653059.ps tmp/10v44k1324653059.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.794 0.610 4.412