R version 2.12.0 (2010-10-15)
Copyright (C) 2010 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
+ ,10.00
+ ,1202.00
+ ,1209.00
+ ,34.40
+ ,38.00
+ ,1183.36
+ ,1444.00
+ ,1307.20
+ ,18.30
+ ,29.95
+ ,10.00
+ ,1180.00
+ ,1207.00
+ ,35.60
+ ,37.00
+ ,1267.36
+ ,1369.00
+ ,1317.20
+ ,18.90
+ ,29.94
+ ,10.00
+ ,1167.00
+ ,1206.00
+ ,32.80
+ ,48.00
+ ,1075.84
+ ,2304.00
+ ,1574.40
+ ,20.60
+ ,29.83
+ ,10.00
+ ,1186.00
+ ,1204.00
+ ,23.30
+ ,81.00
+ ,542.89
+ ,6561.00
+ ,1887.30
+ ,20.00
+ ,29.85
+ ,9.00
+ ,1168.00
+ ,1201.00
+ ,20.00
+ ,58.00
+ ,400.00
+ ,3364.00
+ ,1160.00
+ ,11.76
+ ,29.92
+ ,10.00
+ ,1142.00
+ ,1199.00
+ ,16.70
+ ,93.00
+ ,278.89
+ ,8649.00
+ ,1553.10
+ ,15.60
+ ,29.95
+ ,6.00
+ ,1147.00
+ ,1198.00
+ ,17.80
+ ,86.00
+ ,316.84
+ ,7396.00
+ ,1530.80
+ ,15.60
+ ,29.94
+ ,10.00
+ ,1183.00
+ ,1196.00
+ ,21.20
+ ,68.00
+ ,449.44
+ ,4624.00
+ ,1441.60
+ ,15.80
+ ,29.94
+ ,10.00
+ ,1149.00
+ ,1195.00
+ ,23.90
+ ,68.00
+ ,571.21
+ ,4624.00
+ ,1625.20
+ ,17.80
+ ,30.00
+ ,10.00
+ ,1197.00
+ ,1193.00
+ ,28.80
+ ,68.00
+ ,829.44
+ ,4624.00
+ ,1958.40
+ ,16.70
+ ,30.03
+ ,10.00
+ ,1210.00
+ ,1191.00
+ ,25.60
+ ,59.00
+ ,655.36
+ ,3481.00
+ ,1510.40
+ ,17.20
+ ,29.99
+ ,10.00
+ ,1206.00
+ ,1190.00
+ ,29.40
+ ,43.00
+ ,864.36
+ ,1849.00
+ ,1264.20
+ ,15.60
+ ,29.89
+ ,10.00
+ ,1196.00
+ ,1188.00
+ ,22.80
+ ,59.00
+ ,519.84
+ ,3481.00
+ ,1345.20
+ ,14.40
+ ,29.98
+ ,6.00
+ ,1190.00
+ ,1187.00
+ ,16.10
+ ,31.00
+ ,259.21
+ ,961.00
+ ,499.10
+ ,-0.60
+ ,30.26
+ ,10.00
+ ,1175.00
+ ,1185.00
+ ,16.10
+ ,49.00
+ ,259.21
+ ,2401.00
+ ,788.90
+ ,5.60
+ ,30.26
+ ,10.00
+ ,1186.00
+ ,1183.00
+ ,20.00
+ ,52.00
+ ,400.00
+ ,2704.00
+ ,1040.00
+ ,10.08
+ ,30.23
+ ,10.00
+ ,1172.00
+ ,1182.00
+ ,20.60
+ ,75.00
+ ,424.36
+ ,5625.00
+ ,1545.00
+ ,16.10
+ ,30.16
+ ,10.00
+ ,1152.00
+ ,1185.00
+ ,18.30
+ ,90.00
+ ,334.89
+ ,8100.00
+ ,1647.00
+ ,16.70
+ ,30.00
+ ,10.00
+ ,1154.00
+ ,1179.00
+ ,21.60
+ ,86.00
+ ,466.56
+ ,7396.00
+ ,1857.60
+ ,18.30
+ ,30.60
+ ,8.00
+ ,1168.00
+ ,1177.00
+ ,22.80
+ ,87.00
+ ,519.84
+ ,7569.00
+ ,1983.60
+ ,20.60
+ ,30.00
+ ,10.00
+ ,1180.00
+ ,1175.00
+ ,22.80
+ ,47.00
+ ,519.84
+ ,2209.00
+ ,1071.60
+ ,11.10
+ ,30.06
+ ,10.00
+ ,1169.00
+ ,1174.00
+ ,17.20
+ ,70.00
+ ,295.84
+ ,4900.00
+ ,1204.00
+ ,11.70
+ ,30.01
+ ,10.00
+ ,1166.00
+ ,1170.00
+ ,22.20
+ ,61.00
+ ,492.84
+ ,3721.00
+ ,1354.20
+ ,14.40
+ ,29.86
+ ,10.00
+ ,1177.00
+ ,1169.00
+ ,20.60
+ ,48.00
+ ,424.36
+ ,2304.00
+ ,988.80
+ ,9.40
+ ,29.82
+ ,10.00
+ ,1168.00
+ ,1167.00
+ ,18.30
+ ,67.00
+ ,334.89
+ ,4489.00
+ ,1226.10
+ ,12.20
+ ,29.83
+ ,10.00
+ ,1160.00
+ ,1166.00
+ ,16.70
+ ,74.00
+ ,278.89
+ ,5476.00
+ ,1235.80
+ ,12.20
+ ,29.83
+ ,10.00
+ ,1147.00
+ ,1164.00
+ ,22.80
+ ,55.00
+ ,519.84
+ ,3025.00
+ ,1254.00
+ ,13.30
+ ,29.71
+ ,10.00
+ ,1161.00
+ ,1162.00
+ ,13.90
+ ,47.00
+ ,193.21
+ ,2209.00
+ ,653.30
+ ,2.80
+ ,29.98
+ ,10.00
+ ,1143.00
+ ,1161.00
+ ,10.00
+ ,65.00
+ ,100.00
+ ,4225.00
+ ,650.00
+ ,3.90
+ ,30.18
+ ,10.00
+ ,1161.00
+ ,1159.00
+ ,16.10
+ ,28.00
+ ,259.21
+ ,784.00
+ ,450.80
+ ,-2.20
+ ,30.88
+ ,10.00
+ ,1161.00
+ ,1158.00
+ ,20.60
+ ,30.00
+ ,424.36
+ ,900.00
+ ,618.00
+ ,5.00
+ ,30.13
+ ,10.00
+ ,1168.00
+ ,1156.00
+ ,19.40
+ ,67.00
+ ,376.36
+ ,4489.00
+ ,1299.80
+ ,13.30
+ ,30.24
+ ,10.00
+ ,1172.00
+ ,1155.00
+ ,25.60
+ ,32.00
+ ,655.36
+ ,1024.00
+ ,819.20
+ ,7.80
+ ,30.24
+ ,10.00)
+ ,dim=c(10
+ ,34)
+ ,dimnames=list(c('15thbird'
+ ,'Sunset'
+ ,'Temp'
+ ,'humidity'
+ ,'Temp^2'
+ ,'Hum^2'
+ ,'TxH'
+ ,'Dew'
+ ,'pressure'
+ ,'visibility')
+ ,1:34))
> y <- array(NA,dim=c(10,34),dimnames=list(c('15thbird','Sunset','Temp','humidity','Temp^2','Hum^2','TxH','Dew','pressure','visibility'),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 visibility
1 1217 1210 31.0 48 961.00 2304 1488.0 19.00 30.00 10
2 1202 1209 34.4 38 1183.36 1444 1307.2 18.30 29.95 10
3 1180 1207 35.6 37 1267.36 1369 1317.2 18.90 29.94 10
4 1167 1206 32.8 48 1075.84 2304 1574.4 20.60 29.83 10
5 1186 1204 23.3 81 542.89 6561 1887.3 20.00 29.85 9
6 1168 1201 20.0 58 400.00 3364 1160.0 11.76 29.92 10
7 1142 1199 16.7 93 278.89 8649 1553.1 15.60 29.95 6
8 1147 1198 17.8 86 316.84 7396 1530.8 15.60 29.94 10
9 1183 1196 21.2 68 449.44 4624 1441.6 15.80 29.94 10
10 1149 1195 23.9 68 571.21 4624 1625.2 17.80 30.00 10
11 1197 1193 28.8 68 829.44 4624 1958.4 16.70 30.03 10
12 1210 1191 25.6 59 655.36 3481 1510.4 17.20 29.99 10
13 1206 1190 29.4 43 864.36 1849 1264.2 15.60 29.89 10
14 1196 1188 22.8 59 519.84 3481 1345.2 14.40 29.98 6
15 1190 1187 16.1 31 259.21 961 499.1 -0.60 30.26 10
16 1175 1185 16.1 49 259.21 2401 788.9 5.60 30.26 10
17 1186 1183 20.0 52 400.00 2704 1040.0 10.08 30.23 10
18 1172 1182 20.6 75 424.36 5625 1545.0 16.10 30.16 10
19 1152 1185 18.3 90 334.89 8100 1647.0 16.70 30.00 10
20 1154 1179 21.6 86 466.56 7396 1857.6 18.30 30.60 8
21 1168 1177 22.8 87 519.84 7569 1983.6 20.60 30.00 10
22 1180 1175 22.8 47 519.84 2209 1071.6 11.10 30.06 10
23 1169 1174 17.2 70 295.84 4900 1204.0 11.70 30.01 10
24 1166 1170 22.2 61 492.84 3721 1354.2 14.40 29.86 10
25 1177 1169 20.6 48 424.36 2304 988.8 9.40 29.82 10
26 1168 1167 18.3 67 334.89 4489 1226.1 12.20 29.83 10
27 1160 1166 16.7 74 278.89 5476 1235.8 12.20 29.83 10
28 1147 1164 22.8 55 519.84 3025 1254.0 13.30 29.71 10
29 1161 1162 13.9 47 193.21 2209 653.3 2.80 29.98 10
30 1143 1161 10.0 65 100.00 4225 650.0 3.90 30.18 10
31 1161 1159 16.1 28 259.21 784 450.8 -2.20 30.88 10
32 1161 1158 20.6 30 424.36 900 618.0 5.00 30.13 10
33 1168 1156 19.4 67 376.36 4489 1299.8 13.30 30.24 10
34 1172 1155 25.6 32 655.36 1024 819.2 7.80 30.24 10
> 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`
304.62446 0.46699 5.73065 1.84043 -0.08699 -0.01811
TxH Dew pressure visibility
0.01785 -1.97442 6.88851 -0.86755
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-34.206 -6.878 2.116 7.405 23.037
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 304.62446 642.32316 0.474 0.6396
Sunset 0.46699 0.23102 2.021 0.0545 .
Temp 5.73065 7.74392 0.740 0.4665
humidity 1.84043 2.35955 0.780 0.4430
`Temp^2` -0.08699 0.09294 -0.936 0.3586
`Hum^2` -0.01811 0.01017 -1.780 0.0877 .
TxH 0.01785 0.05000 0.357 0.7243
Dew -1.97442 3.09478 -0.638 0.5295
pressure 6.88851 15.52350 0.444 0.6612
visibility -0.86755 2.94953 -0.294 0.7712
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.92 on 24 degrees of freedom
Multiple R-squared: 0.5751, Adjusted R-squared: 0.4158
F-statistic: 3.61 on 9 and 24 DF, p-value: 0.005707
> 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.9909313 0.01813736 0.009068678
[2,] 0.9873782 0.02524364 0.012621820
[3,] 0.9835635 0.03287292 0.016436461
[4,] 0.9713261 0.05734786 0.028673931
[5,] 0.9412333 0.11753345 0.058766726
[6,] 0.8758306 0.24833874 0.124169369
[7,] 0.8641566 0.27168676 0.135843378
[8,] 0.7941626 0.41167489 0.205837446
[9,] 0.8469216 0.30615689 0.153078444
> postscript(file="/var/www/rcomp/tmp/1wq591331048937.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/www/rcomp/tmp/2129c1331048937.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/www/rcomp/tmp/3ia521331048937.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/www/rcomp/tmp/4fftq1331048937.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/www/rcomp/tmp/5o47x1331048937.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.6199992 9.9641490 -9.1144713 -26.0494221 10.5500948 -18.0379007
7 8 9 10 11 12
-6.5445507 -9.9512219 3.9480105 -34.2063734 0.7852219 23.0371087
13 14 15 16 17 18
17.7234961 8.7139322 11.8257016 -2.2207534 4.1473633 3.2188460
19 20 21 22 23 24
4.8966805 -9.6138104 12.5313410 -0.8790489 6.7635654 -6.9887943
25 26 27 28 29 30
2.8792071 6.0366783 7.6193938 -25.1886190 0.4105799 1.3530962
31 32 33 34
-8.3505492 -4.4874503 6.5095558 -0.9010557
> postscript(file="/var/www/rcomp/tmp/6l04c1331048937.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.6199992 NA
1 9.9641490 19.6199992
2 -9.1144713 9.9641490
3 -26.0494221 -9.1144713
4 10.5500948 -26.0494221
5 -18.0379007 10.5500948
6 -6.5445507 -18.0379007
7 -9.9512219 -6.5445507
8 3.9480105 -9.9512219
9 -34.2063734 3.9480105
10 0.7852219 -34.2063734
11 23.0371087 0.7852219
12 17.7234961 23.0371087
13 8.7139322 17.7234961
14 11.8257016 8.7139322
15 -2.2207534 11.8257016
16 4.1473633 -2.2207534
17 3.2188460 4.1473633
18 4.8966805 3.2188460
19 -9.6138104 4.8966805
20 12.5313410 -9.6138104
21 -0.8790489 12.5313410
22 6.7635654 -0.8790489
23 -6.9887943 6.7635654
24 2.8792071 -6.9887943
25 6.0366783 2.8792071
26 7.6193938 6.0366783
27 -25.1886190 7.6193938
28 0.4105799 -25.1886190
29 1.3530962 0.4105799
30 -8.3505492 1.3530962
31 -4.4874503 -8.3505492
32 6.5095558 -4.4874503
33 -0.9010557 6.5095558
34 NA -0.9010557
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 9.9641490 19.6199992
[2,] -9.1144713 9.9641490
[3,] -26.0494221 -9.1144713
[4,] 10.5500948 -26.0494221
[5,] -18.0379007 10.5500948
[6,] -6.5445507 -18.0379007
[7,] -9.9512219 -6.5445507
[8,] 3.9480105 -9.9512219
[9,] -34.2063734 3.9480105
[10,] 0.7852219 -34.2063734
[11,] 23.0371087 0.7852219
[12,] 17.7234961 23.0371087
[13,] 8.7139322 17.7234961
[14,] 11.8257016 8.7139322
[15,] -2.2207534 11.8257016
[16,] 4.1473633 -2.2207534
[17,] 3.2188460 4.1473633
[18,] 4.8966805 3.2188460
[19,] -9.6138104 4.8966805
[20,] 12.5313410 -9.6138104
[21,] -0.8790489 12.5313410
[22,] 6.7635654 -0.8790489
[23,] -6.9887943 6.7635654
[24,] 2.8792071 -6.9887943
[25,] 6.0366783 2.8792071
[26,] 7.6193938 6.0366783
[27,] -25.1886190 7.6193938
[28,] 0.4105799 -25.1886190
[29,] 1.3530962 0.4105799
[30,] -8.3505492 1.3530962
[31,] -4.4874503 -8.3505492
[32,] 6.5095558 -4.4874503
[33,] -0.9010557 6.5095558
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 9.9641490 19.6199992
2 -9.1144713 9.9641490
3 -26.0494221 -9.1144713
4 10.5500948 -26.0494221
5 -18.0379007 10.5500948
6 -6.5445507 -18.0379007
7 -9.9512219 -6.5445507
8 3.9480105 -9.9512219
9 -34.2063734 3.9480105
10 0.7852219 -34.2063734
11 23.0371087 0.7852219
12 17.7234961 23.0371087
13 8.7139322 17.7234961
14 11.8257016 8.7139322
15 -2.2207534 11.8257016
16 4.1473633 -2.2207534
17 3.2188460 4.1473633
18 4.8966805 3.2188460
19 -9.6138104 4.8966805
20 12.5313410 -9.6138104
21 -0.8790489 12.5313410
22 6.7635654 -0.8790489
23 -6.9887943 6.7635654
24 2.8792071 -6.9887943
25 6.0366783 2.8792071
26 7.6193938 6.0366783
27 -25.1886190 7.6193938
28 0.4105799 -25.1886190
29 1.3530962 0.4105799
30 -8.3505492 1.3530962
31 -4.4874503 -8.3505492
32 6.5095558 -4.4874503
33 -0.9010557 6.5095558
> 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/rcomp/tmp/710v81331048937.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/www/rcomp/tmp/8wlos1331048937.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/www/rcomp/tmp/9bu2d1331048937.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/www/rcomp/tmp/10wbqf1331048937.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11mvi21331048937.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/rcomp/tmp/12abtz1331048937.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/rcomp/tmp/13g0w11331048938.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/rcomp/tmp/1411d21331048938.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/rcomp/tmp/1518bx1331048938.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/rcomp/tmp/16uzq61331048938.tab")
+ }
>
> try(system("convert tmp/1wq591331048937.ps tmp/1wq591331048937.png",intern=TRUE))
character(0)
> try(system("convert tmp/2129c1331048937.ps tmp/2129c1331048937.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ia521331048937.ps tmp/3ia521331048937.png",intern=TRUE))
character(0)
> try(system("convert tmp/4fftq1331048937.ps tmp/4fftq1331048937.png",intern=TRUE))
character(0)
> try(system("convert tmp/5o47x1331048937.ps tmp/5o47x1331048937.png",intern=TRUE))
character(0)
> try(system("convert tmp/6l04c1331048937.ps tmp/6l04c1331048937.png",intern=TRUE))
character(0)
> try(system("convert tmp/710v81331048937.ps tmp/710v81331048937.png",intern=TRUE))
character(0)
> try(system("convert tmp/8wlos1331048937.ps tmp/8wlos1331048937.png",intern=TRUE))
character(0)
> try(system("convert tmp/9bu2d1331048937.ps tmp/9bu2d1331048937.png",intern=TRUE))
character(0)
> try(system("convert tmp/10wbqf1331048937.ps tmp/10wbqf1331048937.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.190 0.370 4.517