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.
Natural language support but running in an English locale
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(-999.00
+ ,6654.00
+ ,3.00
+ ,6.30
+ ,1.00
+ ,3.00
+ ,-999.00
+ ,3.39
+ ,1.00
+ ,-999.00
+ ,0.92
+ ,3.00
+ ,2.10
+ ,2547.00
+ ,4.00
+ ,9.10
+ ,10.55
+ ,4.00
+ ,15.80
+ ,0.02
+ ,1.00
+ ,5.20
+ ,160.00
+ ,4.00
+ ,10.90
+ ,3.30
+ ,1.00
+ ,8.30
+ ,52.16
+ ,1.00
+ ,11.00
+ ,0.43
+ ,4.00
+ ,3.20
+ ,465.00
+ ,5.00
+ ,7.60
+ ,0.55
+ ,2.00
+ ,-999.00
+ ,187.10
+ ,5.00
+ ,6.30
+ ,0.08
+ ,1.00
+ ,8.60
+ ,3.00
+ ,2.00
+ ,6.60
+ ,0.79
+ ,2.00
+ ,9.50
+ ,0.20
+ ,2.00
+ ,4.80
+ ,1.41
+ ,1.00
+ ,12.00
+ ,60.00
+ ,1.00
+ ,-999.00
+ ,529.00
+ ,5.00
+ ,3.30
+ ,27.66
+ ,5.00
+ ,11.00
+ ,0.12
+ ,2.00
+ ,-999.00
+ ,207.00
+ ,1.00
+ ,4.70
+ ,85.00
+ ,1.00
+ ,-999.00
+ ,36.33
+ ,1.00
+ ,10.40
+ ,0.10
+ ,3.00
+ ,7.40
+ ,1.04
+ ,4.00
+ ,2.10
+ ,521.00
+ ,5.00
+ ,-999.00
+ ,100.00
+ ,1.00
+ ,-999.00
+ ,35.00
+ ,4.00
+ ,7.70
+ ,0.01
+ ,4.00
+ ,17.90
+ ,0.01
+ ,1.00
+ ,6.10
+ ,62.00
+ ,1.00
+ ,8.20
+ ,0.12
+ ,1.00
+ ,8.40
+ ,1.35
+ ,3.00
+ ,11.90
+ ,0.02
+ ,3.00
+ ,10.80
+ ,0.05
+ ,3.00
+ ,13.80
+ ,1.70
+ ,1.00
+ ,14.30
+ ,3.50
+ ,1.00
+ ,-999.00
+ ,250.00
+ ,5.00
+ ,15.20
+ ,0.48
+ ,2.00
+ ,10.00
+ ,10.00
+ ,4.00
+ ,11.90
+ ,1.62
+ ,2.00
+ ,6.50
+ ,192.00
+ ,4.00
+ ,7.50
+ ,2.50
+ ,5.00
+ ,-999.00
+ ,4.29
+ ,2.00
+ ,10.60
+ ,0.28
+ ,3.00
+ ,7.40
+ ,4.24
+ ,1.00
+ ,8.40
+ ,6.80
+ ,2.00
+ ,5.70
+ ,0.75
+ ,2.00
+ ,4.90
+ ,3.60
+ ,3.00
+ ,-999.00
+ ,14.83
+ ,5.00
+ ,3.20
+ ,55.50
+ ,5.00
+ ,-999.00
+ ,1.40
+ ,2.00
+ ,8.10
+ ,0.06
+ ,2.00
+ ,11.00
+ ,0.90
+ ,2.00
+ ,4.90
+ ,2.00
+ ,3.00
+ ,13.20
+ ,0.10
+ ,2.00
+ ,9.70
+ ,4.19
+ ,4.00
+ ,12.80
+ ,3.50
+ ,1.00
+ ,-999.00
+ ,4.05
+ ,1.00)
+ ,dim=c(3
+ ,62)
+ ,dimnames=list(c('SWS'
+ ,'Wb'
+ ,'D')
+ ,1:62))
> y <- array(NA,dim=c(3,62),dimnames=list(c('SWS','Wb','D'),1:62))
> 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
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
SWS Wb D
1 -999.0 6654.00 3
2 6.3 1.00 3
3 -999.0 3.39 1
4 -999.0 0.92 3
5 2.1 2547.00 4
6 9.1 10.55 4
7 15.8 0.02 1
8 5.2 160.00 4
9 10.9 3.30 1
10 8.3 52.16 1
11 11.0 0.43 4
12 3.2 465.00 5
13 7.6 0.55 2
14 -999.0 187.10 5
15 6.3 0.08 1
16 8.6 3.00 2
17 6.6 0.79 2
18 9.5 0.20 2
19 4.8 1.41 1
20 12.0 60.00 1
21 -999.0 529.00 5
22 3.3 27.66 5
23 11.0 0.12 2
24 -999.0 207.00 1
25 4.7 85.00 1
26 -999.0 36.33 1
27 10.4 0.10 3
28 7.4 1.04 4
29 2.1 521.00 5
30 -999.0 100.00 1
31 -999.0 35.00 4
32 7.7 0.01 4
33 17.9 0.01 1
34 6.1 62.00 1
35 8.2 0.12 1
36 8.4 1.35 3
37 11.9 0.02 3
38 10.8 0.05 3
39 13.8 1.70 1
40 14.3 3.50 1
41 -999.0 250.00 5
42 15.2 0.48 2
43 10.0 10.00 4
44 11.9 1.62 2
45 6.5 192.00 4
46 7.5 2.50 5
47 -999.0 4.29 2
48 10.6 0.28 3
49 7.4 4.24 1
50 8.4 6.80 2
51 5.7 0.75 2
52 4.9 3.60 3
53 -999.0 14.83 5
54 3.2 55.50 5
55 -999.0 1.40 2
56 8.1 0.06 2
57 11.0 0.90 2
58 4.9 2.00 3
59 13.2 0.10 2
60 9.7 4.19 4
61 12.8 3.50 1
62 -999.0 4.05 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Wb D
-168.2383 -0.1052 -11.3715
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-819.0 186.3 198.9 213.1 483.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -168.23833 111.19493 -1.513 0.1356
Wb -0.10521 0.06038 -1.743 0.0866 .
D -11.37149 37.66918 -0.302 0.7638
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 420.2 on 59 degrees of freedom
Multiple R-squared: 0.05339, Adjusted R-squared: 0.0213
F-statistic: 1.664 on 2 and 59 DF, p-value: 0.1982
> 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.7111452 0.5777095 0.2888548
[2,] 0.8966049 0.2067903 0.1033951
[3,] 0.8284958 0.3430085 0.1715042
[4,] 0.8305627 0.3388747 0.1694373
[5,] 0.7969649 0.4060701 0.2030351
[6,] 0.7142138 0.5715724 0.2857862
[7,] 0.6601647 0.6796707 0.3398353
[8,] 0.5869421 0.8261158 0.4130579
[9,] 0.8131893 0.3736214 0.1868107
[10,] 0.7595781 0.4808438 0.2404219
[11,] 0.6980832 0.6038336 0.3019168
[12,] 0.6299790 0.7400421 0.3700210
[13,] 0.5582600 0.8834800 0.4417400
[14,] 0.4837045 0.9674091 0.5162955
[15,] 0.4172348 0.8344696 0.5827652
[16,] 0.4916752 0.9833504 0.5083248
[17,] 0.4449937 0.8899874 0.5550063
[18,] 0.3784249 0.7568497 0.6215751
[19,] 0.5646510 0.8706981 0.4353490
[20,] 0.5035850 0.9928300 0.4964150
[21,] 0.6858012 0.6283977 0.3141988
[22,] 0.6314504 0.7370993 0.3685496
[23,] 0.5746948 0.8506104 0.4253052
[24,] 0.6264888 0.7470223 0.3735112
[25,] 0.7463673 0.5072653 0.2536327
[26,] 0.8698283 0.2603433 0.1301717
[27,] 0.8340934 0.3318132 0.1659066
[28,] 0.7934140 0.4131720 0.2065860
[29,] 0.7545153 0.4909695 0.2454847
[30,] 0.7022372 0.5955256 0.2977628
[31,] 0.6443478 0.7113044 0.3556522
[32,] 0.5828399 0.8343203 0.4171601
[33,] 0.5187773 0.9624454 0.4812227
[34,] 0.4555867 0.9111734 0.5444133
[35,] 0.3951068 0.7902136 0.6048932
[36,] 0.5097246 0.9805508 0.4902754
[37,] 0.4484198 0.8968396 0.5515802
[38,] 0.3856126 0.7712253 0.6143874
[39,] 0.3277576 0.6555152 0.6722424
[40,] 0.2623900 0.5247800 0.7376100
[41,] 0.2112491 0.4224981 0.7887509
[42,] 0.3748012 0.7496024 0.6251988
[43,] 0.3122735 0.6245469 0.6877265
[44,] 0.2435958 0.4871916 0.7564042
[45,] 0.1875752 0.3751504 0.8124248
[46,] 0.1417378 0.2834756 0.8582622
[47,] 0.1058226 0.2116453 0.8941774
[48,] 0.3116561 0.6233122 0.6883439
[49,] 0.2540865 0.5081729 0.7459135
[50,] 0.6399281 0.7201438 0.3600719
[51,] 0.4716334 0.9432667 0.5283666
> postscript(file="/var/www/html/freestat/rcomp/tmp/1cn1t1292966783.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/html/freestat/rcomp/tmp/2cn1t1292966783.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/html/freestat/rcomp/tmp/35x1w1292966783.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/html/freestat/rcomp/tmp/45x1w1292966783.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/html/freestat/rcomp/tmp/55x1w1292966783.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 = 62
Frequency = 1
1 2 3 4 5 6 7
-96.56439 208.75802 -819.03351 -796.55040 483.80008 223.93428 195.41193
8 9 10 11 12 13 14
235.75827 190.85703 193.39770 224.76954 277.21951 198.63918 -754.21899
15 16 17 18 19 20 21
185.91824 199.89695 197.66443 200.50236 184.55817 197.92256 -718.24690
22 23 24 25 26 27 28
231.30596 201.99394 -797.61122 193.25287 -815.56781 212.76333 221.23371
29 30 31 32 33 34 35
282.01140 -808.86894 -781.59327 221.42535 197.51088 192.23299 187.82245
36 37 38 39 40 41 42
210.89484 214.25491 213.15806 193.58869 194.27807 -747.60114 206.23182
43 44 45 46 47 48 49
224.77642 203.05176 240.42506 232.85881 -807.56732 212.98226 187.45592
50 51 52 53 54 55 56
200.09676 196.76022 207.63157 -772.34392 234.13507 -807.87139 199.08763
57 58 59 60 61 62
202.07601 207.46323 204.19184 223.86513 192.77807 -818.96407
> postscript(file="/var/www/html/freestat/rcomp/tmp/6f6ih1292966783.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 = 62
Frequency = 1
lag(myerror, k = 1) myerror
0 -96.56439 NA
1 208.75802 -96.56439
2 -819.03351 208.75802
3 -796.55040 -819.03351
4 483.80008 -796.55040
5 223.93428 483.80008
6 195.41193 223.93428
7 235.75827 195.41193
8 190.85703 235.75827
9 193.39770 190.85703
10 224.76954 193.39770
11 277.21951 224.76954
12 198.63918 277.21951
13 -754.21899 198.63918
14 185.91824 -754.21899
15 199.89695 185.91824
16 197.66443 199.89695
17 200.50236 197.66443
18 184.55817 200.50236
19 197.92256 184.55817
20 -718.24690 197.92256
21 231.30596 -718.24690
22 201.99394 231.30596
23 -797.61122 201.99394
24 193.25287 -797.61122
25 -815.56781 193.25287
26 212.76333 -815.56781
27 221.23371 212.76333
28 282.01140 221.23371
29 -808.86894 282.01140
30 -781.59327 -808.86894
31 221.42535 -781.59327
32 197.51088 221.42535
33 192.23299 197.51088
34 187.82245 192.23299
35 210.89484 187.82245
36 214.25491 210.89484
37 213.15806 214.25491
38 193.58869 213.15806
39 194.27807 193.58869
40 -747.60114 194.27807
41 206.23182 -747.60114
42 224.77642 206.23182
43 203.05176 224.77642
44 240.42506 203.05176
45 232.85881 240.42506
46 -807.56732 232.85881
47 212.98226 -807.56732
48 187.45592 212.98226
49 200.09676 187.45592
50 196.76022 200.09676
51 207.63157 196.76022
52 -772.34392 207.63157
53 234.13507 -772.34392
54 -807.87139 234.13507
55 199.08763 -807.87139
56 202.07601 199.08763
57 207.46323 202.07601
58 204.19184 207.46323
59 223.86513 204.19184
60 192.77807 223.86513
61 -818.96407 192.77807
62 NA -818.96407
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 208.7580 -96.56439
[2,] -819.0335 208.75802
[3,] -796.5504 -819.03351
[4,] 483.8001 -796.55040
[5,] 223.9343 483.80008
[6,] 195.4119 223.93428
[7,] 235.7583 195.41193
[8,] 190.8570 235.75827
[9,] 193.3977 190.85703
[10,] 224.7695 193.39770
[11,] 277.2195 224.76954
[12,] 198.6392 277.21951
[13,] -754.2190 198.63918
[14,] 185.9182 -754.21899
[15,] 199.8970 185.91824
[16,] 197.6644 199.89695
[17,] 200.5024 197.66443
[18,] 184.5582 200.50236
[19,] 197.9226 184.55817
[20,] -718.2469 197.92256
[21,] 231.3060 -718.24690
[22,] 201.9939 231.30596
[23,] -797.6112 201.99394
[24,] 193.2529 -797.61122
[25,] -815.5678 193.25287
[26,] 212.7633 -815.56781
[27,] 221.2337 212.76333
[28,] 282.0114 221.23371
[29,] -808.8689 282.01140
[30,] -781.5933 -808.86894
[31,] 221.4253 -781.59327
[32,] 197.5109 221.42535
[33,] 192.2330 197.51088
[34,] 187.8225 192.23299
[35,] 210.8948 187.82245
[36,] 214.2549 210.89484
[37,] 213.1581 214.25491
[38,] 193.5887 213.15806
[39,] 194.2781 193.58869
[40,] -747.6011 194.27807
[41,] 206.2318 -747.60114
[42,] 224.7764 206.23182
[43,] 203.0518 224.77642
[44,] 240.4251 203.05176
[45,] 232.8588 240.42506
[46,] -807.5673 232.85881
[47,] 212.9823 -807.56732
[48,] 187.4559 212.98226
[49,] 200.0968 187.45592
[50,] 196.7602 200.09676
[51,] 207.6316 196.76022
[52,] -772.3439 207.63157
[53,] 234.1351 -772.34392
[54,] -807.8714 234.13507
[55,] 199.0876 -807.87139
[56,] 202.0760 199.08763
[57,] 207.4632 202.07601
[58,] 204.1918 207.46323
[59,] 223.8651 204.19184
[60,] 192.7781 223.86513
[61,] -818.9641 192.77807
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 208.7580 -96.56439
2 -819.0335 208.75802
3 -796.5504 -819.03351
4 483.8001 -796.55040
5 223.9343 483.80008
6 195.4119 223.93428
7 235.7583 195.41193
8 190.8570 235.75827
9 193.3977 190.85703
10 224.7695 193.39770
11 277.2195 224.76954
12 198.6392 277.21951
13 -754.2190 198.63918
14 185.9182 -754.21899
15 199.8970 185.91824
16 197.6644 199.89695
17 200.5024 197.66443
18 184.5582 200.50236
19 197.9226 184.55817
20 -718.2469 197.92256
21 231.3060 -718.24690
22 201.9939 231.30596
23 -797.6112 201.99394
24 193.2529 -797.61122
25 -815.5678 193.25287
26 212.7633 -815.56781
27 221.2337 212.76333
28 282.0114 221.23371
29 -808.8689 282.01140
30 -781.5933 -808.86894
31 221.4253 -781.59327
32 197.5109 221.42535
33 192.2330 197.51088
34 187.8225 192.23299
35 210.8948 187.82245
36 214.2549 210.89484
37 213.1581 214.25491
38 193.5887 213.15806
39 194.2781 193.58869
40 -747.6011 194.27807
41 206.2318 -747.60114
42 224.7764 206.23182
43 203.0518 224.77642
44 240.4251 203.05176
45 232.8588 240.42506
46 -807.5673 232.85881
47 212.9823 -807.56732
48 187.4559 212.98226
49 200.0968 187.45592
50 196.7602 200.09676
51 207.6316 196.76022
52 -772.3439 207.63157
53 234.1351 -772.34392
54 -807.8714 234.13507
55 199.0876 -807.87139
56 202.0760 199.08763
57 207.4632 202.07601
58 204.1918 207.46323
59 223.8651 204.19184
60 192.7781 223.86513
61 -818.9641 192.77807
> 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/freestat/rcomp/tmp/7qxh21292966783.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/html/freestat/rcomp/tmp/8qxh21292966783.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/html/freestat/rcomp/tmp/9qxh21292966783.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/html/freestat/rcomp/tmp/1016g51292966783.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/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/1147fb1292966783.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/freestat/rcomp/tmp/12q7dh1292966783.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/freestat/rcomp/tmp/134zbp1292966783.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/freestat/rcomp/tmp/147iad1292966783.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/freestat/rcomp/tmp/15b0q11292966783.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/freestat/rcomp/tmp/16wjp71292966783.tab")
+ }
>
> try(system("convert tmp/1cn1t1292966783.ps tmp/1cn1t1292966783.png",intern=TRUE))
character(0)
> try(system("convert tmp/2cn1t1292966783.ps tmp/2cn1t1292966783.png",intern=TRUE))
character(0)
> try(system("convert tmp/35x1w1292966783.ps tmp/35x1w1292966783.png",intern=TRUE))
character(0)
> try(system("convert tmp/45x1w1292966783.ps tmp/45x1w1292966783.png",intern=TRUE))
character(0)
> try(system("convert tmp/55x1w1292966783.ps tmp/55x1w1292966783.png",intern=TRUE))
character(0)
> try(system("convert tmp/6f6ih1292966783.ps tmp/6f6ih1292966783.png",intern=TRUE))
character(0)
> try(system("convert tmp/7qxh21292966783.ps tmp/7qxh21292966783.png",intern=TRUE))
character(0)
> try(system("convert tmp/8qxh21292966783.ps tmp/8qxh21292966783.png",intern=TRUE))
character(0)
> try(system("convert tmp/9qxh21292966783.ps tmp/9qxh21292966783.png",intern=TRUE))
character(0)
> try(system("convert tmp/1016g51292966783.ps tmp/1016g51292966783.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.896 2.511 4.231