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(785.8,35,819.3,31.3,849.4,30,880.4,31.3,900.1,33,937.2,31.3,948.9,29,952.6,28.7,947.3,28,974.2,29.7,1000.8,30.7,1032.8,24,1050.7,29,1057.3,33,1075.4,28,1118.4,28.7,1179.8,31.7,1227,34,1257.8,35.3,1251.5,27,1236.3,31.3,1170.6,38.7,1213.1,37.3,1265.5,37.3,1300.8,37.7,1348.4,34.7,1371.9,34.7,1403.3,33.7,1451.8,38.3,1474.2,38,1438.2,38.3,1513.6,42.7,1562.2,41.7,1546.2,39.7,1527.5,39.3,1418.7,39.3,1448.5,37.7,1492.1,38.3,1395.4,37.7,1403.7,37,1316.6,34.3,1274.5,29.7,1264.4,34.7,1323.9,32,1332.1,30.3,1250.2,28.3,1096.7,31.3,1080.8,17.7,1039.2,15.7,792,14.3,746.6,13.3,688.8,11,715.8,2.7,672.9,3.3,629.5,3.7,681.2,1.4,755.4,7.1,760.6,8.1,765.9,12.4,836.8,12.4,904.9,9.2),dim=c(2,61),dimnames=list(c('Herdiv','handact'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Herdiv','handact'),1:61))
> 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 = 'Include Monthly 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
Herdiv handact M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 785.8 35.0 1 0 0 0 0 0 0 0 0 0 0
2 819.3 31.3 0 1 0 0 0 0 0 0 0 0 0
3 849.4 30.0 0 0 1 0 0 0 0 0 0 0 0
4 880.4 31.3 0 0 0 1 0 0 0 0 0 0 0
5 900.1 33.0 0 0 0 0 1 0 0 0 0 0 0
6 937.2 31.3 0 0 0 0 0 1 0 0 0 0 0
7 948.9 29.0 0 0 0 0 0 0 1 0 0 0 0
8 952.6 28.7 0 0 0 0 0 0 0 1 0 0 0
9 947.3 28.0 0 0 0 0 0 0 0 0 1 0 0
10 974.2 29.7 0 0 0 0 0 0 0 0 0 1 0
11 1000.8 30.7 0 0 0 0 0 0 0 0 0 0 1
12 1032.8 24.0 0 0 0 0 0 0 0 0 0 0 0
13 1050.7 29.0 1 0 0 0 0 0 0 0 0 0 0
14 1057.3 33.0 0 1 0 0 0 0 0 0 0 0 0
15 1075.4 28.0 0 0 1 0 0 0 0 0 0 0 0
16 1118.4 28.7 0 0 0 1 0 0 0 0 0 0 0
17 1179.8 31.7 0 0 0 0 1 0 0 0 0 0 0
18 1227.0 34.0 0 0 0 0 0 1 0 0 0 0 0
19 1257.8 35.3 0 0 0 0 0 0 1 0 0 0 0
20 1251.5 27.0 0 0 0 0 0 0 0 1 0 0 0
21 1236.3 31.3 0 0 0 0 0 0 0 0 1 0 0
22 1170.6 38.7 0 0 0 0 0 0 0 0 0 1 0
23 1213.1 37.3 0 0 0 0 0 0 0 0 0 0 1
24 1265.5 37.3 0 0 0 0 0 0 0 0 0 0 0
25 1300.8 37.7 1 0 0 0 0 0 0 0 0 0 0
26 1348.4 34.7 0 1 0 0 0 0 0 0 0 0 0
27 1371.9 34.7 0 0 1 0 0 0 0 0 0 0 0
28 1403.3 33.7 0 0 0 1 0 0 0 0 0 0 0
29 1451.8 38.3 0 0 0 0 1 0 0 0 0 0 0
30 1474.2 38.0 0 0 0 0 0 1 0 0 0 0 0
31 1438.2 38.3 0 0 0 0 0 0 1 0 0 0 0
32 1513.6 42.7 0 0 0 0 0 0 0 1 0 0 0
33 1562.2 41.7 0 0 0 0 0 0 0 0 1 0 0
34 1546.2 39.7 0 0 0 0 0 0 0 0 0 1 0
35 1527.5 39.3 0 0 0 0 0 0 0 0 0 0 1
36 1418.7 39.3 0 0 0 0 0 0 0 0 0 0 0
37 1448.5 37.7 1 0 0 0 0 0 0 0 0 0 0
38 1492.1 38.3 0 1 0 0 0 0 0 0 0 0 0
39 1395.4 37.7 0 0 1 0 0 0 0 0 0 0 0
40 1403.7 37.0 0 0 0 1 0 0 0 0 0 0 0
41 1316.6 34.3 0 0 0 0 1 0 0 0 0 0 0
42 1274.5 29.7 0 0 0 0 0 1 0 0 0 0 0
43 1264.4 34.7 0 0 0 0 0 0 1 0 0 0 0
44 1323.9 32.0 0 0 0 0 0 0 0 1 0 0 0
45 1332.1 30.3 0 0 0 0 0 0 0 0 1 0 0
46 1250.2 28.3 0 0 0 0 0 0 0 0 0 1 0
47 1096.7 31.3 0 0 0 0 0 0 0 0 0 0 1
48 1080.8 17.7 0 0 0 0 0 0 0 0 0 0 0
49 1039.2 15.7 1 0 0 0 0 0 0 0 0 0 0
50 792.0 14.3 0 1 0 0 0 0 0 0 0 0 0
51 746.6 13.3 0 0 1 0 0 0 0 0 0 0 0
52 688.8 11.0 0 0 0 1 0 0 0 0 0 0 0
53 715.8 2.7 0 0 0 0 1 0 0 0 0 0 0
54 672.9 3.3 0 0 0 0 0 1 0 0 0 0 0
55 629.5 3.7 0 0 0 0 0 0 1 0 0 0 0
56 681.2 1.4 0 0 0 0 0 0 0 1 0 0 0
57 755.4 7.1 0 0 0 0 0 0 0 0 1 0 0
58 760.6 8.1 0 0 0 0 0 0 0 0 0 1 0
59 765.9 12.4 0 0 0 0 0 0 0 0 0 0 1
60 836.8 12.4 0 0 0 0 0 0 0 0 0 0 0
61 904.9 9.2 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) handact M1 M2 M3 M4
608.35 19.84 -63.27 -108.02 -90.76 -71.64
M5 M6 M7 M8 M9 M10
-51.00 -31.98 -60.03 13.28 9.19 -41.31
M11
-86.66
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-453.62 -66.13 9.19 121.31 231.97
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 608.354 91.470 6.651 2.52e-08 ***
handact 19.838 1.961 10.117 1.73e-13 ***
M1 -63.269 102.607 -0.617 0.540
M2 -108.023 107.452 -1.005 0.320
M3 -90.759 107.260 -0.846 0.402
M4 -71.644 107.226 -0.668 0.507
M5 -50.999 107.201 -0.476 0.636
M6 -31.979 107.162 -0.298 0.767
M7 -60.026 107.215 -0.560 0.578
M8 13.276 107.140 0.124 0.902
M9 9.189 107.182 0.086 0.932
M10 -41.313 107.276 -0.385 0.702
M11 -86.662 107.434 -0.807 0.424
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 169.4 on 48 degrees of freedom
Multiple R-squared: 0.6831, Adjusted R-squared: 0.6039
F-statistic: 8.622 on 12 and 48 DF, p-value: 2.046e-08
> 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.9178210 1.643581e-01 8.217904e-02
[2,] 0.9377856 1.244288e-01 6.221442e-02
[3,] 0.9920301 1.593979e-02 7.969893e-03
[4,] 0.9977894 4.421190e-03 2.210595e-03
[5,] 0.9981394 3.721165e-03 1.860582e-03
[6,] 0.9983694 3.261252e-03 1.630626e-03
[7,] 0.9996514 6.971157e-04 3.485579e-04
[8,] 0.9995608 8.783408e-04 4.391704e-04
[9,] 0.9995731 8.537618e-04 4.268809e-04
[10,] 0.9999461 1.078346e-04 5.391732e-05
[11,] 0.9999665 6.705687e-05 3.352844e-05
[12,] 0.9999676 6.478024e-05 3.239012e-05
[13,] 0.9999833 3.346137e-05 1.673068e-05
[14,] 0.9999669 6.618951e-05 3.309476e-05
[15,] 0.9999302 1.395710e-04 6.978551e-05
[16,] 0.9998453 3.094298e-04 1.547149e-04
[17,] 0.9997206 5.588906e-04 2.794453e-04
[18,] 0.9992822 1.435695e-03 7.178473e-04
[19,] 0.9989023 2.195370e-03 1.097685e-03
[20,] 0.9996923 6.153239e-04 3.076620e-04
[21,] 0.9994265 1.147094e-03 5.735468e-04
[22,] 0.9993377 1.324536e-03 6.622680e-04
[23,] 0.9995089 9.821299e-04 4.910650e-04
[24,] 0.9991871 1.625893e-03 8.129466e-04
[25,] 0.9994846 1.030894e-03 5.154471e-04
[26,] 0.9990005 1.998925e-03 9.994624e-04
[27,] 0.9969489 6.102279e-03 3.051139e-03
[28,] 0.9906477 1.870461e-02 9.352305e-03
[29,] 0.9735164 5.296723e-02 2.648362e-02
[30,] 0.9317371 1.365259e-01 6.826295e-02
> postscript(file="/var/www/html/rcomp/tmp/1j2ye1258554214.ps",horizontal=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/2b7sc1258554214.ps",horizontal=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/3yxcz1258554214.ps",horizontal=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/4w4oi1258554214.ps",horizontal=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/58w6v1258554214.ps",horizontal=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 = 61
Frequency = 1
1 2 3 4 5 6
-453.616195 -301.961252 -263.335896 -277.240517 -311.910062 -260.105570
7 8 9 10 11 12
-174.730410 -238.380949 -225.708164 -182.030410 -129.919006 -51.666653
13 14 15 16 17 18
-69.688120 -97.685873 2.340129 12.338316 -6.420646 -23.868204
19 20 21 22 23 24
9.190112 94.243672 -2.173605 -164.172522 -48.549888 -82.812219
25 26 27 28 29 30
7.821172 159.689506 165.925446 198.048253 134.648472 143.979746
31 32 33 34 35 36
130.076074 44.886877 117.411066 191.589466 226.174087 30.711756
37 38 39 40 41 42
155.521172 231.972661 129.911409 132.982812 78.800522 108.935250
43 44 45 46 47 48
27.692919 67.453610 113.464407 121.742807 -45.921814 121.312825
49 50 51 52 53 54
182.657445 7.984959 -34.841088 -66.128864 104.881715 31.058778
55 56 57 58 59 60
7.771305 31.796791 -2.993704 32.870659 -1.783379 -17.545709
61
177.304526
> postscript(file="/var/www/html/rcomp/tmp/6zd8j1258554214.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -453.616195 NA
1 -301.961252 -453.616195
2 -263.335896 -301.961252
3 -277.240517 -263.335896
4 -311.910062 -277.240517
5 -260.105570 -311.910062
6 -174.730410 -260.105570
7 -238.380949 -174.730410
8 -225.708164 -238.380949
9 -182.030410 -225.708164
10 -129.919006 -182.030410
11 -51.666653 -129.919006
12 -69.688120 -51.666653
13 -97.685873 -69.688120
14 2.340129 -97.685873
15 12.338316 2.340129
16 -6.420646 12.338316
17 -23.868204 -6.420646
18 9.190112 -23.868204
19 94.243672 9.190112
20 -2.173605 94.243672
21 -164.172522 -2.173605
22 -48.549888 -164.172522
23 -82.812219 -48.549888
24 7.821172 -82.812219
25 159.689506 7.821172
26 165.925446 159.689506
27 198.048253 165.925446
28 134.648472 198.048253
29 143.979746 134.648472
30 130.076074 143.979746
31 44.886877 130.076074
32 117.411066 44.886877
33 191.589466 117.411066
34 226.174087 191.589466
35 30.711756 226.174087
36 155.521172 30.711756
37 231.972661 155.521172
38 129.911409 231.972661
39 132.982812 129.911409
40 78.800522 132.982812
41 108.935250 78.800522
42 27.692919 108.935250
43 67.453610 27.692919
44 113.464407 67.453610
45 121.742807 113.464407
46 -45.921814 121.742807
47 121.312825 -45.921814
48 182.657445 121.312825
49 7.984959 182.657445
50 -34.841088 7.984959
51 -66.128864 -34.841088
52 104.881715 -66.128864
53 31.058778 104.881715
54 7.771305 31.058778
55 31.796791 7.771305
56 -2.993704 31.796791
57 32.870659 -2.993704
58 -1.783379 32.870659
59 -17.545709 -1.783379
60 177.304526 -17.545709
61 NA 177.304526
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -301.961252 -453.616195
[2,] -263.335896 -301.961252
[3,] -277.240517 -263.335896
[4,] -311.910062 -277.240517
[5,] -260.105570 -311.910062
[6,] -174.730410 -260.105570
[7,] -238.380949 -174.730410
[8,] -225.708164 -238.380949
[9,] -182.030410 -225.708164
[10,] -129.919006 -182.030410
[11,] -51.666653 -129.919006
[12,] -69.688120 -51.666653
[13,] -97.685873 -69.688120
[14,] 2.340129 -97.685873
[15,] 12.338316 2.340129
[16,] -6.420646 12.338316
[17,] -23.868204 -6.420646
[18,] 9.190112 -23.868204
[19,] 94.243672 9.190112
[20,] -2.173605 94.243672
[21,] -164.172522 -2.173605
[22,] -48.549888 -164.172522
[23,] -82.812219 -48.549888
[24,] 7.821172 -82.812219
[25,] 159.689506 7.821172
[26,] 165.925446 159.689506
[27,] 198.048253 165.925446
[28,] 134.648472 198.048253
[29,] 143.979746 134.648472
[30,] 130.076074 143.979746
[31,] 44.886877 130.076074
[32,] 117.411066 44.886877
[33,] 191.589466 117.411066
[34,] 226.174087 191.589466
[35,] 30.711756 226.174087
[36,] 155.521172 30.711756
[37,] 231.972661 155.521172
[38,] 129.911409 231.972661
[39,] 132.982812 129.911409
[40,] 78.800522 132.982812
[41,] 108.935250 78.800522
[42,] 27.692919 108.935250
[43,] 67.453610 27.692919
[44,] 113.464407 67.453610
[45,] 121.742807 113.464407
[46,] -45.921814 121.742807
[47,] 121.312825 -45.921814
[48,] 182.657445 121.312825
[49,] 7.984959 182.657445
[50,] -34.841088 7.984959
[51,] -66.128864 -34.841088
[52,] 104.881715 -66.128864
[53,] 31.058778 104.881715
[54,] 7.771305 31.058778
[55,] 31.796791 7.771305
[56,] -2.993704 31.796791
[57,] 32.870659 -2.993704
[58,] -1.783379 32.870659
[59,] -17.545709 -1.783379
[60,] 177.304526 -17.545709
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -301.961252 -453.616195
2 -263.335896 -301.961252
3 -277.240517 -263.335896
4 -311.910062 -277.240517
5 -260.105570 -311.910062
6 -174.730410 -260.105570
7 -238.380949 -174.730410
8 -225.708164 -238.380949
9 -182.030410 -225.708164
10 -129.919006 -182.030410
11 -51.666653 -129.919006
12 -69.688120 -51.666653
13 -97.685873 -69.688120
14 2.340129 -97.685873
15 12.338316 2.340129
16 -6.420646 12.338316
17 -23.868204 -6.420646
18 9.190112 -23.868204
19 94.243672 9.190112
20 -2.173605 94.243672
21 -164.172522 -2.173605
22 -48.549888 -164.172522
23 -82.812219 -48.549888
24 7.821172 -82.812219
25 159.689506 7.821172
26 165.925446 159.689506
27 198.048253 165.925446
28 134.648472 198.048253
29 143.979746 134.648472
30 130.076074 143.979746
31 44.886877 130.076074
32 117.411066 44.886877
33 191.589466 117.411066
34 226.174087 191.589466
35 30.711756 226.174087
36 155.521172 30.711756
37 231.972661 155.521172
38 129.911409 231.972661
39 132.982812 129.911409
40 78.800522 132.982812
41 108.935250 78.800522
42 27.692919 108.935250
43 67.453610 27.692919
44 113.464407 67.453610
45 121.742807 113.464407
46 -45.921814 121.742807
47 121.312825 -45.921814
48 182.657445 121.312825
49 7.984959 182.657445
50 -34.841088 7.984959
51 -66.128864 -34.841088
52 104.881715 -66.128864
53 31.058778 104.881715
54 7.771305 31.058778
55 31.796791 7.771305
56 -2.993704 31.796791
57 32.870659 -2.993704
58 -1.783379 32.870659
59 -17.545709 -1.783379
60 177.304526 -17.545709
> 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/7fv0z1258554214.ps",horizontal=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/8wsuq1258554214.ps",horizontal=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/94gbc1258554214.ps",horizontal=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/10380a1258554214.ps",horizontal=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/11wu2v1258554214.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/124apq1258554215.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/13wq6j1258554215.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/14o6hu1258554215.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/15i44g1258554215.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/16zy121258554215.tab")
+ }
>
> system("convert tmp/1j2ye1258554214.ps tmp/1j2ye1258554214.png")
> system("convert tmp/2b7sc1258554214.ps tmp/2b7sc1258554214.png")
> system("convert tmp/3yxcz1258554214.ps tmp/3yxcz1258554214.png")
> system("convert tmp/4w4oi1258554214.ps tmp/4w4oi1258554214.png")
> system("convert tmp/58w6v1258554214.ps tmp/58w6v1258554214.png")
> system("convert tmp/6zd8j1258554214.ps tmp/6zd8j1258554214.png")
> system("convert tmp/7fv0z1258554214.ps tmp/7fv0z1258554214.png")
> system("convert tmp/8wsuq1258554214.ps tmp/8wsuq1258554214.png")
> system("convert tmp/94gbc1258554214.ps tmp/94gbc1258554214.png")
> system("convert tmp/10380a1258554214.ps tmp/10380a1258554214.png")
>
>
> proc.time()
user system elapsed
2.382 1.566 2.969