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(-999.000
+ ,6654.000
+ ,3.000
+ ,6.300
+ ,1.000
+ ,3.000
+ ,-999.000
+ ,3.385
+ ,1.000
+ ,-999.000
+ ,0.920
+ ,3.000
+ ,2.100
+ ,2547.000
+ ,4.000
+ ,9.100
+ ,10.550
+ ,4.000
+ ,15.800
+ ,0.023
+ ,1.000
+ ,5.200
+ ,160.000
+ ,4.000
+ ,10.900
+ ,3.300
+ ,1.000
+ ,8.300
+ ,52.160
+ ,1.000
+ ,11.000
+ ,0.425
+ ,4.000
+ ,3.200
+ ,465.000
+ ,5.000
+ ,7.600
+ ,0.550
+ ,2.000
+ ,-999.000
+ ,187.100
+ ,5.000
+ ,6.300
+ ,0.075
+ ,1.000
+ ,8.600
+ ,3.000
+ ,2.000
+ ,6.600
+ ,0.785
+ ,2.000
+ ,9.500
+ ,0.200
+ ,2.000
+ ,4.800
+ ,1.410
+ ,1.000
+ ,12.000
+ ,60.000
+ ,1.000
+ ,-999.000
+ ,529.000
+ ,5.000
+ ,3.300
+ ,27.660
+ ,5.000
+ ,11.000
+ ,0.120
+ ,2.000
+ ,-999.000
+ ,207.000
+ ,1.000
+ ,4.700
+ ,85.000
+ ,1.000
+ ,-999.000
+ ,36.330
+ ,1.000
+ ,10.400
+ ,0.101
+ ,3.000
+ ,7.400
+ ,1.040
+ ,4.000
+ ,2.100
+ ,521.000
+ ,5.000
+ ,-999.000
+ ,100.000
+ ,1.000
+ ,-999.000
+ ,35.000
+ ,4.000
+ ,7.700
+ ,0.005
+ ,4.000
+ ,17.900
+ ,0.010
+ ,1.000
+ ,6.100
+ ,62.000
+ ,1.000
+ ,8.200
+ ,0.122
+ ,1.000
+ ,8.400
+ ,1.350
+ ,3.000
+ ,11.900
+ ,0.023
+ ,3.000
+ ,10.800
+ ,0.048
+ ,3.000
+ ,13.800
+ ,1.700
+ ,1.000
+ ,14.300
+ ,3.500
+ ,1.000
+ ,-999.000
+ ,250.000
+ ,5.000
+ ,15.200
+ ,0.480
+ ,2.000
+ ,10.000
+ ,10.000
+ ,4.000
+ ,11.900
+ ,1.620
+ ,2.000
+ ,6.500
+ ,192.000
+ ,4.000
+ ,7.500
+ ,2.500
+ ,5.000
+ ,-999.000
+ ,4.288
+ ,2.000
+ ,10.600
+ ,0.280
+ ,3.000
+ ,7.400
+ ,4.235
+ ,1.000
+ ,8.400
+ ,6.800
+ ,2.000
+ ,5.700
+ ,0.750
+ ,2.000
+ ,4.900
+ ,3.600
+ ,3.000
+ ,-999.000
+ ,14.830
+ ,5.000
+ ,3.200
+ ,55.500
+ ,5.000
+ ,-999.000
+ ,1.400
+ ,2.000
+ ,8.100
+ ,0.060
+ ,2.000
+ ,11.000
+ ,0.900
+ ,2.000
+ ,4.900
+ ,2.000
+ ,3.000
+ ,13.200
+ ,0.104
+ ,2.000
+ ,9.700
+ ,4.190
+ ,4.000
+ ,12.800
+ ,3.500
+ ,1.000
+ ,-999.000
+ ,4.050
+ ,1.000)
+ ,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
> 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\r
1 -999.0 6654.000 3
2 6.3 1.000 3
3 -999.0 3.385 1
4 -999.0 0.920 3
5 2.1 2547.000 4
6 9.1 10.550 4
7 15.8 0.023 1
8 5.2 160.000 4
9 10.9 3.300 1
10 8.3 52.160 1
11 11.0 0.425 4
12 3.2 465.000 5
13 7.6 0.550 2
14 -999.0 187.100 5
15 6.3 0.075 1
16 8.6 3.000 2
17 6.6 0.785 2
18 9.5 0.200 2
19 4.8 1.410 1
20 12.0 60.000 1
21 -999.0 529.000 5
22 3.3 27.660 5
23 11.0 0.120 2
24 -999.0 207.000 1
25 4.7 85.000 1
26 -999.0 36.330 1
27 10.4 0.101 3
28 7.4 1.040 4
29 2.1 521.000 5
30 -999.0 100.000 1
31 -999.0 35.000 4
32 7.7 0.005 4
33 17.9 0.010 1
34 6.1 62.000 1
35 8.2 0.122 1
36 8.4 1.350 3
37 11.9 0.023 3
38 10.8 0.048 3
39 13.8 1.700 1
40 14.3 3.500 1
41 -999.0 250.000 5
42 15.2 0.480 2
43 10.0 10.000 4
44 11.9 1.620 2
45 6.5 192.000 4
46 7.5 2.500 5
47 -999.0 4.288 2
48 10.6 0.280 3
49 7.4 4.235 1
50 8.4 6.800 2
51 5.7 0.750 2
52 4.9 3.600 3
53 -999.0 14.830 5
54 3.2 55.500 5
55 -999.0 1.400 2
56 8.1 0.060 2
57 11.0 0.900 2
58 4.9 2.000 3
59 13.2 0.104 2
60 9.7 4.190 4
61 12.8 3.500 1
62 -999.0 4.050 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Wb `D\r`
-168.2384 -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.23838 111.19493 -1.513 0.1356
Wb -0.10521 0.06038 -1.743 0.0866 .
`D\r` -11.37149 37.66919 -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.5777096 0.2888548
[2,] 0.8966051 0.2067898 0.1033949
[3,] 0.8284961 0.3430078 0.1715039
[4,] 0.8305630 0.3388739 0.1694370
[5,] 0.7969654 0.4060693 0.2030346
[6,] 0.7142142 0.5715716 0.2857858
[7,] 0.6601651 0.6796698 0.3398349
[8,] 0.5869426 0.8261148 0.4130574
[9,] 0.8131895 0.3736210 0.1868105
[10,] 0.7595781 0.4808439 0.2404219
[11,] 0.6980832 0.6038336 0.3019168
[12,] 0.6299787 0.7400426 0.3700213
[13,] 0.5582597 0.8834806 0.4417403
[14,] 0.4837042 0.9674085 0.5162958
[15,] 0.4172345 0.8344691 0.5827655
[16,] 0.4916750 0.9833500 0.5083250
[17,] 0.4449935 0.8899870 0.5550065
[18,] 0.3784247 0.7568494 0.6215753
[19,] 0.5646508 0.8706983 0.4353492
[20,] 0.5035848 0.9928303 0.4964152
[21,] 0.6858010 0.6283980 0.3141990
[22,] 0.6314503 0.7370995 0.3685497
[23,] 0.5746947 0.8506105 0.4253053
[24,] 0.6264877 0.7470246 0.3735123
[25,] 0.7463667 0.5072667 0.2536333
[26,] 0.8698277 0.2603446 0.1301723
[27,] 0.8340923 0.3318155 0.1659077
[28,] 0.7934127 0.4131746 0.2065873
[29,] 0.7545137 0.4909727 0.2454863
[30,] 0.7022357 0.5955286 0.2977643
[31,] 0.6443463 0.7113075 0.3556537
[32,] 0.5828388 0.8343224 0.4171612
[33,] 0.5187759 0.9624482 0.4812241
[34,] 0.4555854 0.9111708 0.5444146
[35,] 0.3951055 0.7902111 0.6048945
[36,] 0.5097253 0.9805494 0.4902747
[37,] 0.4484205 0.8968409 0.5515795
[38,] 0.3856134 0.7712269 0.6143866
[39,] 0.3277583 0.6555166 0.6722417
[40,] 0.2623895 0.5247789 0.7376105
[41,] 0.2112498 0.4224996 0.7887502
[42,] 0.3748013 0.7496025 0.6251987
[43,] 0.3122737 0.6245474 0.6877263
[44,] 0.2435971 0.4871941 0.7564029
[45,] 0.1875761 0.3751521 0.8124239
[46,] 0.1417387 0.2834775 0.8582613
[47,] 0.1058241 0.2116483 0.8941759
[48,] 0.3116498 0.6232996 0.6883502
[49,] 0.2540510 0.5081020 0.7459490
[50,] 0.6399597 0.7200806 0.3600403
[51,] 0.4716912 0.9433823 0.5283088
> postscript(file="/var/www/rcomp/tmp/1h5li1292934786.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/2h5li1292934786.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/3se231292934786.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/4se231292934786.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/5se231292934786.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.56478 208.75807 -819.03399 -796.55035 483.79997 223.93433 195.41229
8 9 10 11 12 13 14
235.75831 190.85707 193.39774 224.76906 277.21954 198.63923 -754.21895
15 16 17 18 19 20 21
185.91776 199.89700 197.66395 200.50240 184.55822 197.92261 -718.24688
22 23 24 25 26 27 28
231.30601 201.99399 -797.61119 193.25291 -815.56777 212.76348 221.23377
29 30 31 32 33 34 35
282.01142 -808.86890 -781.59323 221.42487 197.51092 192.23303 187.82271
36 37 38 39 40 41 42
210.89489 214.25527 213.15790 193.58873 194.27811 -747.60110 206.23186
43 44 45 46 47 48 49
224.77647 203.05181 240.42510 232.85887 -807.56749 212.98231 187.45544
50 51 52 53 54 55 56
200.09681 196.76027 207.63162 -772.34387 234.13512 -807.87134 199.08768
57 58 59 60 61 62
202.07605 207.46328 204.19230 223.86518 192.77811 -818.96402
> postscript(file="/var/www/rcomp/tmp/6k5j61292934786.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.56478 NA
1 208.75807 -96.56478
2 -819.03399 208.75807
3 -796.55035 -819.03399
4 483.79997 -796.55035
5 223.93433 483.79997
6 195.41229 223.93433
7 235.75831 195.41229
8 190.85707 235.75831
9 193.39774 190.85707
10 224.76906 193.39774
11 277.21954 224.76906
12 198.63923 277.21954
13 -754.21895 198.63923
14 185.91776 -754.21895
15 199.89700 185.91776
16 197.66395 199.89700
17 200.50240 197.66395
18 184.55822 200.50240
19 197.92261 184.55822
20 -718.24688 197.92261
21 231.30601 -718.24688
22 201.99399 231.30601
23 -797.61119 201.99399
24 193.25291 -797.61119
25 -815.56777 193.25291
26 212.76348 -815.56777
27 221.23377 212.76348
28 282.01142 221.23377
29 -808.86890 282.01142
30 -781.59323 -808.86890
31 221.42487 -781.59323
32 197.51092 221.42487
33 192.23303 197.51092
34 187.82271 192.23303
35 210.89489 187.82271
36 214.25527 210.89489
37 213.15790 214.25527
38 193.58873 213.15790
39 194.27811 193.58873
40 -747.60110 194.27811
41 206.23186 -747.60110
42 224.77647 206.23186
43 203.05181 224.77647
44 240.42510 203.05181
45 232.85887 240.42510
46 -807.56749 232.85887
47 212.98231 -807.56749
48 187.45544 212.98231
49 200.09681 187.45544
50 196.76027 200.09681
51 207.63162 196.76027
52 -772.34387 207.63162
53 234.13512 -772.34387
54 -807.87134 234.13512
55 199.08768 -807.87134
56 202.07605 199.08768
57 207.46328 202.07605
58 204.19230 207.46328
59 223.86518 204.19230
60 192.77811 223.86518
61 -818.96402 192.77811
62 NA -818.96402
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 208.7581 -96.56478
[2,] -819.0340 208.75807
[3,] -796.5504 -819.03399
[4,] 483.8000 -796.55035
[5,] 223.9343 483.79997
[6,] 195.4123 223.93433
[7,] 235.7583 195.41229
[8,] 190.8571 235.75831
[9,] 193.3977 190.85707
[10,] 224.7691 193.39774
[11,] 277.2195 224.76906
[12,] 198.6392 277.21954
[13,] -754.2190 198.63923
[14,] 185.9178 -754.21895
[15,] 199.8970 185.91776
[16,] 197.6640 199.89700
[17,] 200.5024 197.66395
[18,] 184.5582 200.50240
[19,] 197.9226 184.55822
[20,] -718.2469 197.92261
[21,] 231.3060 -718.24688
[22,] 201.9940 231.30601
[23,] -797.6112 201.99399
[24,] 193.2529 -797.61119
[25,] -815.5678 193.25291
[26,] 212.7635 -815.56777
[27,] 221.2338 212.76348
[28,] 282.0114 221.23377
[29,] -808.8689 282.01142
[30,] -781.5932 -808.86890
[31,] 221.4249 -781.59323
[32,] 197.5109 221.42487
[33,] 192.2330 197.51092
[34,] 187.8227 192.23303
[35,] 210.8949 187.82271
[36,] 214.2553 210.89489
[37,] 213.1579 214.25527
[38,] 193.5887 213.15790
[39,] 194.2781 193.58873
[40,] -747.6011 194.27811
[41,] 206.2319 -747.60110
[42,] 224.7765 206.23186
[43,] 203.0518 224.77647
[44,] 240.4251 203.05181
[45,] 232.8589 240.42510
[46,] -807.5675 232.85887
[47,] 212.9823 -807.56749
[48,] 187.4554 212.98231
[49,] 200.0968 187.45544
[50,] 196.7603 200.09681
[51,] 207.6316 196.76027
[52,] -772.3439 207.63162
[53,] 234.1351 -772.34387
[54,] -807.8713 234.13512
[55,] 199.0877 -807.87134
[56,] 202.0761 199.08768
[57,] 207.4633 202.07605
[58,] 204.1923 207.46328
[59,] 223.8652 204.19230
[60,] 192.7781 223.86518
[61,] -818.9640 192.77811
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 208.7581 -96.56478
2 -819.0340 208.75807
3 -796.5504 -819.03399
4 483.8000 -796.55035
5 223.9343 483.79997
6 195.4123 223.93433
7 235.7583 195.41229
8 190.8571 235.75831
9 193.3977 190.85707
10 224.7691 193.39774
11 277.2195 224.76906
12 198.6392 277.21954
13 -754.2190 198.63923
14 185.9178 -754.21895
15 199.8970 185.91776
16 197.6640 199.89700
17 200.5024 197.66395
18 184.5582 200.50240
19 197.9226 184.55822
20 -718.2469 197.92261
21 231.3060 -718.24688
22 201.9940 231.30601
23 -797.6112 201.99399
24 193.2529 -797.61119
25 -815.5678 193.25291
26 212.7635 -815.56777
27 221.2338 212.76348
28 282.0114 221.23377
29 -808.8689 282.01142
30 -781.5932 -808.86890
31 221.4249 -781.59323
32 197.5109 221.42487
33 192.2330 197.51092
34 187.8227 192.23303
35 210.8949 187.82271
36 214.2553 210.89489
37 213.1579 214.25527
38 193.5887 213.15790
39 194.2781 193.58873
40 -747.6011 194.27811
41 206.2319 -747.60110
42 224.7765 206.23186
43 203.0518 224.77647
44 240.4251 203.05181
45 232.8589 240.42510
46 -807.5675 232.85887
47 212.9823 -807.56749
48 187.4554 212.98231
49 200.0968 187.45544
50 196.7603 200.09681
51 207.6316 196.76027
52 -772.3439 207.63162
53 234.1351 -772.34387
54 -807.8713 234.13512
55 199.0877 -807.87134
56 202.0761 199.08768
57 207.4633 202.07605
58 204.1923 207.46328
59 223.8652 204.19230
60 192.7781 223.86518
61 -818.9640 192.77811
> 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/7dxi91292934786.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/8dxi91292934786.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/9dxi91292934786.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/1066hu1292934786.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/1197y01292934786.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/12iamc1292934786.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/139hcx1292934786.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/14uzt21292934786.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/1559an1292934786.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/16j0qw1292934786.tab")
+ }
>
> try(system("convert tmp/1h5li1292934786.ps tmp/1h5li1292934786.png",intern=TRUE))
character(0)
> try(system("convert tmp/2h5li1292934786.ps tmp/2h5li1292934786.png",intern=TRUE))
character(0)
> try(system("convert tmp/3se231292934786.ps tmp/3se231292934786.png",intern=TRUE))
character(0)
> try(system("convert tmp/4se231292934786.ps tmp/4se231292934786.png",intern=TRUE))
character(0)
> try(system("convert tmp/5se231292934786.ps tmp/5se231292934786.png",intern=TRUE))
character(0)
> try(system("convert tmp/6k5j61292934786.ps tmp/6k5j61292934786.png",intern=TRUE))
character(0)
> try(system("convert tmp/7dxi91292934786.ps tmp/7dxi91292934786.png",intern=TRUE))
character(0)
> try(system("convert tmp/8dxi91292934786.ps tmp/8dxi91292934786.png",intern=TRUE))
character(0)
> try(system("convert tmp/9dxi91292934786.ps tmp/9dxi91292934786.png",intern=TRUE))
character(0)
> try(system("convert tmp/1066hu1292934786.ps tmp/1066hu1292934786.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.170 0.840 3.986