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(24,33,22,34,25,36,24,36,29,38,26,42,26,35,21,25,23,24,22,22,21,27,16,17,19,30,16,30,25,34,27,37,23,36,22,33,23,33,20,33,24,37,23,40,20,35,21,37,22,43,17,42,21,33,19,39,23,40,22,37,15,44,23,42,21,43,18,40,18,30,18,30,18,31,10,18,13,24,10,22,9,26,9,28,6,23,11,17,9,12,10,9,9,19,16,21,10,18,7,18,7,15,14,24,11,18,10,19,6,30,8,33,13,35,12,36,15,47,16,46,16,43),dim=c(2,61),dimnames=list(c('S.','E.S.'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('S.','E.S.'),1:61))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No 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
S. E.S. M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 24 33 1 0 0 0 0 0 0 0 0 0 0
2 22 34 0 1 0 0 0 0 0 0 0 0 0
3 25 36 0 0 1 0 0 0 0 0 0 0 0
4 24 36 0 0 0 1 0 0 0 0 0 0 0
5 29 38 0 0 0 0 1 0 0 0 0 0 0
6 26 42 0 0 0 0 0 1 0 0 0 0 0
7 26 35 0 0 0 0 0 0 1 0 0 0 0
8 21 25 0 0 0 0 0 0 0 1 0 0 0
9 23 24 0 0 0 0 0 0 0 0 1 0 0
10 22 22 0 0 0 0 0 0 0 0 0 1 0
11 21 27 0 0 0 0 0 0 0 0 0 0 1
12 16 17 0 0 0 0 0 0 0 0 0 0 0
13 19 30 1 0 0 0 0 0 0 0 0 0 0
14 16 30 0 1 0 0 0 0 0 0 0 0 0
15 25 34 0 0 1 0 0 0 0 0 0 0 0
16 27 37 0 0 0 1 0 0 0 0 0 0 0
17 23 36 0 0 0 0 1 0 0 0 0 0 0
18 22 33 0 0 0 0 0 1 0 0 0 0 0
19 23 33 0 0 0 0 0 0 1 0 0 0 0
20 20 33 0 0 0 0 0 0 0 1 0 0 0
21 24 37 0 0 0 0 0 0 0 0 1 0 0
22 23 40 0 0 0 0 0 0 0 0 0 1 0
23 20 35 0 0 0 0 0 0 0 0 0 0 1
24 21 37 0 0 0 0 0 0 0 0 0 0 0
25 22 43 1 0 0 0 0 0 0 0 0 0 0
26 17 42 0 1 0 0 0 0 0 0 0 0 0
27 21 33 0 0 1 0 0 0 0 0 0 0 0
28 19 39 0 0 0 1 0 0 0 0 0 0 0
29 23 40 0 0 0 0 1 0 0 0 0 0 0
30 22 37 0 0 0 0 0 1 0 0 0 0 0
31 15 44 0 0 0 0 0 0 1 0 0 0 0
32 23 42 0 0 0 0 0 0 0 1 0 0 0
33 21 43 0 0 0 0 0 0 0 0 1 0 0
34 18 40 0 0 0 0 0 0 0 0 0 1 0
35 18 30 0 0 0 0 0 0 0 0 0 0 1
36 18 30 0 0 0 0 0 0 0 0 0 0 0
37 18 31 1 0 0 0 0 0 0 0 0 0 0
38 10 18 0 1 0 0 0 0 0 0 0 0 0
39 13 24 0 0 1 0 0 0 0 0 0 0 0
40 10 22 0 0 0 1 0 0 0 0 0 0 0
41 9 26 0 0 0 0 1 0 0 0 0 0 0
42 9 28 0 0 0 0 0 1 0 0 0 0 0
43 6 23 0 0 0 0 0 0 1 0 0 0 0
44 11 17 0 0 0 0 0 0 0 1 0 0 0
45 9 12 0 0 0 0 0 0 0 0 1 0 0
46 10 9 0 0 0 0 0 0 0 0 0 1 0
47 9 19 0 0 0 0 0 0 0 0 0 0 1
48 16 21 0 0 0 0 0 0 0 0 0 0 0
49 10 18 1 0 0 0 0 0 0 0 0 0 0
50 7 18 0 1 0 0 0 0 0 0 0 0 0
51 7 15 0 0 1 0 0 0 0 0 0 0 0
52 14 24 0 0 0 1 0 0 0 0 0 0 0
53 11 18 0 0 0 0 1 0 0 0 0 0 0
54 10 19 0 0 0 0 0 1 0 0 0 0 0
55 6 30 0 0 0 0 0 0 1 0 0 0 0
56 8 33 0 0 0 0 0 0 0 1 0 0 0
57 13 35 0 0 0 0 0 0 0 0 1 0 0
58 12 36 0 0 0 0 0 0 0 0 0 1 0
59 15 47 0 0 0 0 0 0 0 0 0 0 1
60 16 46 0 0 0 0 0 0 0 0 0 0 0
61 16 43 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) E.S. M1 M2 M3 M4
5.2937 0.4009 -0.3558 -2.2784 1.5216 0.8388
M5 M6 M7 M8 M9 M10
1.0388 -0.2414 -3.3224 -0.7198 0.6000 -0.0793
M11
-1.3612
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-9.8026 -3.2310 0.6802 3.4362 9.9983
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.2937 3.3253 1.592 0.118
E.S. 0.4009 0.0760 5.275 3.15e-06 ***
M1 -0.3558 3.2649 -0.109 0.914
M2 -2.2784 3.4055 -0.669 0.507
M3 1.5216 3.4055 0.447 0.657
M4 0.8388 3.4045 0.246 0.806
M5 1.0388 3.4045 0.305 0.762
M6 -0.2414 3.4050 -0.071 0.944
M7 -3.3224 3.4095 -0.974 0.335
M8 -0.7198 3.4028 -0.212 0.833
M9 0.6000 3.4028 0.176 0.861
M10 -0.0793 3.4033 -0.023 0.982
M11 -1.3612 3.4045 -0.400 0.691
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.38 on 48 degrees of freedom
Multiple R-squared: 0.3965, Adjusted R-squared: 0.2456
F-statistic: 2.628 on 12 and 48 DF, p-value: 0.008815
> 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.03434451 0.06868902 0.96565549
[2,] 0.02601929 0.05203858 0.97398071
[3,] 0.06380903 0.12761806 0.93619097
[4,] 0.05456973 0.10913946 0.94543027
[5,] 0.10201046 0.20402092 0.89798954
[6,] 0.10160548 0.20321095 0.89839452
[7,] 0.08434114 0.16868228 0.91565886
[8,] 0.06293881 0.12587762 0.93706119
[9,] 0.03904101 0.07808202 0.96095899
[10,] 0.02445451 0.04890902 0.97554549
[11,] 0.02085731 0.04171463 0.97914269
[12,] 0.02643452 0.05286904 0.97356548
[13,] 0.05828302 0.11656605 0.94171698
[14,] 0.07452091 0.14904182 0.92547909
[15,] 0.10739732 0.21479463 0.89260268
[16,] 0.33894278 0.67788556 0.66105722
[17,] 0.55232670 0.89534660 0.44767330
[18,] 0.66499108 0.67001783 0.33500892
[19,] 0.72219180 0.55561641 0.27780820
[20,] 0.82801792 0.34396416 0.17198208
[21,] 0.78270118 0.43459764 0.21729882
[22,] 0.85231621 0.29536759 0.14768379
[23,] 0.88336129 0.23327742 0.11663871
[24,] 0.95195665 0.09608669 0.04804335
[25,] 0.97276491 0.05447017 0.02723509
[26,] 0.98034787 0.03930426 0.01965213
[27,] 0.97419073 0.05161855 0.02580927
[28,] 0.95157938 0.09684124 0.04842062
[29,] 0.96512197 0.06975606 0.03487803
[30,] 0.90389933 0.19220133 0.09610067
> postscript(file="/var/www/html/rcomp/tmp/1616f1260368420.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/22vgv1260368420.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/3ei9c1260368420.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/4i1c91260368420.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/5ltdq1260368420.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 = 61
Frequency = 1
1 2 3 4 5 6 7
5.8333333 5.3551281 3.7533881 3.4361721 7.4344321 4.1111262 9.9982600
8 9 10 11 12 13 14
6.4043499 7.4853939 7.9664379 6.2440019 3.8914838 2.0359433 0.9586080
15 16 17 18 19 20 21
4.5551281 6.0353021 2.2361721 3.7189560 7.8000000 2.1973901 3.2740841
22 23 24 25 26 27 28
1.7507782 2.0370421 0.8740841 -0.1753665 -2.8518318 0.9559981 -2.7664379
29 30 31 32 33 34 35
0.6326922 2.1154761 -4.6095698 1.5895602 -2.1311358 -3.2492218 2.0413920
36 37 38 39 40 41 42
0.6801740 0.6350733 -0.2309522 -3.4361721 -4.9516482 -7.7551281 -7.2766941
43 44 45 46 47 48 49
-5.1913002 -0.3886902 -1.7041663 1.1777476 -2.5490382 2.2880038 -2.1536169
50 51 52 53 54 55 56
-3.2309522 -5.8283422 -1.7533881 -2.5481682 -2.6688642 -7.9973901 -9.8026099
57 58 59 60 61
-6.9241759 -7.6457419 -7.7733977 -7.7337457 -6.1753665
> postscript(file="/var/www/html/rcomp/tmp/6olh31260368420.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 5.8333333 NA
1 5.3551281 5.8333333
2 3.7533881 5.3551281
3 3.4361721 3.7533881
4 7.4344321 3.4361721
5 4.1111262 7.4344321
6 9.9982600 4.1111262
7 6.4043499 9.9982600
8 7.4853939 6.4043499
9 7.9664379 7.4853939
10 6.2440019 7.9664379
11 3.8914838 6.2440019
12 2.0359433 3.8914838
13 0.9586080 2.0359433
14 4.5551281 0.9586080
15 6.0353021 4.5551281
16 2.2361721 6.0353021
17 3.7189560 2.2361721
18 7.8000000 3.7189560
19 2.1973901 7.8000000
20 3.2740841 2.1973901
21 1.7507782 3.2740841
22 2.0370421 1.7507782
23 0.8740841 2.0370421
24 -0.1753665 0.8740841
25 -2.8518318 -0.1753665
26 0.9559981 -2.8518318
27 -2.7664379 0.9559981
28 0.6326922 -2.7664379
29 2.1154761 0.6326922
30 -4.6095698 2.1154761
31 1.5895602 -4.6095698
32 -2.1311358 1.5895602
33 -3.2492218 -2.1311358
34 2.0413920 -3.2492218
35 0.6801740 2.0413920
36 0.6350733 0.6801740
37 -0.2309522 0.6350733
38 -3.4361721 -0.2309522
39 -4.9516482 -3.4361721
40 -7.7551281 -4.9516482
41 -7.2766941 -7.7551281
42 -5.1913002 -7.2766941
43 -0.3886902 -5.1913002
44 -1.7041663 -0.3886902
45 1.1777476 -1.7041663
46 -2.5490382 1.1777476
47 2.2880038 -2.5490382
48 -2.1536169 2.2880038
49 -3.2309522 -2.1536169
50 -5.8283422 -3.2309522
51 -1.7533881 -5.8283422
52 -2.5481682 -1.7533881
53 -2.6688642 -2.5481682
54 -7.9973901 -2.6688642
55 -9.8026099 -7.9973901
56 -6.9241759 -9.8026099
57 -7.6457419 -6.9241759
58 -7.7733977 -7.6457419
59 -7.7337457 -7.7733977
60 -6.1753665 -7.7337457
61 NA -6.1753665
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 5.3551281 5.8333333
[2,] 3.7533881 5.3551281
[3,] 3.4361721 3.7533881
[4,] 7.4344321 3.4361721
[5,] 4.1111262 7.4344321
[6,] 9.9982600 4.1111262
[7,] 6.4043499 9.9982600
[8,] 7.4853939 6.4043499
[9,] 7.9664379 7.4853939
[10,] 6.2440019 7.9664379
[11,] 3.8914838 6.2440019
[12,] 2.0359433 3.8914838
[13,] 0.9586080 2.0359433
[14,] 4.5551281 0.9586080
[15,] 6.0353021 4.5551281
[16,] 2.2361721 6.0353021
[17,] 3.7189560 2.2361721
[18,] 7.8000000 3.7189560
[19,] 2.1973901 7.8000000
[20,] 3.2740841 2.1973901
[21,] 1.7507782 3.2740841
[22,] 2.0370421 1.7507782
[23,] 0.8740841 2.0370421
[24,] -0.1753665 0.8740841
[25,] -2.8518318 -0.1753665
[26,] 0.9559981 -2.8518318
[27,] -2.7664379 0.9559981
[28,] 0.6326922 -2.7664379
[29,] 2.1154761 0.6326922
[30,] -4.6095698 2.1154761
[31,] 1.5895602 -4.6095698
[32,] -2.1311358 1.5895602
[33,] -3.2492218 -2.1311358
[34,] 2.0413920 -3.2492218
[35,] 0.6801740 2.0413920
[36,] 0.6350733 0.6801740
[37,] -0.2309522 0.6350733
[38,] -3.4361721 -0.2309522
[39,] -4.9516482 -3.4361721
[40,] -7.7551281 -4.9516482
[41,] -7.2766941 -7.7551281
[42,] -5.1913002 -7.2766941
[43,] -0.3886902 -5.1913002
[44,] -1.7041663 -0.3886902
[45,] 1.1777476 -1.7041663
[46,] -2.5490382 1.1777476
[47,] 2.2880038 -2.5490382
[48,] -2.1536169 2.2880038
[49,] -3.2309522 -2.1536169
[50,] -5.8283422 -3.2309522
[51,] -1.7533881 -5.8283422
[52,] -2.5481682 -1.7533881
[53,] -2.6688642 -2.5481682
[54,] -7.9973901 -2.6688642
[55,] -9.8026099 -7.9973901
[56,] -6.9241759 -9.8026099
[57,] -7.6457419 -6.9241759
[58,] -7.7733977 -7.6457419
[59,] -7.7337457 -7.7733977
[60,] -6.1753665 -7.7337457
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 5.3551281 5.8333333
2 3.7533881 5.3551281
3 3.4361721 3.7533881
4 7.4344321 3.4361721
5 4.1111262 7.4344321
6 9.9982600 4.1111262
7 6.4043499 9.9982600
8 7.4853939 6.4043499
9 7.9664379 7.4853939
10 6.2440019 7.9664379
11 3.8914838 6.2440019
12 2.0359433 3.8914838
13 0.9586080 2.0359433
14 4.5551281 0.9586080
15 6.0353021 4.5551281
16 2.2361721 6.0353021
17 3.7189560 2.2361721
18 7.8000000 3.7189560
19 2.1973901 7.8000000
20 3.2740841 2.1973901
21 1.7507782 3.2740841
22 2.0370421 1.7507782
23 0.8740841 2.0370421
24 -0.1753665 0.8740841
25 -2.8518318 -0.1753665
26 0.9559981 -2.8518318
27 -2.7664379 0.9559981
28 0.6326922 -2.7664379
29 2.1154761 0.6326922
30 -4.6095698 2.1154761
31 1.5895602 -4.6095698
32 -2.1311358 1.5895602
33 -3.2492218 -2.1311358
34 2.0413920 -3.2492218
35 0.6801740 2.0413920
36 0.6350733 0.6801740
37 -0.2309522 0.6350733
38 -3.4361721 -0.2309522
39 -4.9516482 -3.4361721
40 -7.7551281 -4.9516482
41 -7.2766941 -7.7551281
42 -5.1913002 -7.2766941
43 -0.3886902 -5.1913002
44 -1.7041663 -0.3886902
45 1.1777476 -1.7041663
46 -2.5490382 1.1777476
47 2.2880038 -2.5490382
48 -2.1536169 2.2880038
49 -3.2309522 -2.1536169
50 -5.8283422 -3.2309522
51 -1.7533881 -5.8283422
52 -2.5481682 -1.7533881
53 -2.6688642 -2.5481682
54 -7.9973901 -2.6688642
55 -9.8026099 -7.9973901
56 -6.9241759 -9.8026099
57 -7.6457419 -6.9241759
58 -7.7733977 -7.6457419
59 -7.7337457 -7.7733977
60 -6.1753665 -7.7337457
> 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/7z15h1260368420.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/8dmo71260368420.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/9hb2k1260368420.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/10l2pv1260368420.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/11flv41260368420.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/12yuap1260368420.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/13sqbf1260368420.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/14eqra1260368420.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/153va21260368420.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/16vdfj1260368420.tab")
+ }
>
> system("convert tmp/1616f1260368420.ps tmp/1616f1260368420.png")
> system("convert tmp/22vgv1260368420.ps tmp/22vgv1260368420.png")
> system("convert tmp/3ei9c1260368420.ps tmp/3ei9c1260368420.png")
> system("convert tmp/4i1c91260368420.ps tmp/4i1c91260368420.png")
> system("convert tmp/5ltdq1260368420.ps tmp/5ltdq1260368420.png")
> system("convert tmp/6olh31260368420.ps tmp/6olh31260368420.png")
> system("convert tmp/7z15h1260368420.ps tmp/7z15h1260368420.png")
> system("convert tmp/8dmo71260368420.ps tmp/8dmo71260368420.png")
> system("convert tmp/9hb2k1260368420.ps tmp/9hb2k1260368420.png")
> system("convert tmp/10l2pv1260368420.ps tmp/10l2pv1260368420.png")
>
>
> proc.time()
user system elapsed
2.434 1.547 4.157