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
+ ,3.823082797
+ ,3.00
+ ,6.30
+ ,0
+ ,3.00
+ ,-999.00
+ ,0.529558673
+ ,1.00
+ ,-999.00
+ ,-0.036212173
+ ,3.00
+ ,2.10
+ ,3.406028945
+ ,4.00
+ ,9.10
+ ,1.02325246
+ ,4.00
+ ,15.80
+ ,-1.638272164
+ ,1.00
+ ,5.20
+ ,2.204119983
+ ,4.00
+ ,10.90
+ ,0.51851394
+ ,1.00
+ ,8.30
+ ,1.717337583
+ ,1.00
+ ,11.00
+ ,-0.37161107
+ ,4.00
+ ,3.20
+ ,2.667452953
+ ,5.00
+ ,7.60
+ ,-0.259637311
+ ,2.00
+ ,-999.00
+ ,2.272073788
+ ,5.00
+ ,6.30
+ ,-1.124938737
+ ,1.00
+ ,8.60
+ ,0.477121255
+ ,2.00
+ ,6.60
+ ,-0.105130343
+ ,2.00
+ ,9.50
+ ,-0.698970004
+ ,2.00
+ ,4.80
+ ,0.149219113
+ ,1.00
+ ,12.00
+ ,1.77815125
+ ,1.00
+ ,-999.00
+ ,2.723455672
+ ,5.00
+ ,3.30
+ ,1.441852176
+ ,5.00
+ ,11.00
+ ,-0.920818754
+ ,2.00
+ ,-999.00
+ ,2.315970345
+ ,1.00
+ ,4.70
+ ,1.929418926
+ ,1.00
+ ,-999.00
+ ,1.560265398
+ ,1.00
+ ,10.40
+ ,-0.995678626
+ ,3.00
+ ,7.40
+ ,0.017033339
+ ,4.00
+ ,2.10
+ ,2.716837723
+ ,5.00
+ ,-999.00
+ ,2
+ ,1.00
+ ,-999.00
+ ,1.544068044
+ ,4.00
+ ,7.70
+ ,-2.301029996
+ ,4.00
+ ,17.90
+ ,-2
+ ,1.00
+ ,6.10
+ ,1.792391689
+ ,1.00
+ ,8.20
+ ,-0.913640169
+ ,1.00
+ ,8.40
+ ,0.130333768
+ ,3.00
+ ,11.90
+ ,-1.638272164
+ ,3.00
+ ,10.80
+ ,-1.318758763
+ ,3.00
+ ,13.80
+ ,0.230448921
+ ,1.00
+ ,14.30
+ ,0.544068044
+ ,1.00
+ ,-999.00
+ ,2.397940009
+ ,5.00
+ ,15.20
+ ,-0.318758763
+ ,2.00
+ ,10.00
+ ,1
+ ,4.00
+ ,11.90
+ ,0.209515015
+ ,2.00
+ ,6.50
+ ,2.283301229
+ ,4.00
+ ,7.50
+ ,0.397940009
+ ,5.00
+ ,-999.00
+ ,0.632254777
+ ,2.00
+ ,10.60
+ ,-0.552841969
+ ,3.00
+ ,7.40
+ ,0.626853415
+ ,1.00
+ ,8.40
+ ,0.832508913
+ ,2.00
+ ,5.70
+ ,-0.124938737
+ ,2.00
+ ,4.90
+ ,0.556302501
+ ,3.00
+ ,-999.00
+ ,1.171141151
+ ,5.00
+ ,3.20
+ ,1.744292983
+ ,5.00
+ ,-999.00
+ ,0.146128036
+ ,2.00
+ ,8.10
+ ,-1.22184875
+ ,2.00
+ ,11.00
+ ,-0.045757491
+ ,2.00
+ ,4.90
+ ,0.301029996
+ ,3.00
+ ,13.20
+ ,-0.982966661
+ ,2.00
+ ,9.70
+ ,0.622214023
+ ,4.00
+ ,12.80
+ ,0.544068044
+ ,1.00
+ ,-999.00
+ ,0.607455023
+ ,1.00)
+ ,dim=c(3
+ ,62)
+ ,dimnames=list(c('SWS'
+ ,'logWb'
+ ,'D')
+ ,1:62))
> y <- array(NA,dim=c(3,62),dimnames=list(c('SWS','logWb','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 logWb D
1 -999.0 3.82308280 3
2 6.3 0.00000000 3
3 -999.0 0.52955867 1
4 -999.0 -0.03621217 3
5 2.1 3.40602895 4
6 9.1 1.02325246 4
7 15.8 -1.63827216 1
8 5.2 2.20411998 4
9 10.9 0.51851394 1
10 8.3 1.71733758 1
11 11.0 -0.37161107 4
12 3.2 2.66745295 5
13 7.6 -0.25963731 2
14 -999.0 2.27207379 5
15 6.3 -1.12493874 1
16 8.6 0.47712126 2
17 6.6 -0.10513034 2
18 9.5 -0.69897000 2
19 4.8 0.14921911 1
20 12.0 1.77815125 1
21 -999.0 2.72345567 5
22 3.3 1.44185218 5
23 11.0 -0.92081875 2
24 -999.0 2.31597035 1
25 4.7 1.92941893 1
26 -999.0 1.56026540 1
27 10.4 -0.99567863 3
28 7.4 0.01703334 4
29 2.1 2.71683772 5
30 -999.0 2.00000000 1
31 -999.0 1.54406804 4
32 7.7 -2.30103000 4
33 17.9 -2.00000000 1
34 6.1 1.79239169 1
35 8.2 -0.91364017 1
36 8.4 0.13033377 3
37 11.9 -1.63827216 3
38 10.8 -1.31875876 3
39 13.8 0.23044892 1
40 14.3 0.54406804 1
41 -999.0 2.39794001 5
42 15.2 -0.31875876 2
43 10.0 1.00000000 4
44 11.9 0.20951501 2
45 6.5 2.28330123 4
46 7.5 0.39794001 5
47 -999.0 0.63225478 2
48 10.6 -0.55284197 3
49 7.4 0.62685342 1
50 8.4 0.83250891 2
51 5.7 -0.12493874 2
52 4.9 0.55630250 3
53 -999.0 1.17114115 5
54 3.2 1.74429298 5
55 -999.0 0.14612804 2
56 8.1 -1.22184875 2
57 11.0 -0.04575749 2
58 4.9 0.30103000 3
59 13.2 -0.98296666 2
60 9.7 0.62221402 4
61 12.8 0.54406804 1
62 -999.0 0.60745502 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) logWb D
-193.77 -129.45 19.18
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-867.44 -65.79 135.98 257.50 560.09
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -193.77 105.22 -1.842 0.07056 .
logWb -129.45 39.53 -3.275 0.00177 **
D 19.18 37.20 0.515 0.60819
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 396.4 on 59 degrees of freedom
Multiple R-squared: 0.1577, Adjusted R-squared: 0.1292
F-statistic: 5.525 on 2 and 59 DF, p-value: 0.006319
> 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.8421806 0.3156388 0.15781940
[2,] 0.9163990 0.1672019 0.08360095
[3,] 0.8784153 0.2431694 0.12158472
[4,] 0.9223707 0.1552586 0.07762928
[5,] 0.9333575 0.1332851 0.06664254
[6,] 0.8931270 0.2137460 0.10687301
[7,] 0.8660992 0.2678015 0.13390076
[8,] 0.8193539 0.3612923 0.18064614
[9,] 0.9091591 0.1816817 0.09084086
[10,] 0.8698024 0.2603953 0.13019764
[11,] 0.8333493 0.3333014 0.16665069
[12,] 0.7819216 0.4361567 0.21807835
[13,] 0.7166461 0.5667077 0.28335386
[14,] 0.6563191 0.6873619 0.34368095
[15,] 0.6382817 0.7234367 0.36171833
[16,] 0.6963499 0.6073002 0.30365011
[17,] 0.6687371 0.6625258 0.33126289
[18,] 0.5942751 0.8114499 0.40572493
[19,] 0.6567928 0.6864143 0.34320717
[20,] 0.6604690 0.6790621 0.33953103
[21,] 0.7462977 0.5074045 0.25370225
[22,] 0.6807070 0.6385859 0.31929296
[23,] 0.6154689 0.7690621 0.38453107
[24,] 0.6354565 0.7290871 0.36454354
[25,] 0.6977328 0.6045344 0.30226718
[26,] 0.8081160 0.3837679 0.19188396
[27,] 0.7598789 0.4802421 0.24012106
[28,] 0.6981814 0.6036373 0.30181864
[29,] 0.6951138 0.6097723 0.30488617
[30,] 0.6256548 0.7486904 0.37434522
[31,] 0.5617087 0.8765827 0.43829134
[32,] 0.4871195 0.9742390 0.51288052
[33,] 0.4109560 0.8219121 0.58904395
[34,] 0.3562099 0.7124198 0.64379008
[35,] 0.3173392 0.6346784 0.68266079
[36,] 0.3985240 0.7970480 0.60147601
[37,] 0.3313116 0.6626232 0.66868840
[38,] 0.2801495 0.5602991 0.71985046
[39,] 0.2309174 0.4618349 0.76908256
[40,] 0.2355621 0.4711242 0.76443790
[41,] 0.1793890 0.3587780 0.82061099
[42,] 0.3127652 0.6255303 0.68723483
[43,] 0.2375453 0.4750907 0.76245466
[44,] 0.2012904 0.4025809 0.79870957
[45,] 0.1812748 0.3625496 0.81872519
[46,] 0.1329031 0.2658063 0.86709685
[47,] 0.1029086 0.2058172 0.89709141
[48,] 0.3244935 0.6489871 0.67550646
[49,] 0.2276096 0.4552191 0.77239045
[50,] 0.5159504 0.9680992 0.48404961
[51,] 0.3582489 0.7164978 0.64175109
> postscript(file="/var/www/html/freestat/rcomp/tmp/1ogrp1293002492.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/2ogrp1293002492.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/3z7qa1293002492.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/4z7qa1293002492.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/5z7qa1293002492.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
-367.845038 142.545974 -755.849856 -867.441795 560.090506 258.633284
7 8 9 10 11 12
-21.681750 407.599958 252.620372 405.211517 79.964216 446.404091
13 14 15 16 17 18
129.410826 -606.978872 35.270728 225.786327 148.412194 74.437954
19 20 21 22 23 24
198.714102 416.784020 -548.546197 287.846569 47.219000 -524.593762
25 26 27 28 29 30
429.066052 -622.421926 17.752532 126.675345 451.697090 -565.497026
31 32 33 34 35 36
-682.045651 -173.104575 -66.408452 412.727485 64.523932 161.518052
37 38 39 40 41 42
-63.933037 -23.671115 218.229532 259.328424 -590.685130 129.357385
43 44 45 46 47 48
256.523187 194.443935 419.150196 156.909157 -761.731196 75.279002
49 50 51 52 53 54
263.145226 271.592274 144.947941 213.160922 -749.497746 326.898396
55 56 57 58 59 60
-824.661691 5.349807 160.498180 180.115167 41.373776 207.317713
61 62
257.828424 -745.765951
> postscript(file="/var/www/html/freestat/rcomp/tmp/6sg8d1293002492.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 -367.845038 NA
1 142.545974 -367.845038
2 -755.849856 142.545974
3 -867.441795 -755.849856
4 560.090506 -867.441795
5 258.633284 560.090506
6 -21.681750 258.633284
7 407.599958 -21.681750
8 252.620372 407.599958
9 405.211517 252.620372
10 79.964216 405.211517
11 446.404091 79.964216
12 129.410826 446.404091
13 -606.978872 129.410826
14 35.270728 -606.978872
15 225.786327 35.270728
16 148.412194 225.786327
17 74.437954 148.412194
18 198.714102 74.437954
19 416.784020 198.714102
20 -548.546197 416.784020
21 287.846569 -548.546197
22 47.219000 287.846569
23 -524.593762 47.219000
24 429.066052 -524.593762
25 -622.421926 429.066052
26 17.752532 -622.421926
27 126.675345 17.752532
28 451.697090 126.675345
29 -565.497026 451.697090
30 -682.045651 -565.497026
31 -173.104575 -682.045651
32 -66.408452 -173.104575
33 412.727485 -66.408452
34 64.523932 412.727485
35 161.518052 64.523932
36 -63.933037 161.518052
37 -23.671115 -63.933037
38 218.229532 -23.671115
39 259.328424 218.229532
40 -590.685130 259.328424
41 129.357385 -590.685130
42 256.523187 129.357385
43 194.443935 256.523187
44 419.150196 194.443935
45 156.909157 419.150196
46 -761.731196 156.909157
47 75.279002 -761.731196
48 263.145226 75.279002
49 271.592274 263.145226
50 144.947941 271.592274
51 213.160922 144.947941
52 -749.497746 213.160922
53 326.898396 -749.497746
54 -824.661691 326.898396
55 5.349807 -824.661691
56 160.498180 5.349807
57 180.115167 160.498180
58 41.373776 180.115167
59 207.317713 41.373776
60 257.828424 207.317713
61 -745.765951 257.828424
62 NA -745.765951
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 142.545974 -367.845038
[2,] -755.849856 142.545974
[3,] -867.441795 -755.849856
[4,] 560.090506 -867.441795
[5,] 258.633284 560.090506
[6,] -21.681750 258.633284
[7,] 407.599958 -21.681750
[8,] 252.620372 407.599958
[9,] 405.211517 252.620372
[10,] 79.964216 405.211517
[11,] 446.404091 79.964216
[12,] 129.410826 446.404091
[13,] -606.978872 129.410826
[14,] 35.270728 -606.978872
[15,] 225.786327 35.270728
[16,] 148.412194 225.786327
[17,] 74.437954 148.412194
[18,] 198.714102 74.437954
[19,] 416.784020 198.714102
[20,] -548.546197 416.784020
[21,] 287.846569 -548.546197
[22,] 47.219000 287.846569
[23,] -524.593762 47.219000
[24,] 429.066052 -524.593762
[25,] -622.421926 429.066052
[26,] 17.752532 -622.421926
[27,] 126.675345 17.752532
[28,] 451.697090 126.675345
[29,] -565.497026 451.697090
[30,] -682.045651 -565.497026
[31,] -173.104575 -682.045651
[32,] -66.408452 -173.104575
[33,] 412.727485 -66.408452
[34,] 64.523932 412.727485
[35,] 161.518052 64.523932
[36,] -63.933037 161.518052
[37,] -23.671115 -63.933037
[38,] 218.229532 -23.671115
[39,] 259.328424 218.229532
[40,] -590.685130 259.328424
[41,] 129.357385 -590.685130
[42,] 256.523187 129.357385
[43,] 194.443935 256.523187
[44,] 419.150196 194.443935
[45,] 156.909157 419.150196
[46,] -761.731196 156.909157
[47,] 75.279002 -761.731196
[48,] 263.145226 75.279002
[49,] 271.592274 263.145226
[50,] 144.947941 271.592274
[51,] 213.160922 144.947941
[52,] -749.497746 213.160922
[53,] 326.898396 -749.497746
[54,] -824.661691 326.898396
[55,] 5.349807 -824.661691
[56,] 160.498180 5.349807
[57,] 180.115167 160.498180
[58,] 41.373776 180.115167
[59,] 207.317713 41.373776
[60,] 257.828424 207.317713
[61,] -745.765951 257.828424
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 142.545974 -367.845038
2 -755.849856 142.545974
3 -867.441795 -755.849856
4 560.090506 -867.441795
5 258.633284 560.090506
6 -21.681750 258.633284
7 407.599958 -21.681750
8 252.620372 407.599958
9 405.211517 252.620372
10 79.964216 405.211517
11 446.404091 79.964216
12 129.410826 446.404091
13 -606.978872 129.410826
14 35.270728 -606.978872
15 225.786327 35.270728
16 148.412194 225.786327
17 74.437954 148.412194
18 198.714102 74.437954
19 416.784020 198.714102
20 -548.546197 416.784020
21 287.846569 -548.546197
22 47.219000 287.846569
23 -524.593762 47.219000
24 429.066052 -524.593762
25 -622.421926 429.066052
26 17.752532 -622.421926
27 126.675345 17.752532
28 451.697090 126.675345
29 -565.497026 451.697090
30 -682.045651 -565.497026
31 -173.104575 -682.045651
32 -66.408452 -173.104575
33 412.727485 -66.408452
34 64.523932 412.727485
35 161.518052 64.523932
36 -63.933037 161.518052
37 -23.671115 -63.933037
38 218.229532 -23.671115
39 259.328424 218.229532
40 -590.685130 259.328424
41 129.357385 -590.685130
42 256.523187 129.357385
43 194.443935 256.523187
44 419.150196 194.443935
45 156.909157 419.150196
46 -761.731196 156.909157
47 75.279002 -761.731196
48 263.145226 75.279002
49 271.592274 263.145226
50 144.947941 271.592274
51 213.160922 144.947941
52 -749.497746 213.160922
53 326.898396 -749.497746
54 -824.661691 326.898396
55 5.349807 -824.661691
56 160.498180 5.349807
57 180.115167 160.498180
58 41.373776 180.115167
59 207.317713 41.373776
60 257.828424 207.317713
61 -745.765951 257.828424
> 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/7kppy1293002492.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/8dh6j1293002492.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/9dh6j1293002492.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/10dh6j1293002492.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/11zz561293002492.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/12ki3u1293002492.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/13gsj31293002492.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/14r1i61293002492.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/15cjzu1293002492.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/168bx31293002492.tab")
+ }
>
> try(system("convert tmp/1ogrp1293002492.ps tmp/1ogrp1293002492.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ogrp1293002492.ps tmp/2ogrp1293002492.png",intern=TRUE))
character(0)
> try(system("convert tmp/3z7qa1293002492.ps tmp/3z7qa1293002492.png",intern=TRUE))
character(0)
> try(system("convert tmp/4z7qa1293002492.ps tmp/4z7qa1293002492.png",intern=TRUE))
character(0)
> try(system("convert tmp/5z7qa1293002492.ps tmp/5z7qa1293002492.png",intern=TRUE))
character(0)
> try(system("convert tmp/6sg8d1293002492.ps tmp/6sg8d1293002492.png",intern=TRUE))
character(0)
> try(system("convert tmp/7kppy1293002492.ps tmp/7kppy1293002492.png",intern=TRUE))
character(0)
> try(system("convert tmp/8dh6j1293002492.ps tmp/8dh6j1293002492.png",intern=TRUE))
character(0)
> try(system("convert tmp/9dh6j1293002492.ps tmp/9dh6j1293002492.png",intern=TRUE))
character(0)
> try(system("convert tmp/10dh6j1293002492.ps tmp/10dh6j1293002492.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.920 2.503 4.285