R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(4.3,29,3.9,31,4,31,4.3,33,4.8,37,4.4,30,4.3,20,4.7,19,4.7,17,4.9,22,5,12,4.2,25,4.3,25,4.8,29,4.8,32,4.8,31,4.2,28,4.6,28,4.8,28,4.5,32,4.4,35,4.3,30,3.9,32,3.7,38,4,37,4.1,28,3.7,34,3.8,35,3.8,32,3.8,39,3.3,37,3.3,38,3.3,35,3.2,25,3.4,25,4.2,26,4.9,13,5.1,19,5.5,17,5.6,21,6.4,23,6.1,18,7.1,12,7.8,7,7.9,4,7.4,14,7.5,16,6.8,13,5.2,13,4.7,10,4.1,19,3.9,13,2.6,14,2.7,25,1.8,28,1,30,0.3,31,1.3,42,1,41,1.1,38),dim=c(2,60),dimnames=list(c('Consumentenprijsindex','Consumentenvertrouwen'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Consumentenprijsindex','Consumentenvertrouwen'),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
Consumentenprijsindex Consumentenvertrouwen M1 M2 M3 M4 M5 M6 M7 M8 M9 M10
1 4.3 29 1 0 0 0 0 0 0 0 0 0
2 3.9 31 0 1 0 0 0 0 0 0 0 0
3 4.0 31 0 0 1 0 0 0 0 0 0 0
4 4.3 33 0 0 0 1 0 0 0 0 0 0
5 4.8 37 0 0 0 0 1 0 0 0 0 0
6 4.4 30 0 0 0 0 0 1 0 0 0 0
7 4.3 20 0 0 0 0 0 0 1 0 0 0
8 4.7 19 0 0 0 0 0 0 0 1 0 0
9 4.7 17 0 0 0 0 0 0 0 0 1 0
10 4.9 22 0 0 0 0 0 0 0 0 0 1
11 5.0 12 0 0 0 0 0 0 0 0 0 0
12 4.2 25 0 0 0 0 0 0 0 0 0 0
13 4.3 25 1 0 0 0 0 0 0 0 0 0
14 4.8 29 0 1 0 0 0 0 0 0 0 0
15 4.8 32 0 0 1 0 0 0 0 0 0 0
16 4.8 31 0 0 0 1 0 0 0 0 0 0
17 4.2 28 0 0 0 0 1 0 0 0 0 0
18 4.6 28 0 0 0 0 0 1 0 0 0 0
19 4.8 28 0 0 0 0 0 0 1 0 0 0
20 4.5 32 0 0 0 0 0 0 0 1 0 0
21 4.4 35 0 0 0 0 0 0 0 0 1 0
22 4.3 30 0 0 0 0 0 0 0 0 0 1
23 3.9 32 0 0 0 0 0 0 0 0 0 0
24 3.7 38 0 0 0 0 0 0 0 0 0 0
25 4.0 37 1 0 0 0 0 0 0 0 0 0
26 4.1 28 0 1 0 0 0 0 0 0 0 0
27 3.7 34 0 0 1 0 0 0 0 0 0 0
28 3.8 35 0 0 0 1 0 0 0 0 0 0
29 3.8 32 0 0 0 0 1 0 0 0 0 0
30 3.8 39 0 0 0 0 0 1 0 0 0 0
31 3.3 37 0 0 0 0 0 0 1 0 0 0
32 3.3 38 0 0 0 0 0 0 0 1 0 0
33 3.3 35 0 0 0 0 0 0 0 0 1 0
34 3.2 25 0 0 0 0 0 0 0 0 0 1
35 3.4 25 0 0 0 0 0 0 0 0 0 0
36 4.2 26 0 0 0 0 0 0 0 0 0 0
37 4.9 13 1 0 0 0 0 0 0 0 0 0
38 5.1 19 0 1 0 0 0 0 0 0 0 0
39 5.5 17 0 0 1 0 0 0 0 0 0 0
40 5.6 21 0 0 0 1 0 0 0 0 0 0
41 6.4 23 0 0 0 0 1 0 0 0 0 0
42 6.1 18 0 0 0 0 0 1 0 0 0 0
43 7.1 12 0 0 0 0 0 0 1 0 0 0
44 7.8 7 0 0 0 0 0 0 0 1 0 0
45 7.9 4 0 0 0 0 0 0 0 0 1 0
46 7.4 14 0 0 0 0 0 0 0 0 0 1
47 7.5 16 0 0 0 0 0 0 0 0 0 0
48 6.8 13 0 0 0 0 0 0 0 0 0 0
49 5.2 13 1 0 0 0 0 0 0 0 0 0
50 4.7 10 0 1 0 0 0 0 0 0 0 0
51 4.1 19 0 0 1 0 0 0 0 0 0 0
52 3.9 13 0 0 0 1 0 0 0 0 0 0
53 2.6 14 0 0 0 0 1 0 0 0 0 0
54 2.7 25 0 0 0 0 0 1 0 0 0 0
55 1.8 28 0 0 0 0 0 0 1 0 0 0
56 1.0 30 0 0 0 0 0 0 0 1 0 0
57 0.3 31 0 0 0 0 0 0 0 0 1 0
58 1.3 42 0 0 0 0 0 0 0 0 0 1
59 1.0 41 0 0 0 0 0 0 0 0 0 0
60 1.1 38 0 0 0 0 0 0 0 0 0 0
M11 t
1 0 1
2 0 2
3 0 3
4 0 4
5 0 5
6 0 6
7 0 7
8 0 8
9 0 9
10 0 10
11 1 11
12 0 12
13 0 13
14 0 14
15 0 15
16 0 16
17 0 17
18 0 18
19 0 19
20 0 20
21 0 21
22 0 22
23 1 23
24 0 24
25 0 25
26 0 26
27 0 27
28 0 28
29 0 29
30 0 30
31 0 31
32 0 32
33 0 33
34 0 34
35 1 35
36 0 36
37 0 37
38 0 38
39 0 39
40 0 40
41 0 41
42 0 42
43 0 43
44 0 44
45 0 45
46 0 46
47 1 47
48 0 48
49 0 49
50 0 50
51 0 51
52 0 52
53 0 53
54 0 54
55 0 55
56 0 56
57 0 57
58 0 58
59 1 59
60 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Consumentenvertrouwen M1
8.77006 -0.12791 -0.41156
M2 M3 M4
-0.39855 -0.05622 0.03680
M5 M6 M7
-0.02460 0.12191 -0.28881
M8 M9 M10
-0.23021 -0.43953 -0.02511
M11 t
-0.23117 -0.03302
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.6049 -0.7110 0.1119 0.6752 2.5594
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.770060 0.863181 10.160 2.45e-13 ***
Consumentenvertrouwen -0.127911 0.017724 -7.217 4.31e-09 ***
M1 -0.411561 0.782238 -0.526 0.60132
M2 -0.398545 0.780763 -0.510 0.61217
M3 -0.056216 0.773942 -0.073 0.94241
M4 0.036800 0.772864 0.048 0.96223
M5 -0.024602 0.771719 -0.032 0.97471
M6 0.121906 0.770194 0.158 0.87493
M7 -0.288810 0.772205 -0.374 0.71012
M8 -0.230212 0.771228 -0.299 0.76667
M9 -0.439525 0.771835 -0.569 0.57182
M10 -0.025106 0.768855 -0.033 0.97409
M11 -0.231165 0.769878 -0.300 0.76533
t -0.033016 0.009545 -3.459 0.00118 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.214 on 46 degrees of freedom
Multiple R-squared: 0.5491, Adjusted R-squared: 0.4216
F-statistic: 4.309 on 13 and 46 DF, p-value: 0.0001097
> 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,] 6.011057e-03 1.202211e-02 0.9939889
[2,] 1.492220e-03 2.984439e-03 0.9985078
[3,] 3.254949e-03 6.509897e-03 0.9967451
[4,] 4.496109e-03 8.992218e-03 0.9955039
[5,] 1.968980e-03 3.937960e-03 0.9980310
[6,] 1.289284e-03 2.578569e-03 0.9987107
[7,] 1.050311e-03 2.100623e-03 0.9989497
[8,] 3.512176e-04 7.024351e-04 0.9996488
[9,] 1.236714e-04 2.473428e-04 0.9998763
[10,] 4.934731e-05 9.869462e-05 0.9999507
[11,] 2.960174e-05 5.920348e-05 0.9999704
[12,] 1.536694e-05 3.073388e-05 0.9999846
[13,] 6.846126e-06 1.369225e-05 0.9999932
[14,] 2.204739e-06 4.409479e-06 0.9999978
[15,] 1.303679e-06 2.607359e-06 0.9999987
[16,] 7.807426e-07 1.561485e-06 0.9999992
[17,] 4.980825e-07 9.961650e-07 0.9999995
[18,] 4.681118e-06 9.362237e-06 0.9999953
[19,] 7.734883e-05 1.546977e-04 0.9999227
[20,] 7.360185e-04 1.472037e-03 0.9992640
[21,] 1.071528e-02 2.143057e-02 0.9892847
[22,] 1.814945e-02 3.629890e-02 0.9818506
[23,] 1.363185e-01 2.726370e-01 0.8636815
[24,] 2.158446e-01 4.316893e-01 0.7841554
[25,] 5.486219e-01 9.027562e-01 0.4513781
[26,] 4.834788e-01 9.669575e-01 0.5165212
[27,] 4.422491e-01 8.844982e-01 0.5577509
> postscript(file="/var/www/html/rcomp/tmp/1buor1259414498.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/2map01259414498.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/3qw5n1259414498.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/4zycu1259414498.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/5a7dg1259414499.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
-0.31607680 -0.44025564 -0.64956950 -0.15374834 0.95231187 -0.45655489
7 8 9 10 11 12
-1.39192896 -1.14542165 -1.15891435 -0.70076472 -1.64079572 -0.97610779
13 14 15 16 17 18
-0.43153110 0.60011122 0.67452911 0.48661853 -0.40269533 -0.11618803
19 20 21 22 23 24
0.52754371 0.71360392 1.23966413 0.11870795 0.21360392 0.58291778
25 26 27 28 29 30
1.19958389 0.16838867 0.22653830 0.39444888 0.10513502 0.88701639
31 32 33 34 35 36
0.57492697 0.67725543 0.53585215 -1.22465693 -0.78558212 -0.05582116
37 38 39 40 41 42
-0.57408201 0.41338147 0.24824645 0.79988878 1.95012782 0.89708222
43 44 45 46 47 48
1.57335048 1.60821546 1.56681218 1.96451471 2.55941069 1.27752932
49 50 51 52 53 54
0.12210602 -0.74162573 -0.49974436 -1.52720784 -2.60487938 -1.21135569
55 56 57 58 59 60
-1.28389221 -1.85365316 -2.18341412 -0.15780101 -0.34663677 -0.82851814
> postscript(file="/var/www/html/rcomp/tmp/6zmo01259414499.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 -0.31607680 NA
1 -0.44025564 -0.31607680
2 -0.64956950 -0.44025564
3 -0.15374834 -0.64956950
4 0.95231187 -0.15374834
5 -0.45655489 0.95231187
6 -1.39192896 -0.45655489
7 -1.14542165 -1.39192896
8 -1.15891435 -1.14542165
9 -0.70076472 -1.15891435
10 -1.64079572 -0.70076472
11 -0.97610779 -1.64079572
12 -0.43153110 -0.97610779
13 0.60011122 -0.43153110
14 0.67452911 0.60011122
15 0.48661853 0.67452911
16 -0.40269533 0.48661853
17 -0.11618803 -0.40269533
18 0.52754371 -0.11618803
19 0.71360392 0.52754371
20 1.23966413 0.71360392
21 0.11870795 1.23966413
22 0.21360392 0.11870795
23 0.58291778 0.21360392
24 1.19958389 0.58291778
25 0.16838867 1.19958389
26 0.22653830 0.16838867
27 0.39444888 0.22653830
28 0.10513502 0.39444888
29 0.88701639 0.10513502
30 0.57492697 0.88701639
31 0.67725543 0.57492697
32 0.53585215 0.67725543
33 -1.22465693 0.53585215
34 -0.78558212 -1.22465693
35 -0.05582116 -0.78558212
36 -0.57408201 -0.05582116
37 0.41338147 -0.57408201
38 0.24824645 0.41338147
39 0.79988878 0.24824645
40 1.95012782 0.79988878
41 0.89708222 1.95012782
42 1.57335048 0.89708222
43 1.60821546 1.57335048
44 1.56681218 1.60821546
45 1.96451471 1.56681218
46 2.55941069 1.96451471
47 1.27752932 2.55941069
48 0.12210602 1.27752932
49 -0.74162573 0.12210602
50 -0.49974436 -0.74162573
51 -1.52720784 -0.49974436
52 -2.60487938 -1.52720784
53 -1.21135569 -2.60487938
54 -1.28389221 -1.21135569
55 -1.85365316 -1.28389221
56 -2.18341412 -1.85365316
57 -0.15780101 -2.18341412
58 -0.34663677 -0.15780101
59 -0.82851814 -0.34663677
60 NA -0.82851814
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.44025564 -0.31607680
[2,] -0.64956950 -0.44025564
[3,] -0.15374834 -0.64956950
[4,] 0.95231187 -0.15374834
[5,] -0.45655489 0.95231187
[6,] -1.39192896 -0.45655489
[7,] -1.14542165 -1.39192896
[8,] -1.15891435 -1.14542165
[9,] -0.70076472 -1.15891435
[10,] -1.64079572 -0.70076472
[11,] -0.97610779 -1.64079572
[12,] -0.43153110 -0.97610779
[13,] 0.60011122 -0.43153110
[14,] 0.67452911 0.60011122
[15,] 0.48661853 0.67452911
[16,] -0.40269533 0.48661853
[17,] -0.11618803 -0.40269533
[18,] 0.52754371 -0.11618803
[19,] 0.71360392 0.52754371
[20,] 1.23966413 0.71360392
[21,] 0.11870795 1.23966413
[22,] 0.21360392 0.11870795
[23,] 0.58291778 0.21360392
[24,] 1.19958389 0.58291778
[25,] 0.16838867 1.19958389
[26,] 0.22653830 0.16838867
[27,] 0.39444888 0.22653830
[28,] 0.10513502 0.39444888
[29,] 0.88701639 0.10513502
[30,] 0.57492697 0.88701639
[31,] 0.67725543 0.57492697
[32,] 0.53585215 0.67725543
[33,] -1.22465693 0.53585215
[34,] -0.78558212 -1.22465693
[35,] -0.05582116 -0.78558212
[36,] -0.57408201 -0.05582116
[37,] 0.41338147 -0.57408201
[38,] 0.24824645 0.41338147
[39,] 0.79988878 0.24824645
[40,] 1.95012782 0.79988878
[41,] 0.89708222 1.95012782
[42,] 1.57335048 0.89708222
[43,] 1.60821546 1.57335048
[44,] 1.56681218 1.60821546
[45,] 1.96451471 1.56681218
[46,] 2.55941069 1.96451471
[47,] 1.27752932 2.55941069
[48,] 0.12210602 1.27752932
[49,] -0.74162573 0.12210602
[50,] -0.49974436 -0.74162573
[51,] -1.52720784 -0.49974436
[52,] -2.60487938 -1.52720784
[53,] -1.21135569 -2.60487938
[54,] -1.28389221 -1.21135569
[55,] -1.85365316 -1.28389221
[56,] -2.18341412 -1.85365316
[57,] -0.15780101 -2.18341412
[58,] -0.34663677 -0.15780101
[59,] -0.82851814 -0.34663677
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.44025564 -0.31607680
2 -0.64956950 -0.44025564
3 -0.15374834 -0.64956950
4 0.95231187 -0.15374834
5 -0.45655489 0.95231187
6 -1.39192896 -0.45655489
7 -1.14542165 -1.39192896
8 -1.15891435 -1.14542165
9 -0.70076472 -1.15891435
10 -1.64079572 -0.70076472
11 -0.97610779 -1.64079572
12 -0.43153110 -0.97610779
13 0.60011122 -0.43153110
14 0.67452911 0.60011122
15 0.48661853 0.67452911
16 -0.40269533 0.48661853
17 -0.11618803 -0.40269533
18 0.52754371 -0.11618803
19 0.71360392 0.52754371
20 1.23966413 0.71360392
21 0.11870795 1.23966413
22 0.21360392 0.11870795
23 0.58291778 0.21360392
24 1.19958389 0.58291778
25 0.16838867 1.19958389
26 0.22653830 0.16838867
27 0.39444888 0.22653830
28 0.10513502 0.39444888
29 0.88701639 0.10513502
30 0.57492697 0.88701639
31 0.67725543 0.57492697
32 0.53585215 0.67725543
33 -1.22465693 0.53585215
34 -0.78558212 -1.22465693
35 -0.05582116 -0.78558212
36 -0.57408201 -0.05582116
37 0.41338147 -0.57408201
38 0.24824645 0.41338147
39 0.79988878 0.24824645
40 1.95012782 0.79988878
41 0.89708222 1.95012782
42 1.57335048 0.89708222
43 1.60821546 1.57335048
44 1.56681218 1.60821546
45 1.96451471 1.56681218
46 2.55941069 1.96451471
47 1.27752932 2.55941069
48 0.12210602 1.27752932
49 -0.74162573 0.12210602
50 -0.49974436 -0.74162573
51 -1.52720784 -0.49974436
52 -2.60487938 -1.52720784
53 -1.21135569 -2.60487938
54 -1.28389221 -1.21135569
55 -1.85365316 -1.28389221
56 -2.18341412 -1.85365316
57 -0.15780101 -2.18341412
58 -0.34663677 -0.15780101
59 -0.82851814 -0.34663677
> 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/71dyf1259414499.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/87y5b1259414499.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/9efd51259414499.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/10dhfd1259414499.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/11afk81259414499.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/128k531259414499.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/13xnah1259414499.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/14sfjy1259414499.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/15caaw1259414499.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/16db6c1259414499.tab")
+ }
>
> system("convert tmp/1buor1259414498.ps tmp/1buor1259414498.png")
> system("convert tmp/2map01259414498.ps tmp/2map01259414498.png")
> system("convert tmp/3qw5n1259414498.ps tmp/3qw5n1259414498.png")
> system("convert tmp/4zycu1259414498.ps tmp/4zycu1259414498.png")
> system("convert tmp/5a7dg1259414499.ps tmp/5a7dg1259414499.png")
> system("convert tmp/6zmo01259414499.ps tmp/6zmo01259414499.png")
> system("convert tmp/71dyf1259414499.ps tmp/71dyf1259414499.png")
> system("convert tmp/87y5b1259414499.ps tmp/87y5b1259414499.png")
> system("convert tmp/9efd51259414499.ps tmp/9efd51259414499.png")
> system("convert tmp/10dhfd1259414499.ps tmp/10dhfd1259414499.png")
>
>
> proc.time()
user system elapsed
2.408 1.556 3.126