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.000
+ ,645.000
+ ,3.000
+ ,2.000
+ ,42.000
+ ,3.000
+ ,-999.000
+ ,60.000
+ ,1.000
+ ,-999.000
+ ,25.000
+ ,3.000
+ ,1.800
+ ,624.000
+ ,4.000
+ ,0.700
+ ,180.000
+ ,4.000
+ ,3.900
+ ,35.000
+ ,1.000
+ ,1.000
+ ,392.000
+ ,4.000
+ ,3.600
+ ,63.000
+ ,1.000
+ ,1.400
+ ,230.000
+ ,1.000
+ ,1.500
+ ,112.000
+ ,4.000
+ ,0.700
+ ,281.000
+ ,5.000
+ ,2.700
+ ,-999.000
+ ,2.000
+ ,-999.000
+ ,365.000
+ ,5.000
+ ,2.100
+ ,42.000
+ ,1.000
+ ,0.000
+ ,28.000
+ ,2.000
+ ,4.100
+ ,42.000
+ ,2.000
+ ,1.200
+ ,120.000
+ ,2.000
+ ,1.300
+ ,-999.000
+ ,1.000
+ ,6.100
+ ,-999.000
+ ,1.000
+ ,0.300
+ ,400.000
+ ,5.000
+ ,0.500
+ ,148.000
+ ,5.000
+ ,3.400
+ ,16.000
+ ,2.000
+ ,-999.000
+ ,252.000
+ ,1.000
+ ,1.500
+ ,310.000
+ ,1.000
+ ,-999.000
+ ,63.000
+ ,1.000
+ ,3.400
+ ,28.000
+ ,3.000
+ ,0.800
+ ,68.000
+ ,4.000
+ ,0.800
+ ,336.000
+ ,5.000
+ ,-999.000
+ ,100.000
+ ,1.000
+ ,-999.000
+ ,33.000
+ ,4.000
+ ,1.400
+ ,21.500
+ ,4.000
+ ,2.000
+ ,50.000
+ ,1.000
+ ,1.900
+ ,267.000
+ ,1.000
+ ,2.400
+ ,30.000
+ ,1.000
+ ,2.800
+ ,45.000
+ ,3.000
+ ,1.300
+ ,19.000
+ ,3.000
+ ,2.000
+ ,30.000
+ ,3.000
+ ,5.600
+ ,12.000
+ ,1.000
+ ,3.100
+ ,120.000
+ ,1.000
+ ,1.000
+ ,440.000
+ ,5.000
+ ,1.800
+ ,140.000
+ ,2.000
+ ,0.900
+ ,170.000
+ ,4.000
+ ,1.800
+ ,17.000
+ ,2.000
+ ,1.900
+ ,115.000
+ ,4.000
+ ,0.900
+ ,31.000
+ ,5.000
+ ,-999.000
+ ,63.000
+ ,2.000
+ ,2.600
+ ,21.000
+ ,3.000
+ ,2.400
+ ,52.000
+ ,1.000
+ ,1.200
+ ,164.000
+ ,2.000
+ ,0.900
+ ,225.000
+ ,2.000
+ ,0.500
+ ,225.000
+ ,3.000
+ ,-999.000
+ ,150.000
+ ,5.000
+ ,0.600
+ ,151.000
+ ,5.000
+ ,-999.000
+ ,90.000
+ ,2.000
+ ,2.200
+ ,-999.000
+ ,2.000
+ ,2.300
+ ,60.000
+ ,2.000
+ ,0.500
+ ,200.000
+ ,3.000
+ ,2.600
+ ,46.000
+ ,2.000
+ ,0.600
+ ,210.000
+ ,4.000
+ ,6.600
+ ,14.000
+ ,1.000
+ ,-999.000
+ ,38.000
+ ,1.000)
+ ,dim=c(3
+ ,62)
+ ,dimnames=list(c('PS'
+ ,'tg'
+ ,'D
')
+ ,1:62))
> y <- array(NA,dim=c(3,62),dimnames=list(c('PS','tg','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
PS tg D\r
1 -999.0 645.0 3
2 2.0 42.0 3
3 -999.0 60.0 1
4 -999.0 25.0 3
5 1.8 624.0 4
6 0.7 180.0 4
7 3.9 35.0 1
8 1.0 392.0 4
9 3.6 63.0 1
10 1.4 230.0 1
11 1.5 112.0 4
12 0.7 281.0 5
13 2.7 -999.0 2
14 -999.0 365.0 5
15 2.1 42.0 1
16 0.0 28.0 2
17 4.1 42.0 2
18 1.2 120.0 2
19 1.3 -999.0 1
20 6.1 -999.0 1
21 0.3 400.0 5
22 0.5 148.0 5
23 3.4 16.0 2
24 -999.0 252.0 1
25 1.5 310.0 1
26 -999.0 63.0 1
27 3.4 28.0 3
28 0.8 68.0 4
29 0.8 336.0 5
30 -999.0 100.0 1
31 -999.0 33.0 4
32 1.4 21.5 4
33 2.0 50.0 1
34 1.9 267.0 1
35 2.4 30.0 1
36 2.8 45.0 3
37 1.3 19.0 3
38 2.0 30.0 3
39 5.6 12.0 1
40 3.1 120.0 1
41 1.0 440.0 5
42 1.8 140.0 2
43 0.9 170.0 4
44 1.8 17.0 2
45 1.9 115.0 4
46 0.9 31.0 5
47 -999.0 63.0 2
48 2.6 21.0 3
49 2.4 52.0 1
50 1.2 164.0 2
51 0.9 225.0 2
52 0.5 225.0 3
53 -999.0 150.0 5
54 0.6 151.0 5
55 -999.0 90.0 2
56 2.2 -999.0 2
57 2.3 60.0 2
58 0.5 200.0 3
59 2.6 46.0 2
60 0.6 210.0 4
61 6.6 14.0 1
62 -999.0 38.0 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) tg `D\r`
-269.5850 -0.2322 35.8887
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-874.04 30.46 172.74 230.80 307.16
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -269.5850 107.8028 -2.501 0.0152 *
tg -0.2322 0.1720 -1.349 0.1824
`D\r` 35.8887 37.7590 0.950 0.3458
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 398.4 on 59 degrees of freedom
Multiple R-squared: 0.03411, Adjusted R-squared: 0.001367
F-statistic: 1.042 on 2 and 59 DF, p-value: 0.3592
> 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.83985214 0.32029572 0.16014786
[2,] 0.95585954 0.08828092 0.04414046
[3,] 0.92727061 0.14545879 0.07272939
[4,] 0.93411603 0.13176794 0.06588397
[5,] 0.92608663 0.14782674 0.07391337
[6,] 0.88730354 0.22539293 0.11269646
[7,] 0.83449920 0.33100160 0.16550080
[8,] 0.76586623 0.46826754 0.23413377
[9,] 0.90088965 0.19822070 0.09911035
[10,] 0.87243039 0.25513923 0.12756961
[11,] 0.83407824 0.33184352 0.16592176
[12,] 0.78970662 0.42058677 0.21029338
[13,] 0.74153142 0.51693715 0.25846858
[14,] 0.67130229 0.65739542 0.32869771
[15,] 0.59565933 0.80868134 0.40434067
[16,] 0.54295169 0.91409661 0.45704831
[17,] 0.47400285 0.94800570 0.52599715
[18,] 0.41265534 0.82531069 0.58734466
[19,] 0.56966970 0.86066060 0.43033030
[20,] 0.53650596 0.92698809 0.46349404
[21,] 0.70055293 0.59889414 0.29944707
[22,] 0.64381009 0.71237983 0.35618991
[23,] 0.57876168 0.84247665 0.42123832
[24,] 0.51311892 0.97376215 0.48688108
[25,] 0.68580921 0.62838159 0.31419079
[26,] 0.88026212 0.23947576 0.11973788
[27,] 0.84252477 0.31495047 0.15747523
[28,] 0.81153255 0.37693491 0.18846745
[29,] 0.78201030 0.43597940 0.21798970
[30,] 0.74031010 0.51937980 0.25968990
[31,] 0.68491711 0.63016579 0.31508289
[32,] 0.62381365 0.75237269 0.37618635
[33,] 0.55950683 0.88098635 0.44049317
[34,] 0.50459014 0.99081972 0.49540986
[35,] 0.45422286 0.90844572 0.54577714
[36,] 0.38813495 0.77626989 0.61186505
[37,] 0.33453861 0.66907721 0.66546139
[38,] 0.27396018 0.54792036 0.72603982
[39,] 0.22506550 0.45013099 0.77493450
[40,] 0.17578513 0.35157027 0.82421487
[41,] 0.13193890 0.26387780 0.86806110
[42,] 0.27275514 0.54551028 0.72724486
[43,] 0.21543133 0.43086266 0.78456867
[44,] 0.16729681 0.33459362 0.83270319
[45,] 0.12895658 0.25791317 0.87104342
[46,] 0.10081681 0.20163362 0.89918319
[47,] 0.07735298 0.15470596 0.92264702
[48,] 0.28430729 0.56861457 0.71569271
[49,] 0.20205438 0.40410876 0.79794562
[50,] 0.43877150 0.87754300 0.56122850
[51,] 0.36896647 0.73793294 0.63103353
> postscript(file="/var/www/html/freestat/rcomp/tmp/159an1292935130.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/259an1292935130.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/3xir81292935130.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/4xir81292935130.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/5xir81292935130.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
-687.340926 173.669458 -751.374377 -831.277181 272.695129 168.518197
7 8 9 10 11 12
245.721741 218.035110 251.922088 288.492015 153.531640 156.077192
13 14 15 16 17 18
-31.415473 -824.121766 245.546828 204.307969 211.658143 226.866252
19 20 21 22 23 24
3.073212 7.873212 183.303667 125.000544 204.922106 -706.800569
25 26 27 28 29 30
307.164435 -750.677912 171.819284 142.616809 168.945731 -742.088167
31 32 33 34 35 36
-865.308625 132.421590 247.304070 297.581759 243.060965 175.165924
37 38 39 40 41 42
167.629887 170.883595 242.082171 264.654938 193.289877 232.109358
43 44 45 46 47 48
166.396645 203.554262 154.628106 98.238380 -786.566597 169.394198
49 50 51 52 53 54
248.168381 237.081084 250.942554 214.653869 -874.035145 125.797010
55 56 57 58 59 60
-780.298405 -31.915473 214.036937 208.849987 211.086764 175.382855
61 62
243.546481 -756.481793
> postscript(file="/var/www/html/freestat/rcomp/tmp/6q9qb1292935130.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 -687.340926 NA
1 173.669458 -687.340926
2 -751.374377 173.669458
3 -831.277181 -751.374377
4 272.695129 -831.277181
5 168.518197 272.695129
6 245.721741 168.518197
7 218.035110 245.721741
8 251.922088 218.035110
9 288.492015 251.922088
10 153.531640 288.492015
11 156.077192 153.531640
12 -31.415473 156.077192
13 -824.121766 -31.415473
14 245.546828 -824.121766
15 204.307969 245.546828
16 211.658143 204.307969
17 226.866252 211.658143
18 3.073212 226.866252
19 7.873212 3.073212
20 183.303667 7.873212
21 125.000544 183.303667
22 204.922106 125.000544
23 -706.800569 204.922106
24 307.164435 -706.800569
25 -750.677912 307.164435
26 171.819284 -750.677912
27 142.616809 171.819284
28 168.945731 142.616809
29 -742.088167 168.945731
30 -865.308625 -742.088167
31 132.421590 -865.308625
32 247.304070 132.421590
33 297.581759 247.304070
34 243.060965 297.581759
35 175.165924 243.060965
36 167.629887 175.165924
37 170.883595 167.629887
38 242.082171 170.883595
39 264.654938 242.082171
40 193.289877 264.654938
41 232.109358 193.289877
42 166.396645 232.109358
43 203.554262 166.396645
44 154.628106 203.554262
45 98.238380 154.628106
46 -786.566597 98.238380
47 169.394198 -786.566597
48 248.168381 169.394198
49 237.081084 248.168381
50 250.942554 237.081084
51 214.653869 250.942554
52 -874.035145 214.653869
53 125.797010 -874.035145
54 -780.298405 125.797010
55 -31.915473 -780.298405
56 214.036937 -31.915473
57 208.849987 214.036937
58 211.086764 208.849987
59 175.382855 211.086764
60 243.546481 175.382855
61 -756.481793 243.546481
62 NA -756.481793
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 173.669458 -687.340926
[2,] -751.374377 173.669458
[3,] -831.277181 -751.374377
[4,] 272.695129 -831.277181
[5,] 168.518197 272.695129
[6,] 245.721741 168.518197
[7,] 218.035110 245.721741
[8,] 251.922088 218.035110
[9,] 288.492015 251.922088
[10,] 153.531640 288.492015
[11,] 156.077192 153.531640
[12,] -31.415473 156.077192
[13,] -824.121766 -31.415473
[14,] 245.546828 -824.121766
[15,] 204.307969 245.546828
[16,] 211.658143 204.307969
[17,] 226.866252 211.658143
[18,] 3.073212 226.866252
[19,] 7.873212 3.073212
[20,] 183.303667 7.873212
[21,] 125.000544 183.303667
[22,] 204.922106 125.000544
[23,] -706.800569 204.922106
[24,] 307.164435 -706.800569
[25,] -750.677912 307.164435
[26,] 171.819284 -750.677912
[27,] 142.616809 171.819284
[28,] 168.945731 142.616809
[29,] -742.088167 168.945731
[30,] -865.308625 -742.088167
[31,] 132.421590 -865.308625
[32,] 247.304070 132.421590
[33,] 297.581759 247.304070
[34,] 243.060965 297.581759
[35,] 175.165924 243.060965
[36,] 167.629887 175.165924
[37,] 170.883595 167.629887
[38,] 242.082171 170.883595
[39,] 264.654938 242.082171
[40,] 193.289877 264.654938
[41,] 232.109358 193.289877
[42,] 166.396645 232.109358
[43,] 203.554262 166.396645
[44,] 154.628106 203.554262
[45,] 98.238380 154.628106
[46,] -786.566597 98.238380
[47,] 169.394198 -786.566597
[48,] 248.168381 169.394198
[49,] 237.081084 248.168381
[50,] 250.942554 237.081084
[51,] 214.653869 250.942554
[52,] -874.035145 214.653869
[53,] 125.797010 -874.035145
[54,] -780.298405 125.797010
[55,] -31.915473 -780.298405
[56,] 214.036937 -31.915473
[57,] 208.849987 214.036937
[58,] 211.086764 208.849987
[59,] 175.382855 211.086764
[60,] 243.546481 175.382855
[61,] -756.481793 243.546481
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 173.669458 -687.340926
2 -751.374377 173.669458
3 -831.277181 -751.374377
4 272.695129 -831.277181
5 168.518197 272.695129
6 245.721741 168.518197
7 218.035110 245.721741
8 251.922088 218.035110
9 288.492015 251.922088
10 153.531640 288.492015
11 156.077192 153.531640
12 -31.415473 156.077192
13 -824.121766 -31.415473
14 245.546828 -824.121766
15 204.307969 245.546828
16 211.658143 204.307969
17 226.866252 211.658143
18 3.073212 226.866252
19 7.873212 3.073212
20 183.303667 7.873212
21 125.000544 183.303667
22 204.922106 125.000544
23 -706.800569 204.922106
24 307.164435 -706.800569
25 -750.677912 307.164435
26 171.819284 -750.677912
27 142.616809 171.819284
28 168.945731 142.616809
29 -742.088167 168.945731
30 -865.308625 -742.088167
31 132.421590 -865.308625
32 247.304070 132.421590
33 297.581759 247.304070
34 243.060965 297.581759
35 175.165924 243.060965
36 167.629887 175.165924
37 170.883595 167.629887
38 242.082171 170.883595
39 264.654938 242.082171
40 193.289877 264.654938
41 232.109358 193.289877
42 166.396645 232.109358
43 203.554262 166.396645
44 154.628106 203.554262
45 98.238380 154.628106
46 -786.566597 98.238380
47 169.394198 -786.566597
48 248.168381 169.394198
49 237.081084 248.168381
50 250.942554 237.081084
51 214.653869 250.942554
52 -874.035145 214.653869
53 125.797010 -874.035145
54 -780.298405 125.797010
55 -31.915473 -780.298405
56 214.036937 -31.915473
57 208.849987 214.036937
58 211.086764 208.849987
59 175.382855 211.086764
60 243.546481 175.382855
61 -756.481793 243.546481
> 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/7j0qw1292935130.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/8j0qw1292935130.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/9j0qw1292935130.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/10bs7z1292935130.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/11xs651292935130.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/12q2n81292935130.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/13xl221292935130.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/14030p1292935130.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/15tu0s1292935130.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/1664xj1292935130.tab")
+ }
>
> try(system("convert tmp/159an1292935130.ps tmp/159an1292935130.png",intern=TRUE))
character(0)
> try(system("convert tmp/259an1292935130.ps tmp/259an1292935130.png",intern=TRUE))
character(0)
> try(system("convert tmp/3xir81292935130.ps tmp/3xir81292935130.png",intern=TRUE))
character(0)
> try(system("convert tmp/4xir81292935130.ps tmp/4xir81292935130.png",intern=TRUE))
character(0)
> try(system("convert tmp/5xir81292935130.ps tmp/5xir81292935130.png",intern=TRUE))
character(0)
> try(system("convert tmp/6q9qb1292935130.ps tmp/6q9qb1292935130.png",intern=TRUE))
character(0)
> try(system("convert tmp/7j0qw1292935130.ps tmp/7j0qw1292935130.png",intern=TRUE))
character(0)
> try(system("convert tmp/8j0qw1292935130.ps tmp/8j0qw1292935130.png",intern=TRUE))
character(0)
> try(system("convert tmp/9j0qw1292935130.ps tmp/9j0qw1292935130.png",intern=TRUE))
character(0)
> try(system("convert tmp/10bs7z1292935130.ps tmp/10bs7z1292935130.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.270 2.694 22.865