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(878,0,859,0,845,0,806,0,786,0,757,0,769,0,771,0,726,0,680,0,661,0,658,0,642,1,613,1,612,1,596,1,592,1,586,1,619,1,618,1,593,1,574,1,564,1,556,1,564,1,539,1,512,1,507,1,491,1,460,1,478,1,495,1,461,1,418,1,409,1,428,1,420,1,403,1,385,1,377,1,376,1,377,1,374,1,374,1,363,1,354,1,342,1,344,1,362,1,360,1,349,1,339,1,346,1,345,1,363,1,356,1,337,1,337,1,344,1,351,1),dim=c(2,60),dimnames=list(c('werkloosheid_tabakssector','verbod_of_niet'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('werkloosheid_tabakssector','verbod_of_niet'),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 = '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
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
werkloosheid_tabakssector verbod_of_niet M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 878 0 1 0 0 0 0 0 0 0 0 0 0
2 859 0 0 1 0 0 0 0 0 0 0 0 0
3 845 0 0 0 1 0 0 0 0 0 0 0 0
4 806 0 0 0 0 1 0 0 0 0 0 0 0
5 786 0 0 0 0 0 1 0 0 0 0 0 0
6 757 0 0 0 0 0 0 1 0 0 0 0 0
7 769 0 0 0 0 0 0 0 1 0 0 0 0
8 771 0 0 0 0 0 0 0 0 1 0 0 0
9 726 0 0 0 0 0 0 0 0 0 1 0 0
10 680 0 0 0 0 0 0 0 0 0 0 1 0
11 661 0 0 0 0 0 0 0 0 0 0 0 1
12 658 0 0 0 0 0 0 0 0 0 0 0 0
13 642 1 1 0 0 0 0 0 0 0 0 0 0
14 613 1 0 1 0 0 0 0 0 0 0 0 0
15 612 1 0 0 1 0 0 0 0 0 0 0 0
16 596 1 0 0 0 1 0 0 0 0 0 0 0
17 592 1 0 0 0 0 1 0 0 0 0 0 0
18 586 1 0 0 0 0 0 1 0 0 0 0 0
19 619 1 0 0 0 0 0 0 1 0 0 0 0
20 618 1 0 0 0 0 0 0 0 1 0 0 0
21 593 1 0 0 0 0 0 0 0 0 1 0 0
22 574 1 0 0 0 0 0 0 0 0 0 1 0
23 564 1 0 0 0 0 0 0 0 0 0 0 1
24 556 1 0 0 0 0 0 0 0 0 0 0 0
25 564 1 1 0 0 0 0 0 0 0 0 0 0
26 539 1 0 1 0 0 0 0 0 0 0 0 0
27 512 1 0 0 1 0 0 0 0 0 0 0 0
28 507 1 0 0 0 1 0 0 0 0 0 0 0
29 491 1 0 0 0 0 1 0 0 0 0 0 0
30 460 1 0 0 0 0 0 1 0 0 0 0 0
31 478 1 0 0 0 0 0 0 1 0 0 0 0
32 495 1 0 0 0 0 0 0 0 1 0 0 0
33 461 1 0 0 0 0 0 0 0 0 1 0 0
34 418 1 0 0 0 0 0 0 0 0 0 1 0
35 409 1 0 0 0 0 0 0 0 0 0 0 1
36 428 1 0 0 0 0 0 0 0 0 0 0 0
37 420 1 1 0 0 0 0 0 0 0 0 0 0
38 403 1 0 1 0 0 0 0 0 0 0 0 0
39 385 1 0 0 1 0 0 0 0 0 0 0 0
40 377 1 0 0 0 1 0 0 0 0 0 0 0
41 376 1 0 0 0 0 1 0 0 0 0 0 0
42 377 1 0 0 0 0 0 1 0 0 0 0 0
43 374 1 0 0 0 0 0 0 1 0 0 0 0
44 374 1 0 0 0 0 0 0 0 1 0 0 0
45 363 1 0 0 0 0 0 0 0 0 1 0 0
46 354 1 0 0 0 0 0 0 0 0 0 1 0
47 342 1 0 0 0 0 0 0 0 0 0 0 1
48 344 1 0 0 0 0 0 0 0 0 0 0 0
49 362 1 1 0 0 0 0 0 0 0 0 0 0
50 360 1 0 1 0 0 0 0 0 0 0 0 0
51 349 1 0 0 1 0 0 0 0 0 0 0 0
52 339 1 0 0 0 1 0 0 0 0 0 0 0
53 346 1 0 0 0 0 1 0 0 0 0 0 0
54 345 1 0 0 0 0 0 1 0 0 0 0 0
55 363 1 0 0 0 0 0 0 1 0 0 0 0
56 356 1 0 0 0 0 0 0 0 1 0 0 0
57 337 1 0 0 0 0 0 0 0 0 1 0 0
58 337 1 0 0 0 0 0 0 0 0 0 1 0
59 344 1 0 0 0 0 0 0 0 0 0 0 1
60 351 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) verbod_of_niet M1 M2 M3
805.3917 -103.7083 27.8757 16.5597 9.4437
M4 M5 M6 M7 M8
0.9278 1.2118 -4.9042 17.7799 27.0639
M9 M10 M11 t
7.3479 -8.9681 -10.4840 -7.0840
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-62.383 -22.433 -2.512 25.254 74.358
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 805.3917 19.7510 40.777 < 2e-16 ***
verbod_of_niet -103.7083 16.9364 -6.123 1.89e-07 ***
M1 27.8757 23.8750 1.168 0.249
M2 16.5597 23.8048 0.696 0.490
M3 9.4437 23.7412 0.398 0.693
M4 0.9278 23.6840 0.039 0.969
M5 1.2118 23.6335 0.051 0.959
M6 -4.9042 23.5896 -0.208 0.836
M7 17.7799 23.5525 0.755 0.454
M8 27.0639 23.5220 1.151 0.256
M9 7.3479 23.4983 0.313 0.756
M10 -8.9681 23.4813 -0.382 0.704
M11 -10.4840 23.4711 -0.447 0.657
t -7.0840 0.3992 -17.746 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 37.11 on 46 degrees of freedom
Multiple R-squared: 0.9584, Adjusted R-squared: 0.9466
F-statistic: 81.47 on 13 and 46 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.1245067 0.24901350 0.875493251
[2,] 0.1829886 0.36597718 0.817011408
[3,] 0.2950243 0.59004854 0.704975732
[4,] 0.3234193 0.64683851 0.676580746
[5,] 0.4172788 0.83455754 0.582721230
[6,] 0.6029071 0.79418578 0.397092888
[7,] 0.7424569 0.51508611 0.257543054
[8,] 0.7918800 0.41623995 0.208119975
[9,] 0.8139435 0.37211300 0.186056502
[10,] 0.8263557 0.34728857 0.173644284
[11,] 0.8371661 0.32566770 0.162833851
[12,] 0.8810009 0.23799819 0.118999093
[13,] 0.8974210 0.20515795 0.102578973
[14,] 0.8658507 0.26829864 0.134149321
[15,] 0.8615978 0.27680439 0.138402197
[16,] 0.9478560 0.10428790 0.052143951
[17,] 0.9841733 0.03165330 0.015826651
[18,] 0.9784794 0.04304127 0.021520634
[19,] 0.9680244 0.06395129 0.031975646
[20,] 0.9865220 0.02695600 0.013477998
[21,] 0.9923323 0.01533532 0.007667661
[22,] 0.9904445 0.01911093 0.009555464
[23,] 0.9845213 0.03095747 0.015478737
[24,] 0.9805737 0.03885259 0.019426293
[25,] 0.9664037 0.06719255 0.033596276
[26,] 0.9604421 0.07911580 0.039557899
[27,] 0.8888667 0.22226663 0.111133314
> postscript(file="/var/www/html/rcomp/tmp/1w78i1229259447.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/2ieg61229259447.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/34qmn1229259447.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/4sfom1229259447.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/5osdu1229259447.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
51.8166667 51.2166667 51.4166667 28.0166667 14.8166667 -0.9833333
7 8 9 10 11 12
-4.5833333 -4.7833333 -22.9833333 -45.5833333 -55.9833333 -62.3833333
13 14 15 16 17 18
4.5333333 -6.0666667 7.1333333 6.7333333 9.5333333 16.7333333
19 20 21 22 23 24
34.1333333 30.9333333 32.7333333 37.1333333 35.7333333 24.3333333
25 26 27 28 29 30
11.5416667 4.9416667 -7.8583333 2.7416667 -6.4583333 -24.2583333
31 32 33 34 35 36
-21.8583333 -7.0583333 -14.2583333 -33.8583333 -34.2583333 -18.6583333
37 38 39 40 41 42
-47.4500000 -46.0500000 -49.8500000 -42.2500000 -36.4500000 -22.2500000
43 44 45 46 47 48
-40.8500000 -43.0500000 -27.2500000 -12.8500000 -16.2500000 -17.6500000
49 50 51 52 53 54
-20.4416667 -4.0416667 -0.8416667 4.7583333 18.5583333 30.7583333
55 56 57 58 59 60
33.1583333 23.9583333 31.7583333 55.1583333 70.7583333 74.3583333
> postscript(file="/var/www/html/rcomp/tmp/6v4im1229259447.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 51.8166667 NA
1 51.2166667 51.8166667
2 51.4166667 51.2166667
3 28.0166667 51.4166667
4 14.8166667 28.0166667
5 -0.9833333 14.8166667
6 -4.5833333 -0.9833333
7 -4.7833333 -4.5833333
8 -22.9833333 -4.7833333
9 -45.5833333 -22.9833333
10 -55.9833333 -45.5833333
11 -62.3833333 -55.9833333
12 4.5333333 -62.3833333
13 -6.0666667 4.5333333
14 7.1333333 -6.0666667
15 6.7333333 7.1333333
16 9.5333333 6.7333333
17 16.7333333 9.5333333
18 34.1333333 16.7333333
19 30.9333333 34.1333333
20 32.7333333 30.9333333
21 37.1333333 32.7333333
22 35.7333333 37.1333333
23 24.3333333 35.7333333
24 11.5416667 24.3333333
25 4.9416667 11.5416667
26 -7.8583333 4.9416667
27 2.7416667 -7.8583333
28 -6.4583333 2.7416667
29 -24.2583333 -6.4583333
30 -21.8583333 -24.2583333
31 -7.0583333 -21.8583333
32 -14.2583333 -7.0583333
33 -33.8583333 -14.2583333
34 -34.2583333 -33.8583333
35 -18.6583333 -34.2583333
36 -47.4500000 -18.6583333
37 -46.0500000 -47.4500000
38 -49.8500000 -46.0500000
39 -42.2500000 -49.8500000
40 -36.4500000 -42.2500000
41 -22.2500000 -36.4500000
42 -40.8500000 -22.2500000
43 -43.0500000 -40.8500000
44 -27.2500000 -43.0500000
45 -12.8500000 -27.2500000
46 -16.2500000 -12.8500000
47 -17.6500000 -16.2500000
48 -20.4416667 -17.6500000
49 -4.0416667 -20.4416667
50 -0.8416667 -4.0416667
51 4.7583333 -0.8416667
52 18.5583333 4.7583333
53 30.7583333 18.5583333
54 33.1583333 30.7583333
55 23.9583333 33.1583333
56 31.7583333 23.9583333
57 55.1583333 31.7583333
58 70.7583333 55.1583333
59 74.3583333 70.7583333
60 NA 74.3583333
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 51.2166667 51.8166667
[2,] 51.4166667 51.2166667
[3,] 28.0166667 51.4166667
[4,] 14.8166667 28.0166667
[5,] -0.9833333 14.8166667
[6,] -4.5833333 -0.9833333
[7,] -4.7833333 -4.5833333
[8,] -22.9833333 -4.7833333
[9,] -45.5833333 -22.9833333
[10,] -55.9833333 -45.5833333
[11,] -62.3833333 -55.9833333
[12,] 4.5333333 -62.3833333
[13,] -6.0666667 4.5333333
[14,] 7.1333333 -6.0666667
[15,] 6.7333333 7.1333333
[16,] 9.5333333 6.7333333
[17,] 16.7333333 9.5333333
[18,] 34.1333333 16.7333333
[19,] 30.9333333 34.1333333
[20,] 32.7333333 30.9333333
[21,] 37.1333333 32.7333333
[22,] 35.7333333 37.1333333
[23,] 24.3333333 35.7333333
[24,] 11.5416667 24.3333333
[25,] 4.9416667 11.5416667
[26,] -7.8583333 4.9416667
[27,] 2.7416667 -7.8583333
[28,] -6.4583333 2.7416667
[29,] -24.2583333 -6.4583333
[30,] -21.8583333 -24.2583333
[31,] -7.0583333 -21.8583333
[32,] -14.2583333 -7.0583333
[33,] -33.8583333 -14.2583333
[34,] -34.2583333 -33.8583333
[35,] -18.6583333 -34.2583333
[36,] -47.4500000 -18.6583333
[37,] -46.0500000 -47.4500000
[38,] -49.8500000 -46.0500000
[39,] -42.2500000 -49.8500000
[40,] -36.4500000 -42.2500000
[41,] -22.2500000 -36.4500000
[42,] -40.8500000 -22.2500000
[43,] -43.0500000 -40.8500000
[44,] -27.2500000 -43.0500000
[45,] -12.8500000 -27.2500000
[46,] -16.2500000 -12.8500000
[47,] -17.6500000 -16.2500000
[48,] -20.4416667 -17.6500000
[49,] -4.0416667 -20.4416667
[50,] -0.8416667 -4.0416667
[51,] 4.7583333 -0.8416667
[52,] 18.5583333 4.7583333
[53,] 30.7583333 18.5583333
[54,] 33.1583333 30.7583333
[55,] 23.9583333 33.1583333
[56,] 31.7583333 23.9583333
[57,] 55.1583333 31.7583333
[58,] 70.7583333 55.1583333
[59,] 74.3583333 70.7583333
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 51.2166667 51.8166667
2 51.4166667 51.2166667
3 28.0166667 51.4166667
4 14.8166667 28.0166667
5 -0.9833333 14.8166667
6 -4.5833333 -0.9833333
7 -4.7833333 -4.5833333
8 -22.9833333 -4.7833333
9 -45.5833333 -22.9833333
10 -55.9833333 -45.5833333
11 -62.3833333 -55.9833333
12 4.5333333 -62.3833333
13 -6.0666667 4.5333333
14 7.1333333 -6.0666667
15 6.7333333 7.1333333
16 9.5333333 6.7333333
17 16.7333333 9.5333333
18 34.1333333 16.7333333
19 30.9333333 34.1333333
20 32.7333333 30.9333333
21 37.1333333 32.7333333
22 35.7333333 37.1333333
23 24.3333333 35.7333333
24 11.5416667 24.3333333
25 4.9416667 11.5416667
26 -7.8583333 4.9416667
27 2.7416667 -7.8583333
28 -6.4583333 2.7416667
29 -24.2583333 -6.4583333
30 -21.8583333 -24.2583333
31 -7.0583333 -21.8583333
32 -14.2583333 -7.0583333
33 -33.8583333 -14.2583333
34 -34.2583333 -33.8583333
35 -18.6583333 -34.2583333
36 -47.4500000 -18.6583333
37 -46.0500000 -47.4500000
38 -49.8500000 -46.0500000
39 -42.2500000 -49.8500000
40 -36.4500000 -42.2500000
41 -22.2500000 -36.4500000
42 -40.8500000 -22.2500000
43 -43.0500000 -40.8500000
44 -27.2500000 -43.0500000
45 -12.8500000 -27.2500000
46 -16.2500000 -12.8500000
47 -17.6500000 -16.2500000
48 -20.4416667 -17.6500000
49 -4.0416667 -20.4416667
50 -0.8416667 -4.0416667
51 4.7583333 -0.8416667
52 18.5583333 4.7583333
53 30.7583333 18.5583333
54 33.1583333 30.7583333
55 23.9583333 33.1583333
56 31.7583333 23.9583333
57 55.1583333 31.7583333
58 70.7583333 55.1583333
59 74.3583333 70.7583333
> 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/7n8mh1229259447.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/869nu1229259447.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/9o8ss1229259447.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/10602j1229259447.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/1147zv1229259447.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/12a1ef1229259447.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/13k7v71229259447.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/14t69s1229259447.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/15y7v91229259447.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/16li811229259447.tab")
+ }
>
> system("convert tmp/1w78i1229259447.ps tmp/1w78i1229259447.png")
> system("convert tmp/2ieg61229259447.ps tmp/2ieg61229259447.png")
> system("convert tmp/34qmn1229259447.ps tmp/34qmn1229259447.png")
> system("convert tmp/4sfom1229259447.ps tmp/4sfom1229259447.png")
> system("convert tmp/5osdu1229259447.ps tmp/5osdu1229259447.png")
> system("convert tmp/6v4im1229259447.ps tmp/6v4im1229259447.png")
> system("convert tmp/7n8mh1229259447.ps tmp/7n8mh1229259447.png")
> system("convert tmp/869nu1229259447.ps tmp/869nu1229259447.png")
> system("convert tmp/9o8ss1229259447.ps tmp/9o8ss1229259447.png")
> system("convert tmp/10602j1229259447.ps tmp/10602j1229259447.png")
>
>
> proc.time()
user system elapsed
2.471 1.589 2.800