R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(0,345,0,334,0,345,0,333,0,336,0,324,0,320,0,330,0,313,0,301,0,288,0,294,0,302,0,294,0,293,0,290,0,283,0,286,0,293,0,334,0,329,0,411,0,416,0,418,0,408,0,402,0,401,0,400,0,389,0,371,0,364,0,350,0,332,0,323,0,316,0,312,0,315,0,314,0,313,0,314,0,317,0,308,0,312,0,306,0,304,0,297,0,284,0,278,1,273,1,265,1,259,1,252,1,245,1,235,1,232,1,229,1,219,1,218,1,215,1,211),dim=c(2,60),dimnames=list(c('rookverbod_of_niet','werklozen_tabakssector'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('rookverbod_of_niet','werklozen_tabakssector'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '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
werklozen_tabakssector rookverbod_of_niet M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 345 0 1 0 0 0 0 0 0 0 0 0 0
2 334 0 0 1 0 0 0 0 0 0 0 0 0
3 345 0 0 0 1 0 0 0 0 0 0 0 0
4 333 0 0 0 0 1 0 0 0 0 0 0 0
5 336 0 0 0 0 0 1 0 0 0 0 0 0
6 324 0 0 0 0 0 0 1 0 0 0 0 0
7 320 0 0 0 0 0 0 0 1 0 0 0 0
8 330 0 0 0 0 0 0 0 0 1 0 0 0
9 313 0 0 0 0 0 0 0 0 0 1 0 0
10 301 0 0 0 0 0 0 0 0 0 0 1 0
11 288 0 0 0 0 0 0 0 0 0 0 0 1
12 294 0 0 0 0 0 0 0 0 0 0 0 0
13 302 0 1 0 0 0 0 0 0 0 0 0 0
14 294 0 0 1 0 0 0 0 0 0 0 0 0
15 293 0 0 0 1 0 0 0 0 0 0 0 0
16 290 0 0 0 0 1 0 0 0 0 0 0 0
17 283 0 0 0 0 0 1 0 0 0 0 0 0
18 286 0 0 0 0 0 0 1 0 0 0 0 0
19 293 0 0 0 0 0 0 0 1 0 0 0 0
20 334 0 0 0 0 0 0 0 0 1 0 0 0
21 329 0 0 0 0 0 0 0 0 0 1 0 0
22 411 0 0 0 0 0 0 0 0 0 0 1 0
23 416 0 0 0 0 0 0 0 0 0 0 0 1
24 418 0 0 0 0 0 0 0 0 0 0 0 0
25 408 0 1 0 0 0 0 0 0 0 0 0 0
26 402 0 0 1 0 0 0 0 0 0 0 0 0
27 401 0 0 0 1 0 0 0 0 0 0 0 0
28 400 0 0 0 0 1 0 0 0 0 0 0 0
29 389 0 0 0 0 0 1 0 0 0 0 0 0
30 371 0 0 0 0 0 0 1 0 0 0 0 0
31 364 0 0 0 0 0 0 0 1 0 0 0 0
32 350 0 0 0 0 0 0 0 0 1 0 0 0
33 332 0 0 0 0 0 0 0 0 0 1 0 0
34 323 0 0 0 0 0 0 0 0 0 0 1 0
35 316 0 0 0 0 0 0 0 0 0 0 0 1
36 312 0 0 0 0 0 0 0 0 0 0 0 0
37 315 0 1 0 0 0 0 0 0 0 0 0 0
38 314 0 0 1 0 0 0 0 0 0 0 0 0
39 313 0 0 0 1 0 0 0 0 0 0 0 0
40 314 0 0 0 0 1 0 0 0 0 0 0 0
41 317 0 0 0 0 0 1 0 0 0 0 0 0
42 308 0 0 0 0 0 0 1 0 0 0 0 0
43 312 0 0 0 0 0 0 0 1 0 0 0 0
44 306 0 0 0 0 0 0 0 0 1 0 0 0
45 304 0 0 0 0 0 0 0 0 0 1 0 0
46 297 0 0 0 0 0 0 0 0 0 0 1 0
47 284 0 0 0 0 0 0 0 0 0 0 0 1
48 278 0 0 0 0 0 0 0 0 0 0 0 0
49 273 1 1 0 0 0 0 0 0 0 0 0 0
50 265 1 0 1 0 0 0 0 0 0 0 0 0
51 259 1 0 0 1 0 0 0 0 0 0 0 0
52 252 1 0 0 0 1 0 0 0 0 0 0 0
53 245 1 0 0 0 0 1 0 0 0 0 0 0
54 235 1 0 0 0 0 0 1 0 0 0 0 0
55 232 1 0 0 0 0 0 0 1 0 0 0 0
56 229 1 0 0 0 0 0 0 0 1 0 0 0
57 219 1 0 0 0 0 0 0 0 0 1 0 0
58 218 1 0 0 0 0 0 0 0 0 0 1 0
59 215 1 0 0 0 0 0 0 0 0 0 0 1
60 211 1 0 0 0 0 0 0 0 0 0 0 0
t
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26
27 27
28 28
29 29
30 30
31 31
32 32
33 33
34 34
35 35
36 36
37 37
38 38
39 39
40 40
41 41
42 42
43 43
44 44
45 45
46 46
47 47
48 48
49 49
50 50
51 51
52 52
53 53
54 54
55 55
56 56
57 57
58 58
59 59
60 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) rookverbod_of_niet M1 M2
324.8917 -88.4583 24.5944 17.9222
M3 M4 M5 M6
18.4500 14.1778 10.5056 1.4333
M7 M8 M9 M10
0.9611 6.6889 -3.5833 7.1444
M11 t
1.0722 -0.1278
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-50.225 -24.308 -6.063 10.440 96.175
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 324.8917 22.3749 14.520 < 2e-16 ***
rookverbod_of_niet -88.4583 18.3921 -4.810 1.66e-05 ***
M1 24.5944 25.9271 0.949 0.348
M2 17.9222 25.8509 0.693 0.492
M3 18.4500 25.7817 0.716 0.478
M4 14.1778 25.7197 0.551 0.584
M5 10.5056 25.6648 0.409 0.684
M6 1.4333 25.6172 0.056 0.956
M7 0.9611 25.5768 0.038 0.970
M8 6.6889 25.5437 0.262 0.795
M9 -3.5833 25.5180 -0.140 0.889
M10 7.1444 25.4996 0.280 0.781
M11 1.0722 25.4885 0.042 0.967
t -0.1278 0.4335 -0.295 0.770
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 40.3 on 46 degrees of freedom
Multiple R-squared: 0.5371, Adjusted R-squared: 0.4062
F-statistic: 4.105 on 13 and 46 DF, p-value: 0.0001797
> 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.004094172 8.188345e-03 9.959058e-01
[2,] 0.001715942 3.431885e-03 9.982841e-01
[3,] 0.006809566 1.361913e-02 9.931904e-01
[4,] 0.210032160 4.200643e-01 7.899678e-01
[5,] 0.853519319 2.929614e-01 1.464807e-01
[6,] 0.999894451 2.110986e-04 1.055493e-04
[7,] 0.999998289 3.421963e-06 1.710981e-06
[8,] 0.999999853 2.937155e-07 1.468577e-07
[9,] 0.999999878 2.431352e-07 1.215676e-07
[10,] 0.999999875 2.495961e-07 1.247981e-07
[11,] 0.999999898 2.046175e-07 1.023088e-07
[12,] 0.999999966 6.701759e-08 3.350879e-08
[13,] 0.999999983 3.412615e-08 1.706308e-08
[14,] 0.999999986 2.824010e-08 1.412005e-08
[15,] 0.999999983 3.454644e-08 1.727322e-08
[16,] 0.999999961 7.846784e-08 3.923392e-08
[17,] 0.999999820 3.590004e-07 1.795002e-07
[18,] 0.999999456 1.088232e-06 5.441162e-07
[19,] 0.999998155 3.690848e-06 1.845424e-06
[20,] 0.999993599 1.280263e-05 6.401315e-06
[21,] 0.999996705 6.589457e-06 3.294728e-06
[22,] 0.999997423 5.153165e-06 2.576583e-06
[23,] 0.999998034 3.932456e-06 1.966228e-06
[24,] 0.999995819 8.362067e-06 4.181033e-06
[25,] 0.999960610 7.878030e-05 3.939015e-05
[26,] 0.999632562 7.348754e-04 3.674377e-04
[27,] 0.996970103 6.059794e-03 3.029897e-03
> postscript(file="/var/www/html/rcomp/tmp/1ws9g1228998626.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/2ag5n1228998626.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/3ib031228998626.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/4thr11228998626.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/5js801228998626.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
-4.3583333 -8.5583333 2.0416667 -5.5583333 1.2416667 -1.5583333
7 8 9 10 11 12
-4.9583333 -0.5583333 -7.1583333 -29.7583333 -36.5583333 -29.3583333
13 14 15 16 17 18
-45.8250000 -47.0250000 -48.4250000 -47.0250000 -50.2250000 -38.0250000
19 20 21 22 23 24
-30.4250000 4.9750000 10.3750000 81.7750000 92.9750000 96.1750000
25 26 27 28 29 30
61.7083333 62.5083333 61.1083333 64.5083333 57.3083333 48.5083333
31 32 33 34 35 36
42.1083333 22.5083333 14.9083333 -4.6916667 -5.4916667 -8.2916667
37 38 39 40 41 42
-29.7583333 -23.9583333 -25.3583333 -19.9583333 -13.1583333 -12.9583333
43 44 45 46 47 48
-8.3583333 -19.9583333 -11.5583333 -29.1583333 -35.9583333 -40.7583333
49 50 51 52 53 54
18.2333333 17.0333333 10.6333333 8.0333333 4.8333333 4.0333333
55 56 57 58 59 60
1.6333333 -6.9666667 -6.5666667 -18.1666667 -14.9666667 -17.7666667
> postscript(file="/var/www/html/rcomp/tmp/6qjvz1228998626.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 -4.3583333 NA
1 -8.5583333 -4.3583333
2 2.0416667 -8.5583333
3 -5.5583333 2.0416667
4 1.2416667 -5.5583333
5 -1.5583333 1.2416667
6 -4.9583333 -1.5583333
7 -0.5583333 -4.9583333
8 -7.1583333 -0.5583333
9 -29.7583333 -7.1583333
10 -36.5583333 -29.7583333
11 -29.3583333 -36.5583333
12 -45.8250000 -29.3583333
13 -47.0250000 -45.8250000
14 -48.4250000 -47.0250000
15 -47.0250000 -48.4250000
16 -50.2250000 -47.0250000
17 -38.0250000 -50.2250000
18 -30.4250000 -38.0250000
19 4.9750000 -30.4250000
20 10.3750000 4.9750000
21 81.7750000 10.3750000
22 92.9750000 81.7750000
23 96.1750000 92.9750000
24 61.7083333 96.1750000
25 62.5083333 61.7083333
26 61.1083333 62.5083333
27 64.5083333 61.1083333
28 57.3083333 64.5083333
29 48.5083333 57.3083333
30 42.1083333 48.5083333
31 22.5083333 42.1083333
32 14.9083333 22.5083333
33 -4.6916667 14.9083333
34 -5.4916667 -4.6916667
35 -8.2916667 -5.4916667
36 -29.7583333 -8.2916667
37 -23.9583333 -29.7583333
38 -25.3583333 -23.9583333
39 -19.9583333 -25.3583333
40 -13.1583333 -19.9583333
41 -12.9583333 -13.1583333
42 -8.3583333 -12.9583333
43 -19.9583333 -8.3583333
44 -11.5583333 -19.9583333
45 -29.1583333 -11.5583333
46 -35.9583333 -29.1583333
47 -40.7583333 -35.9583333
48 18.2333333 -40.7583333
49 17.0333333 18.2333333
50 10.6333333 17.0333333
51 8.0333333 10.6333333
52 4.8333333 8.0333333
53 4.0333333 4.8333333
54 1.6333333 4.0333333
55 -6.9666667 1.6333333
56 -6.5666667 -6.9666667
57 -18.1666667 -6.5666667
58 -14.9666667 -18.1666667
59 -17.7666667 -14.9666667
60 NA -17.7666667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -8.5583333 -4.3583333
[2,] 2.0416667 -8.5583333
[3,] -5.5583333 2.0416667
[4,] 1.2416667 -5.5583333
[5,] -1.5583333 1.2416667
[6,] -4.9583333 -1.5583333
[7,] -0.5583333 -4.9583333
[8,] -7.1583333 -0.5583333
[9,] -29.7583333 -7.1583333
[10,] -36.5583333 -29.7583333
[11,] -29.3583333 -36.5583333
[12,] -45.8250000 -29.3583333
[13,] -47.0250000 -45.8250000
[14,] -48.4250000 -47.0250000
[15,] -47.0250000 -48.4250000
[16,] -50.2250000 -47.0250000
[17,] -38.0250000 -50.2250000
[18,] -30.4250000 -38.0250000
[19,] 4.9750000 -30.4250000
[20,] 10.3750000 4.9750000
[21,] 81.7750000 10.3750000
[22,] 92.9750000 81.7750000
[23,] 96.1750000 92.9750000
[24,] 61.7083333 96.1750000
[25,] 62.5083333 61.7083333
[26,] 61.1083333 62.5083333
[27,] 64.5083333 61.1083333
[28,] 57.3083333 64.5083333
[29,] 48.5083333 57.3083333
[30,] 42.1083333 48.5083333
[31,] 22.5083333 42.1083333
[32,] 14.9083333 22.5083333
[33,] -4.6916667 14.9083333
[34,] -5.4916667 -4.6916667
[35,] -8.2916667 -5.4916667
[36,] -29.7583333 -8.2916667
[37,] -23.9583333 -29.7583333
[38,] -25.3583333 -23.9583333
[39,] -19.9583333 -25.3583333
[40,] -13.1583333 -19.9583333
[41,] -12.9583333 -13.1583333
[42,] -8.3583333 -12.9583333
[43,] -19.9583333 -8.3583333
[44,] -11.5583333 -19.9583333
[45,] -29.1583333 -11.5583333
[46,] -35.9583333 -29.1583333
[47,] -40.7583333 -35.9583333
[48,] 18.2333333 -40.7583333
[49,] 17.0333333 18.2333333
[50,] 10.6333333 17.0333333
[51,] 8.0333333 10.6333333
[52,] 4.8333333 8.0333333
[53,] 4.0333333 4.8333333
[54,] 1.6333333 4.0333333
[55,] -6.9666667 1.6333333
[56,] -6.5666667 -6.9666667
[57,] -18.1666667 -6.5666667
[58,] -14.9666667 -18.1666667
[59,] -17.7666667 -14.9666667
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -8.5583333 -4.3583333
2 2.0416667 -8.5583333
3 -5.5583333 2.0416667
4 1.2416667 -5.5583333
5 -1.5583333 1.2416667
6 -4.9583333 -1.5583333
7 -0.5583333 -4.9583333
8 -7.1583333 -0.5583333
9 -29.7583333 -7.1583333
10 -36.5583333 -29.7583333
11 -29.3583333 -36.5583333
12 -45.8250000 -29.3583333
13 -47.0250000 -45.8250000
14 -48.4250000 -47.0250000
15 -47.0250000 -48.4250000
16 -50.2250000 -47.0250000
17 -38.0250000 -50.2250000
18 -30.4250000 -38.0250000
19 4.9750000 -30.4250000
20 10.3750000 4.9750000
21 81.7750000 10.3750000
22 92.9750000 81.7750000
23 96.1750000 92.9750000
24 61.7083333 96.1750000
25 62.5083333 61.7083333
26 61.1083333 62.5083333
27 64.5083333 61.1083333
28 57.3083333 64.5083333
29 48.5083333 57.3083333
30 42.1083333 48.5083333
31 22.5083333 42.1083333
32 14.9083333 22.5083333
33 -4.6916667 14.9083333
34 -5.4916667 -4.6916667
35 -8.2916667 -5.4916667
36 -29.7583333 -8.2916667
37 -23.9583333 -29.7583333
38 -25.3583333 -23.9583333
39 -19.9583333 -25.3583333
40 -13.1583333 -19.9583333
41 -12.9583333 -13.1583333
42 -8.3583333 -12.9583333
43 -19.9583333 -8.3583333
44 -11.5583333 -19.9583333
45 -29.1583333 -11.5583333
46 -35.9583333 -29.1583333
47 -40.7583333 -35.9583333
48 18.2333333 -40.7583333
49 17.0333333 18.2333333
50 10.6333333 17.0333333
51 8.0333333 10.6333333
52 4.8333333 8.0333333
53 4.0333333 4.8333333
54 1.6333333 4.0333333
55 -6.9666667 1.6333333
56 -6.5666667 -6.9666667
57 -18.1666667 -6.5666667
58 -14.9666667 -18.1666667
59 -17.7666667 -14.9666667
> 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/7jx331228998626.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/8mo9v1228998626.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/91ozi1228998626.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/10yc661228998626.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/11wwmn1228998626.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/12rmji1228998626.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/13fi2y1228998626.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/14x2jz1228998626.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/15vyny1228998627.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/16327n1228998627.tab")
+ }
>
> system("convert tmp/1ws9g1228998626.ps tmp/1ws9g1228998626.png")
> system("convert tmp/2ag5n1228998626.ps tmp/2ag5n1228998626.png")
> system("convert tmp/3ib031228998626.ps tmp/3ib031228998626.png")
> system("convert tmp/4thr11228998626.ps tmp/4thr11228998626.png")
> system("convert tmp/5js801228998626.ps tmp/5js801228998626.png")
> system("convert tmp/6qjvz1228998626.ps tmp/6qjvz1228998626.png")
> system("convert tmp/7jx331228998626.ps tmp/7jx331228998626.png")
> system("convert tmp/8mo9v1228998626.ps tmp/8mo9v1228998626.png")
> system("convert tmp/91ozi1228998626.ps tmp/91ozi1228998626.png")
> system("convert tmp/10yc661228998626.ps tmp/10yc661228998626.png")
>
>
> proc.time()
user system elapsed
2.404 1.596 7.813