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(302,262,218,175,100,77,43,47,49,69,152,205,246,294,242,181,107,56,49,47,47,71,151,244,280,230,185,148,98,61,46,45,55,48,115,185,276,220,181,151,83,55,49,42,46,74,103,200,237,247,215,182,80,46,65,40,44,63,85,185,247,231,167,117,79,45,40,38,41,69,152,232,282,255,161,107,53,40,39,34,35,56,97,210,260,257,210,125,80,42,35,31,32,50,92,189,256,250,198,136,73,39,32,30,31,45),dim=c(1,106),dimnames=list(c('Gasverbruik'),1:106))
> y <- array(NA,dim=c(1,106),dimnames=list(c('Gasverbruik'),1:106))
> 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
> 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 302 1 0 0 0 0 0 0 0 0 0 0 1
2 262 0 1 0 0 0 0 0 0 0 0 0 2
3 218 0 0 1 0 0 0 0 0 0 0 0 3
4 175 0 0 0 1 0 0 0 0 0 0 0 4
5 100 0 0 0 0 1 0 0 0 0 0 0 5
6 77 0 0 0 0 0 1 0 0 0 0 0 6
7 43 0 0 0 0 0 0 1 0 0 0 0 7
8 47 0 0 0 0 0 0 0 1 0 0 0 8
9 49 0 0 0 0 0 0 0 0 1 0 0 9
10 69 0 0 0 0 0 0 0 0 0 1 0 10
11 152 0 0 0 0 0 0 0 0 0 0 1 11
12 205 0 0 0 0 0 0 0 0 0 0 0 12
13 246 1 0 0 0 0 0 0 0 0 0 0 13
14 294 0 1 0 0 0 0 0 0 0 0 0 14
15 242 0 0 1 0 0 0 0 0 0 0 0 15
16 181 0 0 0 1 0 0 0 0 0 0 0 16
17 107 0 0 0 0 1 0 0 0 0 0 0 17
18 56 0 0 0 0 0 1 0 0 0 0 0 18
19 49 0 0 0 0 0 0 1 0 0 0 0 19
20 47 0 0 0 0 0 0 0 1 0 0 0 20
21 47 0 0 0 0 0 0 0 0 1 0 0 21
22 71 0 0 0 0 0 0 0 0 0 1 0 22
23 151 0 0 0 0 0 0 0 0 0 0 1 23
24 244 0 0 0 0 0 0 0 0 0 0 0 24
25 280 1 0 0 0 0 0 0 0 0 0 0 25
26 230 0 1 0 0 0 0 0 0 0 0 0 26
27 185 0 0 1 0 0 0 0 0 0 0 0 27
28 148 0 0 0 1 0 0 0 0 0 0 0 28
29 98 0 0 0 0 1 0 0 0 0 0 0 29
30 61 0 0 0 0 0 1 0 0 0 0 0 30
31 46 0 0 0 0 0 0 1 0 0 0 0 31
32 45 0 0 0 0 0 0 0 1 0 0 0 32
33 55 0 0 0 0 0 0 0 0 1 0 0 33
34 48 0 0 0 0 0 0 0 0 0 1 0 34
35 115 0 0 0 0 0 0 0 0 0 0 1 35
36 185 0 0 0 0 0 0 0 0 0 0 0 36
37 276 1 0 0 0 0 0 0 0 0 0 0 37
38 220 0 1 0 0 0 0 0 0 0 0 0 38
39 181 0 0 1 0 0 0 0 0 0 0 0 39
40 151 0 0 0 1 0 0 0 0 0 0 0 40
41 83 0 0 0 0 1 0 0 0 0 0 0 41
42 55 0 0 0 0 0 1 0 0 0 0 0 42
43 49 0 0 0 0 0 0 1 0 0 0 0 43
44 42 0 0 0 0 0 0 0 1 0 0 0 44
45 46 0 0 0 0 0 0 0 0 1 0 0 45
46 74 0 0 0 0 0 0 0 0 0 1 0 46
47 103 0 0 0 0 0 0 0 0 0 0 1 47
48 200 0 0 0 0 0 0 0 0 0 0 0 48
49 237 1 0 0 0 0 0 0 0 0 0 0 49
50 247 0 1 0 0 0 0 0 0 0 0 0 50
51 215 0 0 1 0 0 0 0 0 0 0 0 51
52 182 0 0 0 1 0 0 0 0 0 0 0 52
53 80 0 0 0 0 1 0 0 0 0 0 0 53
54 46 0 0 0 0 0 1 0 0 0 0 0 54
55 65 0 0 0 0 0 0 1 0 0 0 0 55
56 40 0 0 0 0 0 0 0 1 0 0 0 56
57 44 0 0 0 0 0 0 0 0 1 0 0 57
58 63 0 0 0 0 0 0 0 0 0 1 0 58
59 85 0 0 0 0 0 0 0 0 0 0 1 59
60 185 0 0 0 0 0 0 0 0 0 0 0 60
61 247 1 0 0 0 0 0 0 0 0 0 0 61
62 231 0 1 0 0 0 0 0 0 0 0 0 62
63 167 0 0 1 0 0 0 0 0 0 0 0 63
64 117 0 0 0 1 0 0 0 0 0 0 0 64
65 79 0 0 0 0 1 0 0 0 0 0 0 65
66 45 0 0 0 0 0 1 0 0 0 0 0 66
67 40 0 0 0 0 0 0 1 0 0 0 0 67
68 38 0 0 0 0 0 0 0 1 0 0 0 68
69 41 0 0 0 0 0 0 0 0 1 0 0 69
70 69 0 0 0 0 0 0 0 0 0 1 0 70
71 152 0 0 0 0 0 0 0 0 0 0 1 71
72 232 0 0 0 0 0 0 0 0 0 0 0 72
73 282 1 0 0 0 0 0 0 0 0 0 0 73
74 255 0 1 0 0 0 0 0 0 0 0 0 74
75 161 0 0 1 0 0 0 0 0 0 0 0 75
76 107 0 0 0 1 0 0 0 0 0 0 0 76
77 53 0 0 0 0 1 0 0 0 0 0 0 77
78 40 0 0 0 0 0 1 0 0 0 0 0 78
79 39 0 0 0 0 0 0 1 0 0 0 0 79
80 34 0 0 0 0 0 0 0 1 0 0 0 80
81 35 0 0 0 0 0 0 0 0 1 0 0 81
82 56 0 0 0 0 0 0 0 0 0 1 0 82
83 97 0 0 0 0 0 0 0 0 0 0 1 83
84 210 0 0 0 0 0 0 0 0 0 0 0 84
85 260 1 0 0 0 0 0 0 0 0 0 0 85
86 257 0 1 0 0 0 0 0 0 0 0 0 86
87 210 0 0 1 0 0 0 0 0 0 0 0 87
88 125 0 0 0 1 0 0 0 0 0 0 0 88
89 80 0 0 0 0 1 0 0 0 0 0 0 89
90 42 0 0 0 0 0 1 0 0 0 0 0 90
91 35 0 0 0 0 0 0 1 0 0 0 0 91
92 31 0 0 0 0 0 0 0 1 0 0 0 92
93 32 0 0 0 0 0 0 0 0 1 0 0 93
94 50 0 0 0 0 0 0 0 0 0 1 0 94
95 92 0 0 0 0 0 0 0 0 0 0 1 95
96 189 0 0 0 0 0 0 0 0 0 0 0 96
97 256 1 0 0 0 0 0 0 0 0 0 0 97
98 250 0 1 0 0 0 0 0 0 0 0 0 98
99 198 0 0 1 0 0 0 0 0 0 0 0 99
100 136 0 0 0 1 0 0 0 0 0 0 0 100
101 73 0 0 0 0 1 0 0 0 0 0 0 101
102 39 0 0 0 0 0 1 0 0 0 0 0 102
103 32 0 0 0 0 0 0 1 0 0 0 0 103
104 30 0 0 0 0 0 0 0 1 0 0 0 104
105 31 0 0 0 0 0 0 0 0 1 0 0 105
106 45 0 0 0 0 0 0 0 0 0 1 0 106
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
222.1020 57.3933 42.1313 -9.6862 -59.9482 -122.8769
M6 M7 M8 M9 M10 M11
-155.0278 -161.7342 -166.3296 -163.1471 -144.5202 -88.1686
t
-0.2936
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-33.078 -6.302 0.467 7.191 38.909
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 222.10197 6.56058 33.854 < 2e-16 ***
M1 57.39334 8.11188 7.075 2.75e-10 ***
M2 42.13134 8.11031 5.195 1.20e-06 ***
M3 -9.68622 8.10909 -1.194 0.235
M4 -59.94822 8.10822 -7.394 6.14e-11 ***
M5 -122.87689 8.10770 -15.156 < 2e-16 ***
M6 -155.02778 8.10752 -19.121 < 2e-16 ***
M7 -161.73422 8.10770 -19.948 < 2e-16 ***
M8 -166.32956 8.10822 -20.514 < 2e-16 ***
M9 -163.14711 8.10909 -20.119 < 2e-16 ***
M10 -144.52022 8.11031 -17.819 < 2e-16 ***
M11 -88.16856 8.34274 -10.568 < 2e-16 ***
t -0.29356 0.05316 -5.522 3.03e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 16.69 on 93 degrees of freedom
Multiple R-squared: 0.9652, Adjusted R-squared: 0.9607
F-statistic: 214.8 on 12 and 93 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.9718529 5.629421e-02 2.814710e-02
[2,] 0.9440383 1.119234e-01 5.596169e-02
[3,] 0.9243393 1.513214e-01 7.566068e-02
[4,] 0.8757141 2.485719e-01 1.242859e-01
[5,] 0.8080736 3.838528e-01 1.919264e-01
[6,] 0.7267144 5.465713e-01 2.732856e-01
[7,] 0.6355239 7.289521e-01 3.644761e-01
[8,] 0.5870809 8.258382e-01 4.129191e-01
[9,] 0.7423193 5.153615e-01 2.576807e-01
[10,] 0.6736250 6.527500e-01 3.263750e-01
[11,] 0.8861243 2.277513e-01 1.138757e-01
[12,] 0.9308737 1.382527e-01 6.912633e-02
[13,] 0.9182308 1.635383e-01 8.176916e-02
[14,] 0.8949666 2.100669e-01 1.050334e-01
[15,] 0.8619454 2.761093e-01 1.380546e-01
[16,] 0.8250291 3.499419e-01 1.749709e-01
[17,] 0.7787858 4.424284e-01 2.212142e-01
[18,] 0.7510067 4.979865e-01 2.489933e-01
[19,] 0.7239860 5.520280e-01 2.760140e-01
[20,] 0.7313170 5.373659e-01 2.686830e-01
[21,] 0.7689340 4.621320e-01 2.310660e-01
[22,] 0.7482157 5.035686e-01 2.517843e-01
[23,] 0.8117266 3.765468e-01 1.882734e-01
[24,] 0.7964941 4.070118e-01 2.035059e-01
[25,] 0.7528312 4.943377e-01 2.471688e-01
[26,] 0.6979609 6.040782e-01 3.020391e-01
[27,] 0.6485884 7.028233e-01 3.514116e-01
[28,] 0.6291438 7.417125e-01 3.708562e-01
[29,] 0.5809547 8.380906e-01 4.190453e-01
[30,] 0.5305664 9.388672e-01 4.694336e-01
[31,] 0.5438677 9.122645e-01 4.561323e-01
[32,] 0.5300272 9.399456e-01 4.699728e-01
[33,] 0.4732522 9.465043e-01 5.267478e-01
[34,] 0.5189814 9.620371e-01 4.810186e-01
[35,] 0.4827566 9.655132e-01 5.172434e-01
[36,] 0.5361332 9.277336e-01 4.638668e-01
[37,] 0.8313686 3.372628e-01 1.686314e-01
[38,] 0.7912421 4.175158e-01 2.087579e-01
[39,] 0.7428618 5.142764e-01 2.571382e-01
[40,] 0.8195286 3.609429e-01 1.804714e-01
[41,] 0.7814993 4.370014e-01 2.185007e-01
[42,] 0.7433487 5.133027e-01 2.566513e-01
[43,] 0.7035963 5.928074e-01 2.964037e-01
[44,] 0.7861881 4.276239e-01 2.138119e-01
[45,] 0.8044106 3.911787e-01 1.955894e-01
[46,] 0.8001949 3.996101e-01 1.998051e-01
[47,] 0.8212847 3.574306e-01 1.787153e-01
[48,] 0.8722504 2.554993e-01 1.277496e-01
[49,] 0.8823988 2.352024e-01 1.176012e-01
[50,] 0.8472511 3.054979e-01 1.527489e-01
[51,] 0.8052820 3.894360e-01 1.947180e-01
[52,] 0.7590500 4.819000e-01 2.409500e-01
[53,] 0.7104058 5.791884e-01 2.895942e-01
[54,] 0.6547171 6.905658e-01 3.452829e-01
[55,] 0.6339512 7.320975e-01 3.660488e-01
[56,] 0.9478256 1.043488e-01 5.217442e-02
[57,] 0.9848398 3.032047e-02 1.516024e-02
[58,] 0.9931478 1.370449e-02 6.852247e-03
[59,] 0.9897631 2.047382e-02 1.023691e-02
[60,] 0.9997376 5.248727e-04 2.624363e-04
[61,] 0.9999724 5.513262e-05 2.756631e-05
[62,] 0.9999999 2.225943e-07 1.112972e-07
[63,] 0.9999997 5.230616e-07 2.615308e-07
[64,] 0.9999989 2.115800e-06 1.057900e-06
[65,] 0.9999963 7.365270e-06 3.682635e-06
[66,] 0.9999879 2.411825e-05 1.205912e-05
[67,] 0.9999559 8.819457e-05 4.409728e-05
[68,] 0.9998469 3.062927e-04 1.531464e-04
[69,] 0.9999659 6.814613e-05 3.407307e-05
[70,] 0.9998526 2.947893e-04 1.473946e-04
[71,] 0.9995251 9.497364e-04 4.748682e-04
[72,] 0.9994406 1.118840e-03 5.594199e-04
[73,] 0.9999810 3.791815e-05 1.895908e-05
[74,] 0.9999629 7.411229e-05 3.705615e-05
[75,] 0.9994103 1.179495e-03 5.897473e-04
> postscript(file="/var/wessaorg/rcomp/tmp/1d63n1324646551.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/22fut1324646551.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/3ad691324646551.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/4rb0q1324646551.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/5ss821324646551.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 = 106
Frequency = 1
1 2 3 4 5 6
22.7982456 -1.6461988 6.4649123 14.0204678 2.2426901 11.6871345
7 8 9 10 11 12
-15.3128655 -6.4239766 -7.3128655 -5.6461988 21.2956871 -13.5793129
13 14 15 16 17 18
-29.6790936 33.8764620 33.9875731 23.5431287 12.7653509 -5.7902047
19 20 21 22 23 24
-5.7902047 -2.9013158 -5.7902047 -0.1235380 23.8183480 28.9433480
25 26 27 28 29 30
7.8435673 -26.6008772 -19.4897661 -5.9342105 7.2880117 2.7324561
31 32 33 34 35 36
-5.2675439 -1.3786550 5.7324561 -19.6008772 -8.6589912 -26.5339912
37 38 39 40 41 42
7.3662281 -33.0782164 -19.9671053 0.5884503 -4.1893275 0.2551170
43 44 45 46 47 48
1.2551170 -0.8559942 0.2551170 9.9217836 -17.1363304 -8.0113304
49 50 51 52 53 54
-28.1111111 -2.5555556 17.5555556 35.1111111 -3.6666667 -5.2222222
55 56 57 58 59 60
20.7777778 0.6666667 1.7777778 2.4444444 -31.6136696 -19.4886696
61 62 63 64 65 66
-14.5884503 -15.0328947 -26.9217836 -26.3662281 -1.1440058 -2.6995614
67 68 69 70 71 72
-0.6995614 2.1893275 2.3004386 11.9671053 38.9089912 31.0339912
73 74 75 76 77 78
23.9342105 12.4897661 -29.3991228 -32.8435673 -23.6213450 -4.1769006
79 80 81 82 83 84
1.8230994 1.7119883 -0.1769006 2.4897661 -12.5683480 12.5566520
85 86 87 88 89 90
5.4568713 18.0124269 23.1235380 -11.3209064 6.9013158 1.3457602
91 92 93 94 95 96
1.3457602 2.2346491 0.3457602 0.0124269 -14.0456871 -4.9206871
97 98 99 100 101 102
4.9795322 14.5350877 14.6461988 3.2017544 3.4239766 1.8684211
103 104 105 106
1.8684211 4.7573099 2.8684211 -1.4649123
> postscript(file="/var/wessaorg/rcomp/tmp/63pmz1324646551.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 = 106
Frequency = 1
lag(myerror, k = 1) myerror
0 22.7982456 NA
1 -1.6461988 22.7982456
2 6.4649123 -1.6461988
3 14.0204678 6.4649123
4 2.2426901 14.0204678
5 11.6871345 2.2426901
6 -15.3128655 11.6871345
7 -6.4239766 -15.3128655
8 -7.3128655 -6.4239766
9 -5.6461988 -7.3128655
10 21.2956871 -5.6461988
11 -13.5793129 21.2956871
12 -29.6790936 -13.5793129
13 33.8764620 -29.6790936
14 33.9875731 33.8764620
15 23.5431287 33.9875731
16 12.7653509 23.5431287
17 -5.7902047 12.7653509
18 -5.7902047 -5.7902047
19 -2.9013158 -5.7902047
20 -5.7902047 -2.9013158
21 -0.1235380 -5.7902047
22 23.8183480 -0.1235380
23 28.9433480 23.8183480
24 7.8435673 28.9433480
25 -26.6008772 7.8435673
26 -19.4897661 -26.6008772
27 -5.9342105 -19.4897661
28 7.2880117 -5.9342105
29 2.7324561 7.2880117
30 -5.2675439 2.7324561
31 -1.3786550 -5.2675439
32 5.7324561 -1.3786550
33 -19.6008772 5.7324561
34 -8.6589912 -19.6008772
35 -26.5339912 -8.6589912
36 7.3662281 -26.5339912
37 -33.0782164 7.3662281
38 -19.9671053 -33.0782164
39 0.5884503 -19.9671053
40 -4.1893275 0.5884503
41 0.2551170 -4.1893275
42 1.2551170 0.2551170
43 -0.8559942 1.2551170
44 0.2551170 -0.8559942
45 9.9217836 0.2551170
46 -17.1363304 9.9217836
47 -8.0113304 -17.1363304
48 -28.1111111 -8.0113304
49 -2.5555556 -28.1111111
50 17.5555556 -2.5555556
51 35.1111111 17.5555556
52 -3.6666667 35.1111111
53 -5.2222222 -3.6666667
54 20.7777778 -5.2222222
55 0.6666667 20.7777778
56 1.7777778 0.6666667
57 2.4444444 1.7777778
58 -31.6136696 2.4444444
59 -19.4886696 -31.6136696
60 -14.5884503 -19.4886696
61 -15.0328947 -14.5884503
62 -26.9217836 -15.0328947
63 -26.3662281 -26.9217836
64 -1.1440058 -26.3662281
65 -2.6995614 -1.1440058
66 -0.6995614 -2.6995614
67 2.1893275 -0.6995614
68 2.3004386 2.1893275
69 11.9671053 2.3004386
70 38.9089912 11.9671053
71 31.0339912 38.9089912
72 23.9342105 31.0339912
73 12.4897661 23.9342105
74 -29.3991228 12.4897661
75 -32.8435673 -29.3991228
76 -23.6213450 -32.8435673
77 -4.1769006 -23.6213450
78 1.8230994 -4.1769006
79 1.7119883 1.8230994
80 -0.1769006 1.7119883
81 2.4897661 -0.1769006
82 -12.5683480 2.4897661
83 12.5566520 -12.5683480
84 5.4568713 12.5566520
85 18.0124269 5.4568713
86 23.1235380 18.0124269
87 -11.3209064 23.1235380
88 6.9013158 -11.3209064
89 1.3457602 6.9013158
90 1.3457602 1.3457602
91 2.2346491 1.3457602
92 0.3457602 2.2346491
93 0.0124269 0.3457602
94 -14.0456871 0.0124269
95 -4.9206871 -14.0456871
96 4.9795322 -4.9206871
97 14.5350877 4.9795322
98 14.6461988 14.5350877
99 3.2017544 14.6461988
100 3.4239766 3.2017544
101 1.8684211 3.4239766
102 1.8684211 1.8684211
103 4.7573099 1.8684211
104 2.8684211 4.7573099
105 -1.4649123 2.8684211
106 NA -1.4649123
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.6461988 22.7982456
[2,] 6.4649123 -1.6461988
[3,] 14.0204678 6.4649123
[4,] 2.2426901 14.0204678
[5,] 11.6871345 2.2426901
[6,] -15.3128655 11.6871345
[7,] -6.4239766 -15.3128655
[8,] -7.3128655 -6.4239766
[9,] -5.6461988 -7.3128655
[10,] 21.2956871 -5.6461988
[11,] -13.5793129 21.2956871
[12,] -29.6790936 -13.5793129
[13,] 33.8764620 -29.6790936
[14,] 33.9875731 33.8764620
[15,] 23.5431287 33.9875731
[16,] 12.7653509 23.5431287
[17,] -5.7902047 12.7653509
[18,] -5.7902047 -5.7902047
[19,] -2.9013158 -5.7902047
[20,] -5.7902047 -2.9013158
[21,] -0.1235380 -5.7902047
[22,] 23.8183480 -0.1235380
[23,] 28.9433480 23.8183480
[24,] 7.8435673 28.9433480
[25,] -26.6008772 7.8435673
[26,] -19.4897661 -26.6008772
[27,] -5.9342105 -19.4897661
[28,] 7.2880117 -5.9342105
[29,] 2.7324561 7.2880117
[30,] -5.2675439 2.7324561
[31,] -1.3786550 -5.2675439
[32,] 5.7324561 -1.3786550
[33,] -19.6008772 5.7324561
[34,] -8.6589912 -19.6008772
[35,] -26.5339912 -8.6589912
[36,] 7.3662281 -26.5339912
[37,] -33.0782164 7.3662281
[38,] -19.9671053 -33.0782164
[39,] 0.5884503 -19.9671053
[40,] -4.1893275 0.5884503
[41,] 0.2551170 -4.1893275
[42,] 1.2551170 0.2551170
[43,] -0.8559942 1.2551170
[44,] 0.2551170 -0.8559942
[45,] 9.9217836 0.2551170
[46,] -17.1363304 9.9217836
[47,] -8.0113304 -17.1363304
[48,] -28.1111111 -8.0113304
[49,] -2.5555556 -28.1111111
[50,] 17.5555556 -2.5555556
[51,] 35.1111111 17.5555556
[52,] -3.6666667 35.1111111
[53,] -5.2222222 -3.6666667
[54,] 20.7777778 -5.2222222
[55,] 0.6666667 20.7777778
[56,] 1.7777778 0.6666667
[57,] 2.4444444 1.7777778
[58,] -31.6136696 2.4444444
[59,] -19.4886696 -31.6136696
[60,] -14.5884503 -19.4886696
[61,] -15.0328947 -14.5884503
[62,] -26.9217836 -15.0328947
[63,] -26.3662281 -26.9217836
[64,] -1.1440058 -26.3662281
[65,] -2.6995614 -1.1440058
[66,] -0.6995614 -2.6995614
[67,] 2.1893275 -0.6995614
[68,] 2.3004386 2.1893275
[69,] 11.9671053 2.3004386
[70,] 38.9089912 11.9671053
[71,] 31.0339912 38.9089912
[72,] 23.9342105 31.0339912
[73,] 12.4897661 23.9342105
[74,] -29.3991228 12.4897661
[75,] -32.8435673 -29.3991228
[76,] -23.6213450 -32.8435673
[77,] -4.1769006 -23.6213450
[78,] 1.8230994 -4.1769006
[79,] 1.7119883 1.8230994
[80,] -0.1769006 1.7119883
[81,] 2.4897661 -0.1769006
[82,] -12.5683480 2.4897661
[83,] 12.5566520 -12.5683480
[84,] 5.4568713 12.5566520
[85,] 18.0124269 5.4568713
[86,] 23.1235380 18.0124269
[87,] -11.3209064 23.1235380
[88,] 6.9013158 -11.3209064
[89,] 1.3457602 6.9013158
[90,] 1.3457602 1.3457602
[91,] 2.2346491 1.3457602
[92,] 0.3457602 2.2346491
[93,] 0.0124269 0.3457602
[94,] -14.0456871 0.0124269
[95,] -4.9206871 -14.0456871
[96,] 4.9795322 -4.9206871
[97,] 14.5350877 4.9795322
[98,] 14.6461988 14.5350877
[99,] 3.2017544 14.6461988
[100,] 3.4239766 3.2017544
[101,] 1.8684211 3.4239766
[102,] 1.8684211 1.8684211
[103,] 4.7573099 1.8684211
[104,] 2.8684211 4.7573099
[105,] -1.4649123 2.8684211
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.6461988 22.7982456
2 6.4649123 -1.6461988
3 14.0204678 6.4649123
4 2.2426901 14.0204678
5 11.6871345 2.2426901
6 -15.3128655 11.6871345
7 -6.4239766 -15.3128655
8 -7.3128655 -6.4239766
9 -5.6461988 -7.3128655
10 21.2956871 -5.6461988
11 -13.5793129 21.2956871
12 -29.6790936 -13.5793129
13 33.8764620 -29.6790936
14 33.9875731 33.8764620
15 23.5431287 33.9875731
16 12.7653509 23.5431287
17 -5.7902047 12.7653509
18 -5.7902047 -5.7902047
19 -2.9013158 -5.7902047
20 -5.7902047 -2.9013158
21 -0.1235380 -5.7902047
22 23.8183480 -0.1235380
23 28.9433480 23.8183480
24 7.8435673 28.9433480
25 -26.6008772 7.8435673
26 -19.4897661 -26.6008772
27 -5.9342105 -19.4897661
28 7.2880117 -5.9342105
29 2.7324561 7.2880117
30 -5.2675439 2.7324561
31 -1.3786550 -5.2675439
32 5.7324561 -1.3786550
33 -19.6008772 5.7324561
34 -8.6589912 -19.6008772
35 -26.5339912 -8.6589912
36 7.3662281 -26.5339912
37 -33.0782164 7.3662281
38 -19.9671053 -33.0782164
39 0.5884503 -19.9671053
40 -4.1893275 0.5884503
41 0.2551170 -4.1893275
42 1.2551170 0.2551170
43 -0.8559942 1.2551170
44 0.2551170 -0.8559942
45 9.9217836 0.2551170
46 -17.1363304 9.9217836
47 -8.0113304 -17.1363304
48 -28.1111111 -8.0113304
49 -2.5555556 -28.1111111
50 17.5555556 -2.5555556
51 35.1111111 17.5555556
52 -3.6666667 35.1111111
53 -5.2222222 -3.6666667
54 20.7777778 -5.2222222
55 0.6666667 20.7777778
56 1.7777778 0.6666667
57 2.4444444 1.7777778
58 -31.6136696 2.4444444
59 -19.4886696 -31.6136696
60 -14.5884503 -19.4886696
61 -15.0328947 -14.5884503
62 -26.9217836 -15.0328947
63 -26.3662281 -26.9217836
64 -1.1440058 -26.3662281
65 -2.6995614 -1.1440058
66 -0.6995614 -2.6995614
67 2.1893275 -0.6995614
68 2.3004386 2.1893275
69 11.9671053 2.3004386
70 38.9089912 11.9671053
71 31.0339912 38.9089912
72 23.9342105 31.0339912
73 12.4897661 23.9342105
74 -29.3991228 12.4897661
75 -32.8435673 -29.3991228
76 -23.6213450 -32.8435673
77 -4.1769006 -23.6213450
78 1.8230994 -4.1769006
79 1.7119883 1.8230994
80 -0.1769006 1.7119883
81 2.4897661 -0.1769006
82 -12.5683480 2.4897661
83 12.5566520 -12.5683480
84 5.4568713 12.5566520
85 18.0124269 5.4568713
86 23.1235380 18.0124269
87 -11.3209064 23.1235380
88 6.9013158 -11.3209064
89 1.3457602 6.9013158
90 1.3457602 1.3457602
91 2.2346491 1.3457602
92 0.3457602 2.2346491
93 0.0124269 0.3457602
94 -14.0456871 0.0124269
95 -4.9206871 -14.0456871
96 4.9795322 -4.9206871
97 14.5350877 4.9795322
98 14.6461988 14.5350877
99 3.2017544 14.6461988
100 3.4239766 3.2017544
101 1.8684211 3.4239766
102 1.8684211 1.8684211
103 4.7573099 1.8684211
104 2.8684211 4.7573099
105 -1.4649123 2.8684211
> 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/7q4wp1324646551.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/8sm161324646551.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/9t5tv1324646551.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/10kopv1324646551.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/111w5t1324646551.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/1275j01324646551.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/13ys311324646551.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/14hp071324646551.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/15txxm1324646551.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/16jysd1324646551.tab")
+ }
>
> try(system("convert tmp/1d63n1324646551.ps tmp/1d63n1324646551.png",intern=TRUE))
character(0)
> try(system("convert tmp/22fut1324646551.ps tmp/22fut1324646551.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ad691324646551.ps tmp/3ad691324646551.png",intern=TRUE))
character(0)
> try(system("convert tmp/4rb0q1324646551.ps tmp/4rb0q1324646551.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ss821324646551.ps tmp/5ss821324646551.png",intern=TRUE))
character(0)
> try(system("convert tmp/63pmz1324646551.ps tmp/63pmz1324646551.png",intern=TRUE))
character(0)
> try(system("convert tmp/7q4wp1324646551.ps tmp/7q4wp1324646551.png",intern=TRUE))
character(0)
> try(system("convert tmp/8sm161324646551.ps tmp/8sm161324646551.png",intern=TRUE))
character(0)
> try(system("convert tmp/9t5tv1324646551.ps tmp/9t5tv1324646551.png",intern=TRUE))
character(0)
> try(system("convert tmp/10kopv1324646551.ps tmp/10kopv1324646551.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.276 0.799 5.102