R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(-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/rcomp/tmp/13g9i1292959398.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/rcomp/tmp/23g9i1292959398.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/rcomp/tmp/33g9i1292959398.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/rcomp/tmp/4vpqk1292959398.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/rcomp/tmp/5vpqk1292959398.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/rcomp/tmp/6vpqk1292959398.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/rcomp/tmp/7oz8o1292959398.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/rcomp/tmp/8h7781292959398.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/rcomp/tmp/9h7781292959398.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/rcomp/tmp/10h7781292959398.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/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/11dh4h1292959398.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/125rm21292959398.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/135k3r1292959399.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/14gt3u1292959399.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/151cji1292959399.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/16flz81292959399.tab")
+ }
>
> try(system("convert tmp/13g9i1292959398.ps tmp/13g9i1292959398.png",intern=TRUE))
character(0)
> try(system("convert tmp/23g9i1292959398.ps tmp/23g9i1292959398.png",intern=TRUE))
character(0)
> try(system("convert tmp/33g9i1292959398.ps tmp/33g9i1292959398.png",intern=TRUE))
character(0)
> try(system("convert tmp/4vpqk1292959398.ps tmp/4vpqk1292959398.png",intern=TRUE))
character(0)
> try(system("convert tmp/5vpqk1292959398.ps tmp/5vpqk1292959398.png",intern=TRUE))
character(0)
> try(system("convert tmp/6vpqk1292959398.ps tmp/6vpqk1292959398.png",intern=TRUE))
character(0)
> try(system("convert tmp/7oz8o1292959398.ps tmp/7oz8o1292959398.png",intern=TRUE))
character(0)
> try(system("convert tmp/8h7781292959398.ps tmp/8h7781292959398.png",intern=TRUE))
character(0)
> try(system("convert tmp/9h7781292959398.ps tmp/9h7781292959398.png",intern=TRUE))
character(0)
> try(system("convert tmp/10h7781292959398.ps tmp/10h7781292959398.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.599 1.680 6.490