R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(1217.00
+ ,1210.00
+ ,31.00
+ ,48.00
+ ,961.00
+ ,2304.00
+ ,1488.00
+ ,19.00
+ ,30.00
+ ,1202.00
+ ,1209.00
+ ,34.40
+ ,38.00
+ ,1183.36
+ ,1444.00
+ ,1307.20
+ ,18.30
+ ,29.95
+ ,1180.00
+ ,1207.00
+ ,35.60
+ ,37.00
+ ,1267.36
+ ,1369.00
+ ,1317.20
+ ,18.90
+ ,29.94
+ ,1167.00
+ ,1206.00
+ ,32.80
+ ,48.00
+ ,1075.84
+ ,2304.00
+ ,1574.40
+ ,20.60
+ ,29.83
+ ,1186.00
+ ,1204.00
+ ,23.30
+ ,81.00
+ ,542.89
+ ,6561.00
+ ,1887.30
+ ,20.00
+ ,29.85
+ ,1168.00
+ ,1201.00
+ ,20.00
+ ,58.00
+ ,400.00
+ ,3364.00
+ ,1160.00
+ ,11.76
+ ,29.92
+ ,1142.00
+ ,1199.00
+ ,16.70
+ ,93.00
+ ,278.89
+ ,8649.00
+ ,1553.10
+ ,15.60
+ ,29.95
+ ,1147.00
+ ,1198.00
+ ,17.80
+ ,86.00
+ ,316.84
+ ,7396.00
+ ,1530.80
+ ,15.60
+ ,29.94
+ ,1183.00
+ ,1196.00
+ ,21.20
+ ,68.00
+ ,449.44
+ ,4624.00
+ ,1441.60
+ ,15.80
+ ,29.94
+ ,1149.00
+ ,1195.00
+ ,23.90
+ ,68.00
+ ,571.21
+ ,4624.00
+ ,1625.20
+ ,17.80
+ ,30.00
+ ,1197.00
+ ,1193.00
+ ,28.80
+ ,68.00
+ ,829.44
+ ,4624.00
+ ,1958.40
+ ,16.70
+ ,30.03
+ ,1210.00
+ ,1191.00
+ ,25.60
+ ,59.00
+ ,655.36
+ ,3481.00
+ ,1510.40
+ ,17.20
+ ,29.99
+ ,1206.00
+ ,1190.00
+ ,29.40
+ ,43.00
+ ,864.36
+ ,1849.00
+ ,1264.20
+ ,15.60
+ ,29.89
+ ,1196.00
+ ,1188.00
+ ,22.80
+ ,59.00
+ ,519.84
+ ,3481.00
+ ,1345.20
+ ,14.40
+ ,29.98
+ ,1190.00
+ ,1187.00
+ ,16.10
+ ,31.00
+ ,259.21
+ ,961.00
+ ,499.10
+ ,-0.60
+ ,30.26
+ ,1175.00
+ ,1185.00
+ ,16.10
+ ,49.00
+ ,259.21
+ ,2401.00
+ ,788.90
+ ,5.60
+ ,30.26
+ ,1186.00
+ ,1183.00
+ ,20.00
+ ,52.00
+ ,400.00
+ ,2704.00
+ ,1040.00
+ ,10.08
+ ,30.23
+ ,1172.00
+ ,1182.00
+ ,20.60
+ ,75.00
+ ,424.36
+ ,5625.00
+ ,1545.00
+ ,16.10
+ ,30.16
+ ,1152.00
+ ,1185.00
+ ,18.30
+ ,90.00
+ ,334.89
+ ,8100.00
+ ,1647.00
+ ,16.70
+ ,30.00
+ ,1154.00
+ ,1179.00
+ ,21.60
+ ,86.00
+ ,466.56
+ ,7396.00
+ ,1857.60
+ ,18.30
+ ,30.60
+ ,1168.00
+ ,1177.00
+ ,22.80
+ ,87.00
+ ,519.84
+ ,7569.00
+ ,1983.60
+ ,20.60
+ ,30.00
+ ,1180.00
+ ,1175.00
+ ,22.80
+ ,47.00
+ ,519.84
+ ,2209.00
+ ,1071.60
+ ,11.10
+ ,30.06
+ ,1169.00
+ ,1174.00
+ ,17.20
+ ,70.00
+ ,295.84
+ ,4900.00
+ ,1204.00
+ ,11.70
+ ,30.01
+ ,1166.00
+ ,1170.00
+ ,22.20
+ ,61.00
+ ,492.84
+ ,3721.00
+ ,1354.20
+ ,14.40
+ ,29.86
+ ,1177.00
+ ,1169.00
+ ,20.60
+ ,48.00
+ ,424.36
+ ,2304.00
+ ,988.80
+ ,9.40
+ ,29.82
+ ,1168.00
+ ,1167.00
+ ,18.30
+ ,67.00
+ ,334.89
+ ,4489.00
+ ,1226.10
+ ,12.20
+ ,29.83
+ ,1160.00
+ ,1166.00
+ ,16.70
+ ,74.00
+ ,278.89
+ ,5476.00
+ ,1235.80
+ ,12.20
+ ,29.83
+ ,1147.00
+ ,1164.00
+ ,22.80
+ ,55.00
+ ,519.84
+ ,3025.00
+ ,1254.00
+ ,13.30
+ ,29.71
+ ,1161.00
+ ,1162.00
+ ,13.90
+ ,47.00
+ ,193.21
+ ,2209.00
+ ,653.30
+ ,2.80
+ ,29.98
+ ,1143.00
+ ,1161.00
+ ,10.00
+ ,65.00
+ ,100.00
+ ,4225.00
+ ,650.00
+ ,3.90
+ ,30.18
+ ,1161.00
+ ,1159.00
+ ,16.10
+ ,28.00
+ ,259.21
+ ,784.00
+ ,450.80
+ ,-2.20
+ ,30.88
+ ,1161.00
+ ,1158.00
+ ,20.60
+ ,30.00
+ ,424.36
+ ,900.00
+ ,618.00
+ ,5.00
+ ,30.13
+ ,1168.00
+ ,1156.00
+ ,19.40
+ ,67.00
+ ,376.36
+ ,4489.00
+ ,1299.80
+ ,13.30
+ ,30.24
+ ,1172.00
+ ,1155.00
+ ,25.60
+ ,32.00
+ ,655.36
+ ,1024.00
+ ,819.20
+ ,7.80
+ ,30.24)
+ ,dim=c(9
+ ,34)
+ ,dimnames=list(c('15thbird'
+ ,'Sunset'
+ ,'Temp'
+ ,'humidity'
+ ,'Temp^2'
+ ,'Hum^2'
+ ,'TxH'
+ ,'Dew'
+ ,'pressure')
+ ,1:34))
> y <- array(NA,dim=c(9,34),dimnames=list(c('15thbird','Sunset','Temp','humidity','Temp^2','Hum^2','TxH','Dew','pressure'),1:34))
> 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 = 'Do not include Seasonal 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
> 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
15thbird Sunset Temp humidity Temp^2 Hum^2 TxH Dew pressure
1 1217 1210 31.0 48 961.00 2304 1488.0 19.00 30.00
2 1202 1209 34.4 38 1183.36 1444 1307.2 18.30 29.95
3 1180 1207 35.6 37 1267.36 1369 1317.2 18.90 29.94
4 1167 1206 32.8 48 1075.84 2304 1574.4 20.60 29.83
5 1186 1204 23.3 81 542.89 6561 1887.3 20.00 29.85
6 1168 1201 20.0 58 400.00 3364 1160.0 11.76 29.92
7 1142 1199 16.7 93 278.89 8649 1553.1 15.60 29.95
8 1147 1198 17.8 86 316.84 7396 1530.8 15.60 29.94
9 1183 1196 21.2 68 449.44 4624 1441.6 15.80 29.94
10 1149 1195 23.9 68 571.21 4624 1625.2 17.80 30.00
11 1197 1193 28.8 68 829.44 4624 1958.4 16.70 30.03
12 1210 1191 25.6 59 655.36 3481 1510.4 17.20 29.99
13 1206 1190 29.4 43 864.36 1849 1264.2 15.60 29.89
14 1196 1188 22.8 59 519.84 3481 1345.2 14.40 29.98
15 1190 1187 16.1 31 259.21 961 499.1 -0.60 30.26
16 1175 1185 16.1 49 259.21 2401 788.9 5.60 30.26
17 1186 1183 20.0 52 400.00 2704 1040.0 10.08 30.23
18 1172 1182 20.6 75 424.36 5625 1545.0 16.10 30.16
19 1152 1185 18.3 90 334.89 8100 1647.0 16.70 30.00
20 1154 1179 21.6 86 466.56 7396 1857.6 18.30 30.60
21 1168 1177 22.8 87 519.84 7569 1983.6 20.60 30.00
22 1180 1175 22.8 47 519.84 2209 1071.6 11.10 30.06
23 1169 1174 17.2 70 295.84 4900 1204.0 11.70 30.01
24 1166 1170 22.2 61 492.84 3721 1354.2 14.40 29.86
25 1177 1169 20.6 48 424.36 2304 988.8 9.40 29.82
26 1168 1167 18.3 67 334.89 4489 1226.1 12.20 29.83
27 1160 1166 16.7 74 278.89 5476 1235.8 12.20 29.83
28 1147 1164 22.8 55 519.84 3025 1254.0 13.30 29.71
29 1161 1162 13.9 47 193.21 2209 653.3 2.80 29.98
30 1143 1161 10.0 65 100.00 4225 650.0 3.90 30.18
31 1161 1159 16.1 28 259.21 784 450.8 -2.20 30.88
32 1161 1158 20.6 30 424.36 900 618.0 5.00 30.13
33 1168 1156 19.4 67 376.36 4489 1299.8 13.30 30.24
34 1172 1155 25.6 32 655.36 1024 819.2 7.80 30.24
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Sunset Temp humidity `Temp^2` `Hum^2`
254.49552 0.48361 6.18207 1.89803 -0.09285 -0.01789
TxH Dew pressure
0.01499 -2.04193 7.40204
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-34.565 -6.371 2.085 7.110 22.882
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 254.495517 607.880416 0.419 0.6790
Sunset 0.483607 0.219877 2.199 0.0373 *
Temp 6.182072 7.450348 0.830 0.4145
humidity 1.898029 2.308048 0.822 0.4187
`Temp^2` -0.092852 0.089099 -1.042 0.3073
`Hum^2` -0.017885 0.009958 -1.796 0.0846 .
TxH 0.014994 0.048151 0.311 0.7581
Dew -2.041931 3.029348 -0.674 0.5065
pressure 7.402039 15.140577 0.489 0.6292
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.65 on 25 degrees of freedom
Multiple R-squared: 0.5736, Adjusted R-squared: 0.4371
F-statistic: 4.204 on 8 and 25 DF, p-value: 0.002664
> 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.9826199 0.03476020 0.01738010
[2,] 0.9878797 0.02424051 0.01212026
[3,] 0.9841154 0.03176915 0.01588457
[4,] 0.9744490 0.05110192 0.02555096
[5,] 0.9592613 0.08147734 0.04073867
[6,] 0.9303780 0.13924398 0.06962199
[7,] 0.8661924 0.26761516 0.13380758
[8,] 0.8019547 0.39609054 0.19804527
[9,] 0.8630095 0.27398096 0.13699048
[10,] 0.9229198 0.15416031 0.07708015
[11,] 0.8366389 0.32672225 0.16336113
> postscript(file="/var/wessaorg/rcomp/tmp/1fad21331049149.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/2juvf1331049149.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/3900g1331049149.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/48k0i1331049149.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/5g2ak1331049149.ps",horizontal=F,onefile=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 = 34
Frequency = 1
1 2 3 4 5 6
19.4544477 9.8159747 -9.1298566 -25.8457223 10.8044911 -18.5755496
7 8 9 10 11 12
-4.6340652 -11.1429322 3.4492002 -34.5652623 0.6229646 22.8824271
13 14 15 16 17 18
17.6241441 11.8932797 11.6546356 -2.4727420 3.7872847 2.6511738
19 20 21 22 23 24
3.7879193 -8.8164806 12.1240433 -1.0221100 6.3675043 -6.9497967
25 26 27 28 29 30
2.9623348 5.9434778 7.3400134 -24.9439752 0.8728755 1.5196252
31 32 33 34
-8.4081737 -4.3841952 6.4197164 -1.0866716
> postscript(file="/var/wessaorg/rcomp/tmp/6r3i41331049149.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 34
Frequency = 1
lag(myerror, k = 1) myerror
0 19.4544477 NA
1 9.8159747 19.4544477
2 -9.1298566 9.8159747
3 -25.8457223 -9.1298566
4 10.8044911 -25.8457223
5 -18.5755496 10.8044911
6 -4.6340652 -18.5755496
7 -11.1429322 -4.6340652
8 3.4492002 -11.1429322
9 -34.5652623 3.4492002
10 0.6229646 -34.5652623
11 22.8824271 0.6229646
12 17.6241441 22.8824271
13 11.8932797 17.6241441
14 11.6546356 11.8932797
15 -2.4727420 11.6546356
16 3.7872847 -2.4727420
17 2.6511738 3.7872847
18 3.7879193 2.6511738
19 -8.8164806 3.7879193
20 12.1240433 -8.8164806
21 -1.0221100 12.1240433
22 6.3675043 -1.0221100
23 -6.9497967 6.3675043
24 2.9623348 -6.9497967
25 5.9434778 2.9623348
26 7.3400134 5.9434778
27 -24.9439752 7.3400134
28 0.8728755 -24.9439752
29 1.5196252 0.8728755
30 -8.4081737 1.5196252
31 -4.3841952 -8.4081737
32 6.4197164 -4.3841952
33 -1.0866716 6.4197164
34 NA -1.0866716
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 9.8159747 19.4544477
[2,] -9.1298566 9.8159747
[3,] -25.8457223 -9.1298566
[4,] 10.8044911 -25.8457223
[5,] -18.5755496 10.8044911
[6,] -4.6340652 -18.5755496
[7,] -11.1429322 -4.6340652
[8,] 3.4492002 -11.1429322
[9,] -34.5652623 3.4492002
[10,] 0.6229646 -34.5652623
[11,] 22.8824271 0.6229646
[12,] 17.6241441 22.8824271
[13,] 11.8932797 17.6241441
[14,] 11.6546356 11.8932797
[15,] -2.4727420 11.6546356
[16,] 3.7872847 -2.4727420
[17,] 2.6511738 3.7872847
[18,] 3.7879193 2.6511738
[19,] -8.8164806 3.7879193
[20,] 12.1240433 -8.8164806
[21,] -1.0221100 12.1240433
[22,] 6.3675043 -1.0221100
[23,] -6.9497967 6.3675043
[24,] 2.9623348 -6.9497967
[25,] 5.9434778 2.9623348
[26,] 7.3400134 5.9434778
[27,] -24.9439752 7.3400134
[28,] 0.8728755 -24.9439752
[29,] 1.5196252 0.8728755
[30,] -8.4081737 1.5196252
[31,] -4.3841952 -8.4081737
[32,] 6.4197164 -4.3841952
[33,] -1.0866716 6.4197164
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 9.8159747 19.4544477
2 -9.1298566 9.8159747
3 -25.8457223 -9.1298566
4 10.8044911 -25.8457223
5 -18.5755496 10.8044911
6 -4.6340652 -18.5755496
7 -11.1429322 -4.6340652
8 3.4492002 -11.1429322
9 -34.5652623 3.4492002
10 0.6229646 -34.5652623
11 22.8824271 0.6229646
12 17.6241441 22.8824271
13 11.8932797 17.6241441
14 11.6546356 11.8932797
15 -2.4727420 11.6546356
16 3.7872847 -2.4727420
17 2.6511738 3.7872847
18 3.7879193 2.6511738
19 -8.8164806 3.7879193
20 12.1240433 -8.8164806
21 -1.0221100 12.1240433
22 6.3675043 -1.0221100
23 -6.9497967 6.3675043
24 2.9623348 -6.9497967
25 5.9434778 2.9623348
26 7.3400134 5.9434778
27 -24.9439752 7.3400134
28 0.8728755 -24.9439752
29 1.5196252 0.8728755
30 -8.4081737 1.5196252
31 -4.3841952 -8.4081737
32 6.4197164 -4.3841952
33 -1.0866716 6.4197164
> 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/wessaorg/rcomp/tmp/77r1g1331049149.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/8lc501331049149.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/9mi361331049149.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/10cymm1331049149.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11zm271331049149.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/wessaorg/rcomp/tmp/1263qn1331049149.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/wessaorg/rcomp/tmp/13507w1331049149.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/wessaorg/rcomp/tmp/143q1a1331049149.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/wessaorg/rcomp/tmp/15h3r31331049149.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/wessaorg/rcomp/tmp/16fw421331049149.tab")
+ }
>
> try(system("convert tmp/1fad21331049149.ps tmp/1fad21331049149.png",intern=TRUE))
character(0)
> try(system("convert tmp/2juvf1331049149.ps tmp/2juvf1331049149.png",intern=TRUE))
character(0)
> try(system("convert tmp/3900g1331049149.ps tmp/3900g1331049149.png",intern=TRUE))
character(0)
> try(system("convert tmp/48k0i1331049149.ps tmp/48k0i1331049149.png",intern=TRUE))
character(0)
> try(system("convert tmp/5g2ak1331049149.ps tmp/5g2ak1331049149.png",intern=TRUE))
character(0)
> try(system("convert tmp/6r3i41331049149.ps tmp/6r3i41331049149.png",intern=TRUE))
character(0)
> try(system("convert tmp/77r1g1331049149.ps tmp/77r1g1331049149.png",intern=TRUE))
character(0)
> try(system("convert tmp/8lc501331049149.ps tmp/8lc501331049149.png",intern=TRUE))
character(0)
> try(system("convert tmp/9mi361331049149.ps tmp/9mi361331049149.png",intern=TRUE))
character(0)
> try(system("convert tmp/10cymm1331049149.ps tmp/10cymm1331049149.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.194 0.662 3.879