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(902.2,0,891.9,0,874,0,930.9,0,944.2,0,935.9,0,937.1,0,885.1,0,892.4,0,987.3,0,946.3,0,799.6,0,875.4,0,846.2,0,880.6,0,885.7,0,868.9,0,882.5,0,789.6,0,773.3,0,804.3,0,817.8,0,836.7,0,721.8,0,760.8,0,841.4,0,1045.6,0,949.2,1,850.1,1,957.4,0,851.8,0,913.9,0,888,0,973.8,0,927.6,1,833,1,879.5,1,797.3,1,834.5,1,735.1,1,835,1,892.8,1,697.2,1,821.1,1,732.7,1,797.6,1,866.3,1,826.3,1,778.6,1,779.2,1,951,1,692.3,1,841.4,1,857.3,1,760.7,1,841.2,0,810.3,0,1007.4,1,931.3,0,931.2,0),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),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 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 902.2 0 1 0 0 0 0 0 0 0 0 0 0
2 891.9 0 0 1 0 0 0 0 0 0 0 0 0
3 874.0 0 0 0 1 0 0 0 0 0 0 0 0
4 930.9 0 0 0 0 1 0 0 0 0 0 0 0
5 944.2 0 0 0 0 0 1 0 0 0 0 0 0
6 935.9 0 0 0 0 0 0 1 0 0 0 0 0
7 937.1 0 0 0 0 0 0 0 1 0 0 0 0
8 885.1 0 0 0 0 0 0 0 0 1 0 0 0
9 892.4 0 0 0 0 0 0 0 0 0 1 0 0
10 987.3 0 0 0 0 0 0 0 0 0 0 1 0
11 946.3 0 0 0 0 0 0 0 0 0 0 0 1
12 799.6 0 0 0 0 0 0 0 0 0 0 0 0
13 875.4 0 1 0 0 0 0 0 0 0 0 0 0
14 846.2 0 0 1 0 0 0 0 0 0 0 0 0
15 880.6 0 0 0 1 0 0 0 0 0 0 0 0
16 885.7 0 0 0 0 1 0 0 0 0 0 0 0
17 868.9 0 0 0 0 0 1 0 0 0 0 0 0
18 882.5 0 0 0 0 0 0 1 0 0 0 0 0
19 789.6 0 0 0 0 0 0 0 1 0 0 0 0
20 773.3 0 0 0 0 0 0 0 0 1 0 0 0
21 804.3 0 0 0 0 0 0 0 0 0 1 0 0
22 817.8 0 0 0 0 0 0 0 0 0 0 1 0
23 836.7 0 0 0 0 0 0 0 0 0 0 0 1
24 721.8 0 0 0 0 0 0 0 0 0 0 0 0
25 760.8 0 1 0 0 0 0 0 0 0 0 0 0
26 841.4 0 0 1 0 0 0 0 0 0 0 0 0
27 1045.6 0 0 0 1 0 0 0 0 0 0 0 0
28 949.2 1 0 0 0 1 0 0 0 0 0 0 0
29 850.1 1 0 0 0 0 1 0 0 0 0 0 0
30 957.4 0 0 0 0 0 0 1 0 0 0 0 0
31 851.8 0 0 0 0 0 0 0 1 0 0 0 0
32 913.9 0 0 0 0 0 0 0 0 1 0 0 0
33 888.0 0 0 0 0 0 0 0 0 0 1 0 0
34 973.8 0 0 0 0 0 0 0 0 0 0 1 0
35 927.6 1 0 0 0 0 0 0 0 0 0 0 1
36 833.0 1 0 0 0 0 0 0 0 0 0 0 0
37 879.5 1 1 0 0 0 0 0 0 0 0 0 0
38 797.3 1 0 1 0 0 0 0 0 0 0 0 0
39 834.5 1 0 0 1 0 0 0 0 0 0 0 0
40 735.1 1 0 0 0 1 0 0 0 0 0 0 0
41 835.0 1 0 0 0 0 1 0 0 0 0 0 0
42 892.8 1 0 0 0 0 0 1 0 0 0 0 0
43 697.2 1 0 0 0 0 0 0 1 0 0 0 0
44 821.1 1 0 0 0 0 0 0 0 1 0 0 0
45 732.7 1 0 0 0 0 0 0 0 0 1 0 0
46 797.6 1 0 0 0 0 0 0 0 0 0 1 0
47 866.3 1 0 0 0 0 0 0 0 0 0 0 1
48 826.3 1 0 0 0 0 0 0 0 0 0 0 0
49 778.6 1 1 0 0 0 0 0 0 0 0 0 0
50 779.2 1 0 1 0 0 0 0 0 0 0 0 0
51 951.0 1 0 0 1 0 0 0 0 0 0 0 0
52 692.3 1 0 0 0 1 0 0 0 0 0 0 0
53 841.4 1 0 0 0 0 1 0 0 0 0 0 0
54 857.3 1 0 0 0 0 0 1 0 0 0 0 0
55 760.7 1 0 0 0 0 0 0 1 0 0 0 0
56 841.2 0 0 0 0 0 0 0 0 1 0 0 0
57 810.3 0 0 0 0 0 0 0 0 0 1 0 0
58 1007.4 1 0 0 0 0 0 0 0 0 0 1 0
59 931.3 0 0 0 0 0 0 0 0 0 0 0 1
60 931.2 0 0 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) X M1 M2 M3 M4
843.872 -53.731 16.920 8.820 94.760 27.006
M5 M6 M7 M8 M9 M10
56.286 82.800 -15.100 13.794 -7.586 94.400
M11
79.260
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-124.848 -40.422 1.005 41.770 132.052
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 843.872 31.331 26.934 <2e-16 ***
X -53.731 18.462 -2.910 0.0055 **
M1 16.920 43.061 0.393 0.6961
M2 8.820 43.061 0.205 0.8386
M3 94.760 43.061 2.201 0.0327 *
M4 27.006 43.219 0.625 0.5351
M5 56.286 43.219 1.302 0.1991
M6 82.800 43.061 1.923 0.0606 .
M7 -15.100 43.061 -0.351 0.7274
M8 13.794 43.219 0.319 0.7510
M9 -7.586 43.219 -0.176 0.8614
M10 94.400 43.061 2.192 0.0333 *
M11 79.260 43.061 1.841 0.0720 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 68.08 on 47 degrees of freedom
Multiple R-squared: 0.3687, Adjusted R-squared: 0.2075
F-statistic: 2.287 on 12 and 47 DF, p-value: 0.02180
> 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.0743658 0.1487316 0.92563421
[2,] 0.0883838 0.1767676 0.91161621
[3,] 0.0589700 0.1179400 0.94103000
[4,] 0.2030683 0.4061365 0.79693174
[5,] 0.2466312 0.4932624 0.75336878
[6,] 0.2192518 0.4385037 0.78074817
[7,] 0.4118674 0.8237348 0.58813258
[8,] 0.4396103 0.8792206 0.56038971
[9,] 0.5395616 0.9208768 0.46043839
[10,] 0.6936072 0.6127855 0.30639277
[11,] 0.6125059 0.7749883 0.38749413
[12,] 0.7347342 0.5305317 0.26526585
[13,] 0.9340698 0.1318604 0.06593021
[14,] 0.9094582 0.1810836 0.09054182
[15,] 0.8694707 0.2610587 0.13052933
[16,] 0.8225154 0.3549693 0.17748464
[17,] 0.7988891 0.4022219 0.20111094
[18,] 0.7900341 0.4199318 0.20996592
[19,] 0.7325432 0.5349136 0.26745678
[20,] 0.6878341 0.6243318 0.31216588
[21,] 0.6056208 0.7887583 0.39437915
[22,] 0.5978286 0.8043427 0.40217136
[23,] 0.5232885 0.9534229 0.47671145
[24,] 0.6031998 0.7936003 0.39680017
[25,] 0.6219124 0.7561752 0.37808762
[26,] 0.4974061 0.9948122 0.50259391
[27,] 0.3754040 0.7508081 0.62459597
[28,] 0.3239554 0.6479109 0.67604457
[29,] 0.2152308 0.4304616 0.78476922
> postscript(file="/var/www/html/rcomp/tmp/1am3t1258570107.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/2a3fh1258570107.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/34qaq1258570107.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/4d0x11258570107.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/58evw1258570107.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
41.407647 39.207647 -64.632353 60.021471 44.041471 9.227647
7 8 9 10 11 12
108.327647 27.433824 56.113824 49.027647 23.167647 -44.272353
13 14 15 16 17 18
14.607647 -6.492353 -58.032353 14.821471 -31.258529 -44.172353
19 20 21 22 23 24
-39.172353 -84.366176 -31.986176 -120.472353 -86.432353 -122.072353
25 26 27 28 29 30
-99.992353 -11.292353 106.967647 132.052353 3.672353 30.727647
31 32 33 34 35 36
23.027647 56.233824 51.713824 35.527647 58.198529 42.858529
37 38 39 40 41 42
72.438529 -1.661471 -50.401471 -82.047647 -11.427647 19.858529
43 44 45 46 47 48
-77.841471 17.164706 -49.855294 -86.941471 -3.101471 36.158529
49 50 51 52 53 54
-28.461471 -19.761471 66.098529 -124.847647 -5.027647 -15.641471
55 56 57 58 59 60
-14.341471 -16.466176 -25.986176 122.858529 8.167647 87.327647
> postscript(file="/var/www/html/rcomp/tmp/6aad01258570107.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 41.407647 NA
1 39.207647 41.407647
2 -64.632353 39.207647
3 60.021471 -64.632353
4 44.041471 60.021471
5 9.227647 44.041471
6 108.327647 9.227647
7 27.433824 108.327647
8 56.113824 27.433824
9 49.027647 56.113824
10 23.167647 49.027647
11 -44.272353 23.167647
12 14.607647 -44.272353
13 -6.492353 14.607647
14 -58.032353 -6.492353
15 14.821471 -58.032353
16 -31.258529 14.821471
17 -44.172353 -31.258529
18 -39.172353 -44.172353
19 -84.366176 -39.172353
20 -31.986176 -84.366176
21 -120.472353 -31.986176
22 -86.432353 -120.472353
23 -122.072353 -86.432353
24 -99.992353 -122.072353
25 -11.292353 -99.992353
26 106.967647 -11.292353
27 132.052353 106.967647
28 3.672353 132.052353
29 30.727647 3.672353
30 23.027647 30.727647
31 56.233824 23.027647
32 51.713824 56.233824
33 35.527647 51.713824
34 58.198529 35.527647
35 42.858529 58.198529
36 72.438529 42.858529
37 -1.661471 72.438529
38 -50.401471 -1.661471
39 -82.047647 -50.401471
40 -11.427647 -82.047647
41 19.858529 -11.427647
42 -77.841471 19.858529
43 17.164706 -77.841471
44 -49.855294 17.164706
45 -86.941471 -49.855294
46 -3.101471 -86.941471
47 36.158529 -3.101471
48 -28.461471 36.158529
49 -19.761471 -28.461471
50 66.098529 -19.761471
51 -124.847647 66.098529
52 -5.027647 -124.847647
53 -15.641471 -5.027647
54 -14.341471 -15.641471
55 -16.466176 -14.341471
56 -25.986176 -16.466176
57 122.858529 -25.986176
58 8.167647 122.858529
59 87.327647 8.167647
60 NA 87.327647
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 39.207647 41.407647
[2,] -64.632353 39.207647
[3,] 60.021471 -64.632353
[4,] 44.041471 60.021471
[5,] 9.227647 44.041471
[6,] 108.327647 9.227647
[7,] 27.433824 108.327647
[8,] 56.113824 27.433824
[9,] 49.027647 56.113824
[10,] 23.167647 49.027647
[11,] -44.272353 23.167647
[12,] 14.607647 -44.272353
[13,] -6.492353 14.607647
[14,] -58.032353 -6.492353
[15,] 14.821471 -58.032353
[16,] -31.258529 14.821471
[17,] -44.172353 -31.258529
[18,] -39.172353 -44.172353
[19,] -84.366176 -39.172353
[20,] -31.986176 -84.366176
[21,] -120.472353 -31.986176
[22,] -86.432353 -120.472353
[23,] -122.072353 -86.432353
[24,] -99.992353 -122.072353
[25,] -11.292353 -99.992353
[26,] 106.967647 -11.292353
[27,] 132.052353 106.967647
[28,] 3.672353 132.052353
[29,] 30.727647 3.672353
[30,] 23.027647 30.727647
[31,] 56.233824 23.027647
[32,] 51.713824 56.233824
[33,] 35.527647 51.713824
[34,] 58.198529 35.527647
[35,] 42.858529 58.198529
[36,] 72.438529 42.858529
[37,] -1.661471 72.438529
[38,] -50.401471 -1.661471
[39,] -82.047647 -50.401471
[40,] -11.427647 -82.047647
[41,] 19.858529 -11.427647
[42,] -77.841471 19.858529
[43,] 17.164706 -77.841471
[44,] -49.855294 17.164706
[45,] -86.941471 -49.855294
[46,] -3.101471 -86.941471
[47,] 36.158529 -3.101471
[48,] -28.461471 36.158529
[49,] -19.761471 -28.461471
[50,] 66.098529 -19.761471
[51,] -124.847647 66.098529
[52,] -5.027647 -124.847647
[53,] -15.641471 -5.027647
[54,] -14.341471 -15.641471
[55,] -16.466176 -14.341471
[56,] -25.986176 -16.466176
[57,] 122.858529 -25.986176
[58,] 8.167647 122.858529
[59,] 87.327647 8.167647
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 39.207647 41.407647
2 -64.632353 39.207647
3 60.021471 -64.632353
4 44.041471 60.021471
5 9.227647 44.041471
6 108.327647 9.227647
7 27.433824 108.327647
8 56.113824 27.433824
9 49.027647 56.113824
10 23.167647 49.027647
11 -44.272353 23.167647
12 14.607647 -44.272353
13 -6.492353 14.607647
14 -58.032353 -6.492353
15 14.821471 -58.032353
16 -31.258529 14.821471
17 -44.172353 -31.258529
18 -39.172353 -44.172353
19 -84.366176 -39.172353
20 -31.986176 -84.366176
21 -120.472353 -31.986176
22 -86.432353 -120.472353
23 -122.072353 -86.432353
24 -99.992353 -122.072353
25 -11.292353 -99.992353
26 106.967647 -11.292353
27 132.052353 106.967647
28 3.672353 132.052353
29 30.727647 3.672353
30 23.027647 30.727647
31 56.233824 23.027647
32 51.713824 56.233824
33 35.527647 51.713824
34 58.198529 35.527647
35 42.858529 58.198529
36 72.438529 42.858529
37 -1.661471 72.438529
38 -50.401471 -1.661471
39 -82.047647 -50.401471
40 -11.427647 -82.047647
41 19.858529 -11.427647
42 -77.841471 19.858529
43 17.164706 -77.841471
44 -49.855294 17.164706
45 -86.941471 -49.855294
46 -3.101471 -86.941471
47 36.158529 -3.101471
48 -28.461471 36.158529
49 -19.761471 -28.461471
50 66.098529 -19.761471
51 -124.847647 66.098529
52 -5.027647 -124.847647
53 -15.641471 -5.027647
54 -14.341471 -15.641471
55 -16.466176 -14.341471
56 -25.986176 -16.466176
57 122.858529 -25.986176
58 8.167647 122.858529
59 87.327647 8.167647
> 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/7qo1i1258570107.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/8xif21258570107.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/95k6w1258570107.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/10i2la1258570107.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/114iu61258570107.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/12lcto1258570107.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/13nn781258570107.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/14ph351258570107.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/15s27q1258570107.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/16t17k1258570107.tab")
+ }
>
> system("convert tmp/1am3t1258570107.ps tmp/1am3t1258570107.png")
> system("convert tmp/2a3fh1258570107.ps tmp/2a3fh1258570107.png")
> system("convert tmp/34qaq1258570107.ps tmp/34qaq1258570107.png")
> system("convert tmp/4d0x11258570107.ps tmp/4d0x11258570107.png")
> system("convert tmp/58evw1258570107.ps tmp/58evw1258570107.png")
> system("convert tmp/6aad01258570107.ps tmp/6aad01258570107.png")
> system("convert tmp/7qo1i1258570107.ps tmp/7qo1i1258570107.png")
> system("convert tmp/8xif21258570107.ps tmp/8xif21258570107.png")
> system("convert tmp/95k6w1258570107.ps tmp/95k6w1258570107.png")
> system("convert tmp/10i2la1258570107.ps tmp/10i2la1258570107.png")
>
>
> proc.time()
user system elapsed
2.429 1.580 2.862