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(12.6,18,15.7,16,13.2,19,20.3,18,12.8,23,8,20,0.9,20,3.6,15,14.1,17,21.7,16,24.5,15,18.9,10,13.9,13,11,10,5.8,19,15.5,21,22.4,17,31.7,16,30.3,17,31.4,14,20.2,18,19.7,17,10.8,14,13.2,15,15.1,16,15.6,11,15.5,15,12.7,13,10.9,17,10,16,9.1,9,10.3,17,16.9,15,22,12,27.6,12,28.9,12,31,12,32.9,4,38.1,7,28.8,4,29,3,21.8,3,28.8,0,25.6,5,28.2,3,20.2,4,17.9,3,16.3,10,13.2,4,8.1,1,4.5,1,-0.1,8,0,5,2.3,4,2.8,0,2.9,2,0.1,7,3.5,6,8.6,9,13.8,10),dim=c(2,60),dimnames=list(c('Rvnp','Svdg'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Rvnp','Svdg'),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
Rvnp Svdg M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 12.6 18 1 0 0 0 0 0 0 0 0 0 0 1
2 15.7 16 0 1 0 0 0 0 0 0 0 0 0 2
3 13.2 19 0 0 1 0 0 0 0 0 0 0 0 3
4 20.3 18 0 0 0 1 0 0 0 0 0 0 0 4
5 12.8 23 0 0 0 0 1 0 0 0 0 0 0 5
6 8.0 20 0 0 0 0 0 1 0 0 0 0 0 6
7 0.9 20 0 0 0 0 0 0 1 0 0 0 0 7
8 3.6 15 0 0 0 0 0 0 0 1 0 0 0 8
9 14.1 17 0 0 0 0 0 0 0 0 1 0 0 9
10 21.7 16 0 0 0 0 0 0 0 0 0 1 0 10
11 24.5 15 0 0 0 0 0 0 0 0 0 0 1 11
12 18.9 10 0 0 0 0 0 0 0 0 0 0 0 12
13 13.9 13 1 0 0 0 0 0 0 0 0 0 0 13
14 11.0 10 0 1 0 0 0 0 0 0 0 0 0 14
15 5.8 19 0 0 1 0 0 0 0 0 0 0 0 15
16 15.5 21 0 0 0 1 0 0 0 0 0 0 0 16
17 22.4 17 0 0 0 0 1 0 0 0 0 0 0 17
18 31.7 16 0 0 0 0 0 1 0 0 0 0 0 18
19 30.3 17 0 0 0 0 0 0 1 0 0 0 0 19
20 31.4 14 0 0 0 0 0 0 0 1 0 0 0 20
21 20.2 18 0 0 0 0 0 0 0 0 1 0 0 21
22 19.7 17 0 0 0 0 0 0 0 0 0 1 0 22
23 10.8 14 0 0 0 0 0 0 0 0 0 0 1 23
24 13.2 15 0 0 0 0 0 0 0 0 0 0 0 24
25 15.1 16 1 0 0 0 0 0 0 0 0 0 0 25
26 15.6 11 0 1 0 0 0 0 0 0 0 0 0 26
27 15.5 15 0 0 1 0 0 0 0 0 0 0 0 27
28 12.7 13 0 0 0 1 0 0 0 0 0 0 0 28
29 10.9 17 0 0 0 0 1 0 0 0 0 0 0 29
30 10.0 16 0 0 0 0 0 1 0 0 0 0 0 30
31 9.1 9 0 0 0 0 0 0 1 0 0 0 0 31
32 10.3 17 0 0 0 0 0 0 0 1 0 0 0 32
33 16.9 15 0 0 0 0 0 0 0 0 1 0 0 33
34 22.0 12 0 0 0 0 0 0 0 0 0 1 0 34
35 27.6 12 0 0 0 0 0 0 0 0 0 0 1 35
36 28.9 12 0 0 0 0 0 0 0 0 0 0 0 36
37 31.0 12 1 0 0 0 0 0 0 0 0 0 0 37
38 32.9 4 0 1 0 0 0 0 0 0 0 0 0 38
39 38.1 7 0 0 1 0 0 0 0 0 0 0 0 39
40 28.8 4 0 0 0 1 0 0 0 0 0 0 0 40
41 29.0 3 0 0 0 0 1 0 0 0 0 0 0 41
42 21.8 3 0 0 0 0 0 1 0 0 0 0 0 42
43 28.8 0 0 0 0 0 0 0 1 0 0 0 0 43
44 25.6 5 0 0 0 0 0 0 0 1 0 0 0 44
45 28.2 3 0 0 0 0 0 0 0 0 1 0 0 45
46 20.2 4 0 0 0 0 0 0 0 0 0 1 0 46
47 17.9 3 0 0 0 0 0 0 0 0 0 0 1 47
48 16.3 10 0 0 0 0 0 0 0 0 0 0 0 48
49 13.2 4 1 0 0 0 0 0 0 0 0 0 0 49
50 8.1 1 0 1 0 0 0 0 0 0 0 0 0 50
51 4.5 1 0 0 1 0 0 0 0 0 0 0 0 51
52 -0.1 8 0 0 0 1 0 0 0 0 0 0 0 52
53 0.0 5 0 0 0 0 1 0 0 0 0 0 0 53
54 2.3 4 0 0 0 0 0 1 0 0 0 0 0 54
55 2.8 0 0 0 0 0 0 0 1 0 0 0 0 55
56 2.9 2 0 0 0 0 0 0 0 1 0 0 0 56
57 0.1 7 0 0 0 0 0 0 0 0 1 0 0 57
58 3.5 6 0 0 0 0 0 0 0 0 0 1 0 58
59 8.6 9 0 0 0 0 0 0 0 0 0 0 1 59
60 13.8 10 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Svdg M1 M2 M3 M4
38.2335 -0.7649 -3.5927 -6.9918 -5.0113 -4.2186
M5 M6 M7 M8 M9 M10
-4.1720 -5.0362 -7.0914 -5.3267 -2.8021 -1.7334
M11 t
-1.2657 -0.3137
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-15.3229 -7.9292 -0.3123 5.1802 22.4666
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 38.2335 10.7038 3.572 0.000844 ***
Svdg -0.7649 0.4147 -1.845 0.071520 .
M1 -3.5927 6.6805 -0.538 0.593311
M2 -6.9918 7.0667 -0.989 0.327643
M3 -5.0113 6.6513 -0.753 0.455031
M4 -4.2186 6.6093 -0.638 0.526453
M5 -4.1720 6.5924 -0.633 0.529969
M6 -5.0362 6.6086 -0.762 0.449911
M7 -7.0914 6.7546 -1.050 0.299268
M8 -5.3267 6.6257 -0.804 0.425559
M9 -2.8021 6.5713 -0.426 0.671793
M10 -1.7334 6.5809 -0.263 0.793421
M11 -1.2657 6.5822 -0.192 0.848365
t -0.3137 0.1476 -2.125 0.038955 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.38 on 46 degrees of freedom
Multiple R-squared: 0.1056, Adjusted R-squared: -0.1472
F-statistic: 0.4177 on 13 and 46 DF, p-value: 0.9556
> 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.1117808 0.2235617 0.8882192
[2,] 0.4432000 0.8863999 0.5568000
[3,] 0.7107610 0.5784780 0.2892390
[4,] 0.7805086 0.4389828 0.2194914
[5,] 0.6765902 0.6468197 0.3234098
[6,] 0.5840496 0.8319007 0.4159504
[7,] 0.6940532 0.6118936 0.3059468
[8,] 0.7570014 0.4859972 0.2429986
[9,] 0.7174047 0.5651906 0.2825953
[10,] 0.6708686 0.6582628 0.3291314
[11,] 0.5997064 0.8005873 0.4002936
[12,] 0.6539661 0.6920677 0.3460339
[13,] 0.6110443 0.7779114 0.3889557
[14,] 0.5746772 0.8506455 0.4253228
[15,] 0.6518868 0.6962264 0.3481132
[16,] 0.6315920 0.7368161 0.3684080
[17,] 0.6356499 0.7287002 0.3643501
[18,] 0.6341347 0.7317307 0.3658653
[19,] 0.6707847 0.6584307 0.3292153
[20,] 0.8207001 0.3585998 0.1792999
[21,] 0.8178018 0.3643965 0.1821982
[22,] 0.7543238 0.4913524 0.2456762
[23,] 0.7466873 0.5066254 0.2533127
[24,] 0.7004823 0.5990354 0.2995177
[25,] 0.6703806 0.6592388 0.3296194
[26,] 0.5332699 0.9334602 0.4667301
[27,] 0.4347550 0.8695101 0.5652450
> postscript(file="/var/www/html/rcomp/tmp/1wjvv1258740670.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/2ol7m1258740670.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/3q8gg1258740670.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/4o96n1258740670.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/5962e1258740670.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
-7.9579922 -2.6751026 -4.5470624 1.3090179 -2.0992230 -8.0161327
7 8 9 10 11 12
-12.7472632 -15.3229420 -5.5039721 0.5760279 2.4570580 -7.9196509
13 14 15 16 17 18
-6.7183806 -8.2004409 -8.1827019 2.5682279 6.6754387 16.3884286
19 20 21 22 23 24
18.1222480 15.4764688 5.1253383 3.1053383 -8.2435312 -6.0305412
25 26 27 28 29 30
0.5408294 0.9288695 2.2218595 -2.5870100 -1.0602008 -1.5472108
31 32 33 34 35 36
-5.4329900 0.4356788 3.2948494 5.3449498 10.7909297 11.1389699
37 38 39 40 41 42
17.1453907 16.6385814 22.4666216 10.3928022 10.0948625 4.0728022
43 44 45 46 47 48
11.1468223 10.3206416 9.1798123 1.1897119 -2.0292580 0.7734308
49 50 51 52 53 54
-3.0098472 -6.6919075 -11.9587167 -11.6830380 -13.6108774 -10.8978874
55 56 57 58 59 60
-11.0888171 -10.9098472 -12.0960279 -10.2160279 -2.9751986 2.0377914
> postscript(file="/var/www/html/rcomp/tmp/698vb1258740670.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 -7.9579922 NA
1 -2.6751026 -7.9579922
2 -4.5470624 -2.6751026
3 1.3090179 -4.5470624
4 -2.0992230 1.3090179
5 -8.0161327 -2.0992230
6 -12.7472632 -8.0161327
7 -15.3229420 -12.7472632
8 -5.5039721 -15.3229420
9 0.5760279 -5.5039721
10 2.4570580 0.5760279
11 -7.9196509 2.4570580
12 -6.7183806 -7.9196509
13 -8.2004409 -6.7183806
14 -8.1827019 -8.2004409
15 2.5682279 -8.1827019
16 6.6754387 2.5682279
17 16.3884286 6.6754387
18 18.1222480 16.3884286
19 15.4764688 18.1222480
20 5.1253383 15.4764688
21 3.1053383 5.1253383
22 -8.2435312 3.1053383
23 -6.0305412 -8.2435312
24 0.5408294 -6.0305412
25 0.9288695 0.5408294
26 2.2218595 0.9288695
27 -2.5870100 2.2218595
28 -1.0602008 -2.5870100
29 -1.5472108 -1.0602008
30 -5.4329900 -1.5472108
31 0.4356788 -5.4329900
32 3.2948494 0.4356788
33 5.3449498 3.2948494
34 10.7909297 5.3449498
35 11.1389699 10.7909297
36 17.1453907 11.1389699
37 16.6385814 17.1453907
38 22.4666216 16.6385814
39 10.3928022 22.4666216
40 10.0948625 10.3928022
41 4.0728022 10.0948625
42 11.1468223 4.0728022
43 10.3206416 11.1468223
44 9.1798123 10.3206416
45 1.1897119 9.1798123
46 -2.0292580 1.1897119
47 0.7734308 -2.0292580
48 -3.0098472 0.7734308
49 -6.6919075 -3.0098472
50 -11.9587167 -6.6919075
51 -11.6830380 -11.9587167
52 -13.6108774 -11.6830380
53 -10.8978874 -13.6108774
54 -11.0888171 -10.8978874
55 -10.9098472 -11.0888171
56 -12.0960279 -10.9098472
57 -10.2160279 -12.0960279
58 -2.9751986 -10.2160279
59 2.0377914 -2.9751986
60 NA 2.0377914
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.6751026 -7.9579922
[2,] -4.5470624 -2.6751026
[3,] 1.3090179 -4.5470624
[4,] -2.0992230 1.3090179
[5,] -8.0161327 -2.0992230
[6,] -12.7472632 -8.0161327
[7,] -15.3229420 -12.7472632
[8,] -5.5039721 -15.3229420
[9,] 0.5760279 -5.5039721
[10,] 2.4570580 0.5760279
[11,] -7.9196509 2.4570580
[12,] -6.7183806 -7.9196509
[13,] -8.2004409 -6.7183806
[14,] -8.1827019 -8.2004409
[15,] 2.5682279 -8.1827019
[16,] 6.6754387 2.5682279
[17,] 16.3884286 6.6754387
[18,] 18.1222480 16.3884286
[19,] 15.4764688 18.1222480
[20,] 5.1253383 15.4764688
[21,] 3.1053383 5.1253383
[22,] -8.2435312 3.1053383
[23,] -6.0305412 -8.2435312
[24,] 0.5408294 -6.0305412
[25,] 0.9288695 0.5408294
[26,] 2.2218595 0.9288695
[27,] -2.5870100 2.2218595
[28,] -1.0602008 -2.5870100
[29,] -1.5472108 -1.0602008
[30,] -5.4329900 -1.5472108
[31,] 0.4356788 -5.4329900
[32,] 3.2948494 0.4356788
[33,] 5.3449498 3.2948494
[34,] 10.7909297 5.3449498
[35,] 11.1389699 10.7909297
[36,] 17.1453907 11.1389699
[37,] 16.6385814 17.1453907
[38,] 22.4666216 16.6385814
[39,] 10.3928022 22.4666216
[40,] 10.0948625 10.3928022
[41,] 4.0728022 10.0948625
[42,] 11.1468223 4.0728022
[43,] 10.3206416 11.1468223
[44,] 9.1798123 10.3206416
[45,] 1.1897119 9.1798123
[46,] -2.0292580 1.1897119
[47,] 0.7734308 -2.0292580
[48,] -3.0098472 0.7734308
[49,] -6.6919075 -3.0098472
[50,] -11.9587167 -6.6919075
[51,] -11.6830380 -11.9587167
[52,] -13.6108774 -11.6830380
[53,] -10.8978874 -13.6108774
[54,] -11.0888171 -10.8978874
[55,] -10.9098472 -11.0888171
[56,] -12.0960279 -10.9098472
[57,] -10.2160279 -12.0960279
[58,] -2.9751986 -10.2160279
[59,] 2.0377914 -2.9751986
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.6751026 -7.9579922
2 -4.5470624 -2.6751026
3 1.3090179 -4.5470624
4 -2.0992230 1.3090179
5 -8.0161327 -2.0992230
6 -12.7472632 -8.0161327
7 -15.3229420 -12.7472632
8 -5.5039721 -15.3229420
9 0.5760279 -5.5039721
10 2.4570580 0.5760279
11 -7.9196509 2.4570580
12 -6.7183806 -7.9196509
13 -8.2004409 -6.7183806
14 -8.1827019 -8.2004409
15 2.5682279 -8.1827019
16 6.6754387 2.5682279
17 16.3884286 6.6754387
18 18.1222480 16.3884286
19 15.4764688 18.1222480
20 5.1253383 15.4764688
21 3.1053383 5.1253383
22 -8.2435312 3.1053383
23 -6.0305412 -8.2435312
24 0.5408294 -6.0305412
25 0.9288695 0.5408294
26 2.2218595 0.9288695
27 -2.5870100 2.2218595
28 -1.0602008 -2.5870100
29 -1.5472108 -1.0602008
30 -5.4329900 -1.5472108
31 0.4356788 -5.4329900
32 3.2948494 0.4356788
33 5.3449498 3.2948494
34 10.7909297 5.3449498
35 11.1389699 10.7909297
36 17.1453907 11.1389699
37 16.6385814 17.1453907
38 22.4666216 16.6385814
39 10.3928022 22.4666216
40 10.0948625 10.3928022
41 4.0728022 10.0948625
42 11.1468223 4.0728022
43 10.3206416 11.1468223
44 9.1798123 10.3206416
45 1.1897119 9.1798123
46 -2.0292580 1.1897119
47 0.7734308 -2.0292580
48 -3.0098472 0.7734308
49 -6.6919075 -3.0098472
50 -11.9587167 -6.6919075
51 -11.6830380 -11.9587167
52 -13.6108774 -11.6830380
53 -10.8978874 -13.6108774
54 -11.0888171 -10.8978874
55 -10.9098472 -11.0888171
56 -12.0960279 -10.9098472
57 -10.2160279 -12.0960279
58 -2.9751986 -10.2160279
59 2.0377914 -2.9751986
> 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/7w4ir1258740670.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/8xw531258740670.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/92fik1258740670.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/100zh41258740670.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/11z1mp1258740670.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/12f9p21258740670.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/132r961258740670.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/14v0c81258740670.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/15m1cy1258740670.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/165lt41258740670.tab")
+ }
> system("convert tmp/1wjvv1258740670.ps tmp/1wjvv1258740670.png")
> system("convert tmp/2ol7m1258740670.ps tmp/2ol7m1258740670.png")
> system("convert tmp/3q8gg1258740670.ps tmp/3q8gg1258740670.png")
> system("convert tmp/4o96n1258740670.ps tmp/4o96n1258740670.png")
> system("convert tmp/5962e1258740670.ps tmp/5962e1258740670.png")
> system("convert tmp/698vb1258740670.ps tmp/698vb1258740670.png")
> system("convert tmp/7w4ir1258740670.ps tmp/7w4ir1258740670.png")
> system("convert tmp/8xw531258740670.ps tmp/8xw531258740670.png")
> system("convert tmp/92fik1258740670.ps tmp/92fik1258740670.png")
> system("convert tmp/100zh41258740670.ps tmp/100zh41258740670.png")
>
>
> proc.time()
user system elapsed
2.472 1.584 2.857