R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(34,0,39,0,40,0,45,0,43,0,42,0,49,0,43,0,50,0,44,0,40,0,41,0,45,0,45,0,48,0,54,0,47,0,35,0,28,0,28,0,34,0,23,0,33,0,38,0,41,0,47,0,46,0,45,0,47,0,49,0,50,0,56,0,50,0,56,0,58,0,59,0,51,0,59,0,60,0,60,0,68,0,62,0,62,0,58,0,56,0,50,0,52,0,36,0,33,0,26,0,28,0,27,0,20,0,16,0,11,0,0,1,3,1,10,1,0,1,3,1),dim=c(2,60),dimnames=list(c('Eco','Val'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Eco','Val'),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
Eco Val M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 34 0 1 0 0 0 0 0 0 0 0 0 0 1
2 39 0 0 1 0 0 0 0 0 0 0 0 0 2
3 40 0 0 0 1 0 0 0 0 0 0 0 0 3
4 45 0 0 0 0 1 0 0 0 0 0 0 0 4
5 43 0 0 0 0 0 1 0 0 0 0 0 0 5
6 42 0 0 0 0 0 0 1 0 0 0 0 0 6
7 49 0 0 0 0 0 0 0 1 0 0 0 0 7
8 43 0 0 0 0 0 0 0 0 1 0 0 0 8
9 50 0 0 0 0 0 0 0 0 0 1 0 0 9
10 44 0 0 0 0 0 0 0 0 0 0 1 0 10
11 40 0 0 0 0 0 0 0 0 0 0 0 1 11
12 41 0 0 0 0 0 0 0 0 0 0 0 0 12
13 45 0 1 0 0 0 0 0 0 0 0 0 0 13
14 45 0 0 1 0 0 0 0 0 0 0 0 0 14
15 48 0 0 0 1 0 0 0 0 0 0 0 0 15
16 54 0 0 0 0 1 0 0 0 0 0 0 0 16
17 47 0 0 0 0 0 1 0 0 0 0 0 0 17
18 35 0 0 0 0 0 0 1 0 0 0 0 0 18
19 28 0 0 0 0 0 0 0 1 0 0 0 0 19
20 28 0 0 0 0 0 0 0 0 1 0 0 0 20
21 34 0 0 0 0 0 0 0 0 0 1 0 0 21
22 23 0 0 0 0 0 0 0 0 0 0 1 0 22
23 33 0 0 0 0 0 0 0 0 0 0 0 1 23
24 38 0 0 0 0 0 0 0 0 0 0 0 0 24
25 41 0 1 0 0 0 0 0 0 0 0 0 0 25
26 47 0 0 1 0 0 0 0 0 0 0 0 0 26
27 46 0 0 0 1 0 0 0 0 0 0 0 0 27
28 45 0 0 0 0 1 0 0 0 0 0 0 0 28
29 47 0 0 0 0 0 1 0 0 0 0 0 0 29
30 49 0 0 0 0 0 0 1 0 0 0 0 0 30
31 50 0 0 0 0 0 0 0 1 0 0 0 0 31
32 56 0 0 0 0 0 0 0 0 1 0 0 0 32
33 50 0 0 0 0 0 0 0 0 0 1 0 0 33
34 56 0 0 0 0 0 0 0 0 0 0 1 0 34
35 58 0 0 0 0 0 0 0 0 0 0 0 1 35
36 59 0 0 0 0 0 0 0 0 0 0 0 0 36
37 51 0 1 0 0 0 0 0 0 0 0 0 0 37
38 59 0 0 1 0 0 0 0 0 0 0 0 0 38
39 60 0 0 0 1 0 0 0 0 0 0 0 0 39
40 60 0 0 0 0 1 0 0 0 0 0 0 0 40
41 68 0 0 0 0 0 1 0 0 0 0 0 0 41
42 62 0 0 0 0 0 0 1 0 0 0 0 0 42
43 62 0 0 0 0 0 0 0 1 0 0 0 0 43
44 58 0 0 0 0 0 0 0 0 1 0 0 0 44
45 56 0 0 0 0 0 0 0 0 0 1 0 0 45
46 50 0 0 0 0 0 0 0 0 0 0 1 0 46
47 52 0 0 0 0 0 0 0 0 0 0 0 1 47
48 36 0 0 0 0 0 0 0 0 0 0 0 0 48
49 33 0 1 0 0 0 0 0 0 0 0 0 0 49
50 26 0 0 1 0 0 0 0 0 0 0 0 0 50
51 28 0 0 0 1 0 0 0 0 0 0 0 0 51
52 27 0 0 0 0 1 0 0 0 0 0 0 0 52
53 20 0 0 0 0 0 1 0 0 0 0 0 0 53
54 16 0 0 0 0 0 0 1 0 0 0 0 0 54
55 11 0 0 0 0 0 0 0 1 0 0 0 0 55
56 0 1 0 0 0 0 0 0 0 1 0 0 0 56
57 3 1 0 0 0 0 0 0 0 0 1 0 0 57
58 10 1 0 0 0 0 0 0 0 0 0 1 0 58
59 0 1 0 0 0 0 0 0 0 0 0 0 1 59
60 3 1 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) Val M1 M2 M3 M4
44.58632 -41.27368 -3.13939 -0.71351 0.51237 2.33825
M5 M6 M7 M8 M9 M10
1.16412 -3.01000 -3.78412 1.49649 3.12237 1.14825
M11 t
1.17412 -0.02588
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-28.3789 -6.0016 0.9095 8.2447 23.3105
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 44.58632 7.02323 6.348 8.68e-08 ***
Val -41.27368 7.51531 -5.492 1.66e-06 ***
M1 -3.13939 8.58234 -0.366 0.716
M2 -0.71351 8.57548 -0.083 0.934
M3 0.51237 8.57014 0.060 0.953
M4 2.33825 8.56633 0.273 0.786
M5 1.16412 8.56404 0.136 0.892
M6 -3.01000 8.56327 -0.352 0.727
M7 -3.78412 8.56404 -0.442 0.661
M8 1.49649 8.47056 0.177 0.861
M9 3.12237 8.46515 0.369 0.714
M10 1.14825 8.46129 0.136 0.893
M11 1.17412 8.45897 0.139 0.890
t -0.02588 0.11434 -0.226 0.822
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.37 on 46 degrees of freedom
Multiple R-squared: 0.4873, Adjusted R-squared: 0.3425
F-statistic: 3.364 on 13 and 46 DF, p-value: 0.001155
> 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.003376799 0.006753599 0.9966232
[2,] 0.016296943 0.032593886 0.9837031
[3,] 0.083108120 0.166216239 0.9168919
[4,] 0.082136680 0.164273360 0.9178633
[5,] 0.074172870 0.148345740 0.9258271
[6,] 0.117600092 0.235200184 0.8823999
[7,] 0.103008972 0.206017944 0.8969910
[8,] 0.083848471 0.167696942 0.9161515
[9,] 0.065247305 0.130494609 0.9347527
[10,] 0.051740459 0.103480917 0.9482595
[11,] 0.038779500 0.077559001 0.9612205
[12,] 0.030854644 0.061709287 0.9691454
[13,] 0.029288550 0.058577099 0.9707115
[14,] 0.037134896 0.074269792 0.9628651
[15,] 0.048621697 0.097243393 0.9513783
[16,] 0.089954826 0.179909651 0.9100452
[17,] 0.142221955 0.284443909 0.8577780
[18,] 0.328925467 0.657850933 0.6710745
[19,] 0.548849038 0.902301924 0.4511510
[20,] 0.712918390 0.574163221 0.2870816
[21,] 0.838588054 0.322823892 0.1614119
[22,] 0.800304047 0.399391906 0.1996960
[23,] 0.791802232 0.416395537 0.2081978
[24,] 0.823230785 0.353538430 0.1767692
[25,] 0.732037543 0.535924913 0.2679625
[26,] 0.610541892 0.778916216 0.3894581
[27,] 0.443063435 0.886126871 0.5569366
> postscript(file="/var/www/html/rcomp/tmp/13qf41228669245.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/2nolh1228669245.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/390wi1228669245.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/4n8yb1228669245.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/5y6ax1228669245.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.4210526 -4.8210526 -5.0210526 -1.8210526 -2.6210526 0.5789474
7 8 9 10 11 12
8.3789474 -2.8757895 2.5242105 -1.4757895 -5.4757895 -3.2757895
13 14 15 16 17 18
3.8894737 1.4894737 3.2894737 7.4894737 1.6894737 -6.1105263
19 20 21 22 23 24
-12.3105263 -17.5652632 -13.1652632 -22.1652632 -12.1652632 -5.9652632
25 26 27 28 29 30
0.2000000 3.8000000 1.6000000 -1.2000000 2.0000000 8.2000000
31 32 33 34 35 36
10.0000000 10.7452632 3.1452632 11.1452632 13.1452632 15.3452632
37 38 39 40 41 42
10.5105263 16.1105263 15.9105263 14.1105263 23.3105263 21.5105263
43 44 45 46 47 48
22.3105263 13.0557895 9.4557895 5.4557895 7.4557895 -7.3442105
49 50 51 52 53 54
-7.1789474 -16.5789474 -15.7789474 -18.5789474 -24.3789474 -24.1789474
55 56 57 58 59 60
-28.3789474 -3.3600000 -1.9600000 7.0400000 -2.9600000 1.2400000
> postscript(file="/var/www/html/rcomp/tmp/6mn261228669245.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.4210526 NA
1 -4.8210526 -7.4210526
2 -5.0210526 -4.8210526
3 -1.8210526 -5.0210526
4 -2.6210526 -1.8210526
5 0.5789474 -2.6210526
6 8.3789474 0.5789474
7 -2.8757895 8.3789474
8 2.5242105 -2.8757895
9 -1.4757895 2.5242105
10 -5.4757895 -1.4757895
11 -3.2757895 -5.4757895
12 3.8894737 -3.2757895
13 1.4894737 3.8894737
14 3.2894737 1.4894737
15 7.4894737 3.2894737
16 1.6894737 7.4894737
17 -6.1105263 1.6894737
18 -12.3105263 -6.1105263
19 -17.5652632 -12.3105263
20 -13.1652632 -17.5652632
21 -22.1652632 -13.1652632
22 -12.1652632 -22.1652632
23 -5.9652632 -12.1652632
24 0.2000000 -5.9652632
25 3.8000000 0.2000000
26 1.6000000 3.8000000
27 -1.2000000 1.6000000
28 2.0000000 -1.2000000
29 8.2000000 2.0000000
30 10.0000000 8.2000000
31 10.7452632 10.0000000
32 3.1452632 10.7452632
33 11.1452632 3.1452632
34 13.1452632 11.1452632
35 15.3452632 13.1452632
36 10.5105263 15.3452632
37 16.1105263 10.5105263
38 15.9105263 16.1105263
39 14.1105263 15.9105263
40 23.3105263 14.1105263
41 21.5105263 23.3105263
42 22.3105263 21.5105263
43 13.0557895 22.3105263
44 9.4557895 13.0557895
45 5.4557895 9.4557895
46 7.4557895 5.4557895
47 -7.3442105 7.4557895
48 -7.1789474 -7.3442105
49 -16.5789474 -7.1789474
50 -15.7789474 -16.5789474
51 -18.5789474 -15.7789474
52 -24.3789474 -18.5789474
53 -24.1789474 -24.3789474
54 -28.3789474 -24.1789474
55 -3.3600000 -28.3789474
56 -1.9600000 -3.3600000
57 7.0400000 -1.9600000
58 -2.9600000 7.0400000
59 1.2400000 -2.9600000
60 NA 1.2400000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -4.8210526 -7.4210526
[2,] -5.0210526 -4.8210526
[3,] -1.8210526 -5.0210526
[4,] -2.6210526 -1.8210526
[5,] 0.5789474 -2.6210526
[6,] 8.3789474 0.5789474
[7,] -2.8757895 8.3789474
[8,] 2.5242105 -2.8757895
[9,] -1.4757895 2.5242105
[10,] -5.4757895 -1.4757895
[11,] -3.2757895 -5.4757895
[12,] 3.8894737 -3.2757895
[13,] 1.4894737 3.8894737
[14,] 3.2894737 1.4894737
[15,] 7.4894737 3.2894737
[16,] 1.6894737 7.4894737
[17,] -6.1105263 1.6894737
[18,] -12.3105263 -6.1105263
[19,] -17.5652632 -12.3105263
[20,] -13.1652632 -17.5652632
[21,] -22.1652632 -13.1652632
[22,] -12.1652632 -22.1652632
[23,] -5.9652632 -12.1652632
[24,] 0.2000000 -5.9652632
[25,] 3.8000000 0.2000000
[26,] 1.6000000 3.8000000
[27,] -1.2000000 1.6000000
[28,] 2.0000000 -1.2000000
[29,] 8.2000000 2.0000000
[30,] 10.0000000 8.2000000
[31,] 10.7452632 10.0000000
[32,] 3.1452632 10.7452632
[33,] 11.1452632 3.1452632
[34,] 13.1452632 11.1452632
[35,] 15.3452632 13.1452632
[36,] 10.5105263 15.3452632
[37,] 16.1105263 10.5105263
[38,] 15.9105263 16.1105263
[39,] 14.1105263 15.9105263
[40,] 23.3105263 14.1105263
[41,] 21.5105263 23.3105263
[42,] 22.3105263 21.5105263
[43,] 13.0557895 22.3105263
[44,] 9.4557895 13.0557895
[45,] 5.4557895 9.4557895
[46,] 7.4557895 5.4557895
[47,] -7.3442105 7.4557895
[48,] -7.1789474 -7.3442105
[49,] -16.5789474 -7.1789474
[50,] -15.7789474 -16.5789474
[51,] -18.5789474 -15.7789474
[52,] -24.3789474 -18.5789474
[53,] -24.1789474 -24.3789474
[54,] -28.3789474 -24.1789474
[55,] -3.3600000 -28.3789474
[56,] -1.9600000 -3.3600000
[57,] 7.0400000 -1.9600000
[58,] -2.9600000 7.0400000
[59,] 1.2400000 -2.9600000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -4.8210526 -7.4210526
2 -5.0210526 -4.8210526
3 -1.8210526 -5.0210526
4 -2.6210526 -1.8210526
5 0.5789474 -2.6210526
6 8.3789474 0.5789474
7 -2.8757895 8.3789474
8 2.5242105 -2.8757895
9 -1.4757895 2.5242105
10 -5.4757895 -1.4757895
11 -3.2757895 -5.4757895
12 3.8894737 -3.2757895
13 1.4894737 3.8894737
14 3.2894737 1.4894737
15 7.4894737 3.2894737
16 1.6894737 7.4894737
17 -6.1105263 1.6894737
18 -12.3105263 -6.1105263
19 -17.5652632 -12.3105263
20 -13.1652632 -17.5652632
21 -22.1652632 -13.1652632
22 -12.1652632 -22.1652632
23 -5.9652632 -12.1652632
24 0.2000000 -5.9652632
25 3.8000000 0.2000000
26 1.6000000 3.8000000
27 -1.2000000 1.6000000
28 2.0000000 -1.2000000
29 8.2000000 2.0000000
30 10.0000000 8.2000000
31 10.7452632 10.0000000
32 3.1452632 10.7452632
33 11.1452632 3.1452632
34 13.1452632 11.1452632
35 15.3452632 13.1452632
36 10.5105263 15.3452632
37 16.1105263 10.5105263
38 15.9105263 16.1105263
39 14.1105263 15.9105263
40 23.3105263 14.1105263
41 21.5105263 23.3105263
42 22.3105263 21.5105263
43 13.0557895 22.3105263
44 9.4557895 13.0557895
45 5.4557895 9.4557895
46 7.4557895 5.4557895
47 -7.3442105 7.4557895
48 -7.1789474 -7.3442105
49 -16.5789474 -7.1789474
50 -15.7789474 -16.5789474
51 -18.5789474 -15.7789474
52 -24.3789474 -18.5789474
53 -24.1789474 -24.3789474
54 -28.3789474 -24.1789474
55 -3.3600000 -28.3789474
56 -1.9600000 -3.3600000
57 7.0400000 -1.9600000
58 -2.9600000 7.0400000
59 1.2400000 -2.9600000
> 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/7kvc21228669245.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/8zhlc1228669245.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/9xnox1228669245.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/10kiof1228669245.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/11ahw61228669245.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/12yas11228669245.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/13p3zv1228669245.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/14atj01228669245.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/15ez9g1228669246.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/16d9if1228669246.tab")
+ }
>
> system("convert tmp/13qf41228669245.ps tmp/13qf41228669245.png")
> system("convert tmp/2nolh1228669245.ps tmp/2nolh1228669245.png")
> system("convert tmp/390wi1228669245.ps tmp/390wi1228669245.png")
> system("convert tmp/4n8yb1228669245.ps tmp/4n8yb1228669245.png")
> system("convert tmp/5y6ax1228669245.ps tmp/5y6ax1228669245.png")
> system("convert tmp/6mn261228669245.ps tmp/6mn261228669245.png")
> system("convert tmp/7kvc21228669245.ps tmp/7kvc21228669245.png")
> system("convert tmp/8zhlc1228669245.ps tmp/8zhlc1228669245.png")
> system("convert tmp/9xnox1228669245.ps tmp/9xnox1228669245.png")
> system("convert tmp/10kiof1228669245.ps tmp/10kiof1228669245.png")
>
>
> proc.time()
user system elapsed
2.451 1.640 10.272