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(462,1919,455,1911,461,1870,461,2263,463,1802,462,1863,456,1989,455,2197,456,2409,472,2502,472,2593,471,2598,465,2053,459,2213,465,2238,468,2359,467,2151,463,2474,460,3079,462,2312,461,2565,476,1972,476,2484,471,2202,453,2151,443,1976,442,2012,444,2114,438,1772,427,1957,424,2070,416,1990,406,2182,431,2008,434,1916,418,2397,412,2114,404,1778,409,1641,412,2186,406,1773,398,1785,397,2217,385,2153,390,1895,413,2475,413,1793,401,2308,397,2051,397,1898,409,2142,419,1874,424,1560,428,1808,430,1575,424,1525,433,1997,456,1753,459,1623,446,2251,441,1890),dim=c(2,61),dimnames=list(c('wkl','bvg'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('wkl','bvg'),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
wkl bvg M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 462 1919 1 0 0 0 0 0 0 0 0 0 0
2 455 1911 0 1 0 0 0 0 0 0 0 0 0
3 461 1870 0 0 1 0 0 0 0 0 0 0 0
4 461 2263 0 0 0 1 0 0 0 0 0 0 0
5 463 1802 0 0 0 0 1 0 0 0 0 0 0
6 462 1863 0 0 0 0 0 1 0 0 0 0 0
7 456 1989 0 0 0 0 0 0 1 0 0 0 0
8 455 2197 0 0 0 0 0 0 0 1 0 0 0
9 456 2409 0 0 0 0 0 0 0 0 1 0 0
10 472 2502 0 0 0 0 0 0 0 0 0 1 0
11 472 2593 0 0 0 0 0 0 0 0 0 0 1
12 471 2598 0 0 0 0 0 0 0 0 0 0 0
13 465 2053 1 0 0 0 0 0 0 0 0 0 0
14 459 2213 0 1 0 0 0 0 0 0 0 0 0
15 465 2238 0 0 1 0 0 0 0 0 0 0 0
16 468 2359 0 0 0 1 0 0 0 0 0 0 0
17 467 2151 0 0 0 0 1 0 0 0 0 0 0
18 463 2474 0 0 0 0 0 1 0 0 0 0 0
19 460 3079 0 0 0 0 0 0 1 0 0 0 0
20 462 2312 0 0 0 0 0 0 0 1 0 0 0
21 461 2565 0 0 0 0 0 0 0 0 1 0 0
22 476 1972 0 0 0 0 0 0 0 0 0 1 0
23 476 2484 0 0 0 0 0 0 0 0 0 0 1
24 471 2202 0 0 0 0 0 0 0 0 0 0 0
25 453 2151 1 0 0 0 0 0 0 0 0 0 0
26 443 1976 0 1 0 0 0 0 0 0 0 0 0
27 442 2012 0 0 1 0 0 0 0 0 0 0 0
28 444 2114 0 0 0 1 0 0 0 0 0 0 0
29 438 1772 0 0 0 0 1 0 0 0 0 0 0
30 427 1957 0 0 0 0 0 1 0 0 0 0 0
31 424 2070 0 0 0 0 0 0 1 0 0 0 0
32 416 1990 0 0 0 0 0 0 0 1 0 0 0
33 406 2182 0 0 0 0 0 0 0 0 1 0 0
34 431 2008 0 0 0 0 0 0 0 0 0 1 0
35 434 1916 0 0 0 0 0 0 0 0 0 0 1
36 418 2397 0 0 0 0 0 0 0 0 0 0 0
37 412 2114 1 0 0 0 0 0 0 0 0 0 0
38 404 1778 0 1 0 0 0 0 0 0 0 0 0
39 409 1641 0 0 1 0 0 0 0 0 0 0 0
40 412 2186 0 0 0 1 0 0 0 0 0 0 0
41 406 1773 0 0 0 0 1 0 0 0 0 0 0
42 398 1785 0 0 0 0 0 1 0 0 0 0 0
43 397 2217 0 0 0 0 0 0 1 0 0 0 0
44 385 2153 0 0 0 0 0 0 0 1 0 0 0
45 390 1895 0 0 0 0 0 0 0 0 1 0 0
46 413 2475 0 0 0 0 0 0 0 0 0 1 0
47 413 1793 0 0 0 0 0 0 0 0 0 0 1
48 401 2308 0 0 0 0 0 0 0 0 0 0 0
49 397 2051 1 0 0 0 0 0 0 0 0 0 0
50 397 1898 0 1 0 0 0 0 0 0 0 0 0
51 409 2142 0 0 1 0 0 0 0 0 0 0 0
52 419 1874 0 0 0 1 0 0 0 0 0 0 0
53 424 1560 0 0 0 0 1 0 0 0 0 0 0
54 428 1808 0 0 0 0 0 1 0 0 0 0 0
55 430 1575 0 0 0 0 0 0 1 0 0 0 0
56 424 1525 0 0 0 0 0 0 0 1 0 0 0
57 433 1997 0 0 0 0 0 0 0 0 1 0 0
58 456 1753 0 0 0 0 0 0 0 0 0 1 0
59 459 1623 0 0 0 0 0 0 0 0 0 0 1
60 446 2251 0 0 0 0 0 0 0 0 0 0 0
61 441 1890 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) bvg M1 M2 M3 M4
356.35906 0.03617 8.56293 4.52299 9.20429 6.34448
M5 M6 M7 M8 M9 M10
17.71688 7.72003 -2.02485 -1.57778 -7.07845 15.76659
M11
19.14397
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-48.644 -21.191 7.718 19.588 34.996
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 356.35906 32.20814 11.064 8.34e-15 ***
bvg 0.03617 0.01276 2.835 0.0067 **
M1 8.56293 16.38610 0.523 0.6037
M2 4.52299 17.32292 0.261 0.7951
M3 9.20429 17.23117 0.534 0.5957
M4 6.34448 16.74971 0.379 0.7065
M5 17.71688 17.94320 0.987 0.3284
M6 7.72003 17.24242 0.448 0.6564
M7 -2.02485 16.70312 -0.121 0.9040
M8 -1.57778 17.05253 -0.093 0.9267
M9 -7.07845 16.66779 -0.425 0.6730
M10 15.76659 16.78322 0.939 0.3522
M11 19.14397 16.92240 1.131 0.2636
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 26.2 on 48 degrees of freedom
Multiple R-squared: 0.2014, Adjusted R-squared: 0.00173
F-statistic: 1.009 on 12 and 48 DF, p-value: 0.456
> 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,] 1.113839e-03 2.227678e-03 0.9988862
[2,] 9.296476e-05 1.859295e-04 0.9999070
[3,] 2.850280e-05 5.700560e-05 0.9999715
[4,] 3.552727e-06 7.105454e-06 0.9999964
[5,] 2.041318e-06 4.082636e-06 0.9999980
[6,] 6.050324e-07 1.210065e-06 0.9999994
[7,] 3.573275e-07 7.146551e-07 0.9999996
[8,] 1.774167e-07 3.548335e-07 0.9999998
[9,] 5.331124e-08 1.066225e-07 0.9999999
[10,] 5.027090e-07 1.005418e-06 0.9999995
[11,] 7.041744e-06 1.408349e-05 0.9999930
[12,] 2.316284e-04 4.632568e-04 0.9997684
[13,] 1.237819e-03 2.475638e-03 0.9987622
[14,] 8.152556e-03 1.630511e-02 0.9918474
[15,] 4.675738e-02 9.351475e-02 0.9532426
[16,] 8.048430e-02 1.609686e-01 0.9195157
[17,] 1.830498e-01 3.660996e-01 0.8169502
[18,] 3.228721e-01 6.457442e-01 0.6771279
[19,] 3.672622e-01 7.345243e-01 0.6327378
[20,] 3.146489e-01 6.292978e-01 0.6853511
[21,] 4.269733e-01 8.539467e-01 0.5730267
[22,] 4.943092e-01 9.886185e-01 0.5056908
[23,] 4.706093e-01 9.412186e-01 0.5293907
[24,] 4.788972e-01 9.577944e-01 0.5211028
[25,] 4.586529e-01 9.173058e-01 0.5413471
[26,] 4.122200e-01 8.244399e-01 0.5877800
[27,] 4.197505e-01 8.395010e-01 0.5802495
[28,] 3.443112e-01 6.886223e-01 0.6556888
[29,] 2.892290e-01 5.784580e-01 0.7107710
[30,] 3.838590e-01 7.677180e-01 0.6161410
> postscript(file="/var/www/html/rcomp/tmp/130wq1258740726.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/2qsc51258740726.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/3mswu1258740726.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/4zdvs1258740726.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/5c5xw1258740726.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
27.6693876 24.9986771 27.8003097 16.4456407 23.7472240 30.5377525
7 8 9 10 11 12
29.7253255 20.7550629 19.5878685 9.3791006 2.7103228 20.6734501
13 14 15 16 17 18
25.8227195 18.0755892 18.4900569 19.9734008 15.1241854 9.4383926
19 20 21 22 23 24
-5.6990644 23.5956089 18.9454787 32.5487581 10.6527618 34.9964394
25 26 27 28 29 30
10.2781413 10.6476814 3.6642882 4.8348463 -0.1677011 -7.8621490
31 32 33 34 35 36
-5.2043768 -10.7579199 -22.2017310 -13.7533319 -10.8031524 -25.0565478
37 38 39 40 41 42
-29.3835996 -21.1908239 -15.9169515 -29.7693336 -32.2038702 -30.6410526
43 44 45 46 47 48
-37.5212441 -47.6534938 -27.8211807 -48.6443320 -27.3543451 -38.8374921
49 50 51 52 53 54
-42.1049422 -32.5311238 -34.0377032 -11.4845541 -6.4998381 -1.4729434
55 56 57 58 59 60
18.6993599 14.0607419 11.4895645 20.4698052 24.7944130 8.2241503
61
7.7182934
> postscript(file="/var/www/html/rcomp/tmp/64njb1258740726.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 27.6693876 NA
1 24.9986771 27.6693876
2 27.8003097 24.9986771
3 16.4456407 27.8003097
4 23.7472240 16.4456407
5 30.5377525 23.7472240
6 29.7253255 30.5377525
7 20.7550629 29.7253255
8 19.5878685 20.7550629
9 9.3791006 19.5878685
10 2.7103228 9.3791006
11 20.6734501 2.7103228
12 25.8227195 20.6734501
13 18.0755892 25.8227195
14 18.4900569 18.0755892
15 19.9734008 18.4900569
16 15.1241854 19.9734008
17 9.4383926 15.1241854
18 -5.6990644 9.4383926
19 23.5956089 -5.6990644
20 18.9454787 23.5956089
21 32.5487581 18.9454787
22 10.6527618 32.5487581
23 34.9964394 10.6527618
24 10.2781413 34.9964394
25 10.6476814 10.2781413
26 3.6642882 10.6476814
27 4.8348463 3.6642882
28 -0.1677011 4.8348463
29 -7.8621490 -0.1677011
30 -5.2043768 -7.8621490
31 -10.7579199 -5.2043768
32 -22.2017310 -10.7579199
33 -13.7533319 -22.2017310
34 -10.8031524 -13.7533319
35 -25.0565478 -10.8031524
36 -29.3835996 -25.0565478
37 -21.1908239 -29.3835996
38 -15.9169515 -21.1908239
39 -29.7693336 -15.9169515
40 -32.2038702 -29.7693336
41 -30.6410526 -32.2038702
42 -37.5212441 -30.6410526
43 -47.6534938 -37.5212441
44 -27.8211807 -47.6534938
45 -48.6443320 -27.8211807
46 -27.3543451 -48.6443320
47 -38.8374921 -27.3543451
48 -42.1049422 -38.8374921
49 -32.5311238 -42.1049422
50 -34.0377032 -32.5311238
51 -11.4845541 -34.0377032
52 -6.4998381 -11.4845541
53 -1.4729434 -6.4998381
54 18.6993599 -1.4729434
55 14.0607419 18.6993599
56 11.4895645 14.0607419
57 20.4698052 11.4895645
58 24.7944130 20.4698052
59 8.2241503 24.7944130
60 7.7182934 8.2241503
61 NA 7.7182934
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 24.9986771 27.6693876
[2,] 27.8003097 24.9986771
[3,] 16.4456407 27.8003097
[4,] 23.7472240 16.4456407
[5,] 30.5377525 23.7472240
[6,] 29.7253255 30.5377525
[7,] 20.7550629 29.7253255
[8,] 19.5878685 20.7550629
[9,] 9.3791006 19.5878685
[10,] 2.7103228 9.3791006
[11,] 20.6734501 2.7103228
[12,] 25.8227195 20.6734501
[13,] 18.0755892 25.8227195
[14,] 18.4900569 18.0755892
[15,] 19.9734008 18.4900569
[16,] 15.1241854 19.9734008
[17,] 9.4383926 15.1241854
[18,] -5.6990644 9.4383926
[19,] 23.5956089 -5.6990644
[20,] 18.9454787 23.5956089
[21,] 32.5487581 18.9454787
[22,] 10.6527618 32.5487581
[23,] 34.9964394 10.6527618
[24,] 10.2781413 34.9964394
[25,] 10.6476814 10.2781413
[26,] 3.6642882 10.6476814
[27,] 4.8348463 3.6642882
[28,] -0.1677011 4.8348463
[29,] -7.8621490 -0.1677011
[30,] -5.2043768 -7.8621490
[31,] -10.7579199 -5.2043768
[32,] -22.2017310 -10.7579199
[33,] -13.7533319 -22.2017310
[34,] -10.8031524 -13.7533319
[35,] -25.0565478 -10.8031524
[36,] -29.3835996 -25.0565478
[37,] -21.1908239 -29.3835996
[38,] -15.9169515 -21.1908239
[39,] -29.7693336 -15.9169515
[40,] -32.2038702 -29.7693336
[41,] -30.6410526 -32.2038702
[42,] -37.5212441 -30.6410526
[43,] -47.6534938 -37.5212441
[44,] -27.8211807 -47.6534938
[45,] -48.6443320 -27.8211807
[46,] -27.3543451 -48.6443320
[47,] -38.8374921 -27.3543451
[48,] -42.1049422 -38.8374921
[49,] -32.5311238 -42.1049422
[50,] -34.0377032 -32.5311238
[51,] -11.4845541 -34.0377032
[52,] -6.4998381 -11.4845541
[53,] -1.4729434 -6.4998381
[54,] 18.6993599 -1.4729434
[55,] 14.0607419 18.6993599
[56,] 11.4895645 14.0607419
[57,] 20.4698052 11.4895645
[58,] 24.7944130 20.4698052
[59,] 8.2241503 24.7944130
[60,] 7.7182934 8.2241503
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 24.9986771 27.6693876
2 27.8003097 24.9986771
3 16.4456407 27.8003097
4 23.7472240 16.4456407
5 30.5377525 23.7472240
6 29.7253255 30.5377525
7 20.7550629 29.7253255
8 19.5878685 20.7550629
9 9.3791006 19.5878685
10 2.7103228 9.3791006
11 20.6734501 2.7103228
12 25.8227195 20.6734501
13 18.0755892 25.8227195
14 18.4900569 18.0755892
15 19.9734008 18.4900569
16 15.1241854 19.9734008
17 9.4383926 15.1241854
18 -5.6990644 9.4383926
19 23.5956089 -5.6990644
20 18.9454787 23.5956089
21 32.5487581 18.9454787
22 10.6527618 32.5487581
23 34.9964394 10.6527618
24 10.2781413 34.9964394
25 10.6476814 10.2781413
26 3.6642882 10.6476814
27 4.8348463 3.6642882
28 -0.1677011 4.8348463
29 -7.8621490 -0.1677011
30 -5.2043768 -7.8621490
31 -10.7579199 -5.2043768
32 -22.2017310 -10.7579199
33 -13.7533319 -22.2017310
34 -10.8031524 -13.7533319
35 -25.0565478 -10.8031524
36 -29.3835996 -25.0565478
37 -21.1908239 -29.3835996
38 -15.9169515 -21.1908239
39 -29.7693336 -15.9169515
40 -32.2038702 -29.7693336
41 -30.6410526 -32.2038702
42 -37.5212441 -30.6410526
43 -47.6534938 -37.5212441
44 -27.8211807 -47.6534938
45 -48.6443320 -27.8211807
46 -27.3543451 -48.6443320
47 -38.8374921 -27.3543451
48 -42.1049422 -38.8374921
49 -32.5311238 -42.1049422
50 -34.0377032 -32.5311238
51 -11.4845541 -34.0377032
52 -6.4998381 -11.4845541
53 -1.4729434 -6.4998381
54 18.6993599 -1.4729434
55 14.0607419 18.6993599
56 11.4895645 14.0607419
57 20.4698052 11.4895645
58 24.7944130 20.4698052
59 8.2241503 24.7944130
60 7.7182934 8.2241503
> 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/7mqiw1258740726.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/8c2yy1258740726.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/9gdp21258740726.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/10uyxe1258740726.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/11yea71258740726.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/1208ck1258740726.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/13clg21258740727.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/14i5ix1258740727.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/1567hq1258740727.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/161q331258740727.tab")
+ }
>
> system("convert tmp/130wq1258740726.ps tmp/130wq1258740726.png")
> system("convert tmp/2qsc51258740726.ps tmp/2qsc51258740726.png")
> system("convert tmp/3mswu1258740726.ps tmp/3mswu1258740726.png")
> system("convert tmp/4zdvs1258740726.ps tmp/4zdvs1258740726.png")
> system("convert tmp/5c5xw1258740726.ps tmp/5c5xw1258740726.png")
> system("convert tmp/64njb1258740726.ps tmp/64njb1258740726.png")
> system("convert tmp/7mqiw1258740726.ps tmp/7mqiw1258740726.png")
> system("convert tmp/8c2yy1258740726.ps tmp/8c2yy1258740726.png")
> system("convert tmp/9gdp21258740726.ps tmp/9gdp21258740726.png")
> system("convert tmp/10uyxe1258740726.ps tmp/10uyxe1258740726.png")
>
>
> proc.time()
user system elapsed
2.383 1.560 2.794