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(441,1919,449,1911,452,1870,462,2263,455,1802,461,1863,461,1989,463,2197,462,2409,456,2502,455,2593,456,2598,472,2053,472,2213,471,2238,465,2359,459,2151,465,2474,468,3079,467,2312,463,2565,460,1972,462,2484,461,2202,476,2151,476,1976,471,2012,453,2114,443,1772,442,1957,444,2070,438,1990,427,2182,424,2008,416,1916,406,2397,431,2114,434,1778,418,1641,412,2186,404,1773,409,1785,412,2217,406,2153,398,1895,397,2475,385,1793,390,2308,413,2051,413,1898,401,2142,397,1874,397,1560,409,1808,419,1575,424,1525,428,1997,430,1753,424,1623,433,2251,456,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 441 1919 1 0 0 0 0 0 0 0 0 0 0
2 449 1911 0 1 0 0 0 0 0 0 0 0 0
3 452 1870 0 0 1 0 0 0 0 0 0 0 0
4 462 2263 0 0 0 1 0 0 0 0 0 0 0
5 455 1802 0 0 0 0 1 0 0 0 0 0 0
6 461 1863 0 0 0 0 0 1 0 0 0 0 0
7 461 1989 0 0 0 0 0 0 1 0 0 0 0
8 463 2197 0 0 0 0 0 0 0 1 0 0 0
9 462 2409 0 0 0 0 0 0 0 0 1 0 0
10 456 2502 0 0 0 0 0 0 0 0 0 1 0
11 455 2593 0 0 0 0 0 0 0 0 0 0 1
12 456 2598 0 0 0 0 0 0 0 0 0 0 0
13 472 2053 1 0 0 0 0 0 0 0 0 0 0
14 472 2213 0 1 0 0 0 0 0 0 0 0 0
15 471 2238 0 0 1 0 0 0 0 0 0 0 0
16 465 2359 0 0 0 1 0 0 0 0 0 0 0
17 459 2151 0 0 0 0 1 0 0 0 0 0 0
18 465 2474 0 0 0 0 0 1 0 0 0 0 0
19 468 3079 0 0 0 0 0 0 1 0 0 0 0
20 467 2312 0 0 0 0 0 0 0 1 0 0 0
21 463 2565 0 0 0 0 0 0 0 0 1 0 0
22 460 1972 0 0 0 0 0 0 0 0 0 1 0
23 462 2484 0 0 0 0 0 0 0 0 0 0 1
24 461 2202 0 0 0 0 0 0 0 0 0 0 0
25 476 2151 1 0 0 0 0 0 0 0 0 0 0
26 476 1976 0 1 0 0 0 0 0 0 0 0 0
27 471 2012 0 0 1 0 0 0 0 0 0 0 0
28 453 2114 0 0 0 1 0 0 0 0 0 0 0
29 443 1772 0 0 0 0 1 0 0 0 0 0 0
30 442 1957 0 0 0 0 0 1 0 0 0 0 0
31 444 2070 0 0 0 0 0 0 1 0 0 0 0
32 438 1990 0 0 0 0 0 0 0 1 0 0 0
33 427 2182 0 0 0 0 0 0 0 0 1 0 0
34 424 2008 0 0 0 0 0 0 0 0 0 1 0
35 416 1916 0 0 0 0 0 0 0 0 0 0 1
36 406 2397 0 0 0 0 0 0 0 0 0 0 0
37 431 2114 1 0 0 0 0 0 0 0 0 0 0
38 434 1778 0 1 0 0 0 0 0 0 0 0 0
39 418 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 404 1773 0 0 0 0 1 0 0 0 0 0 0
42 409 1785 0 0 0 0 0 1 0 0 0 0 0
43 412 2217 0 0 0 0 0 0 1 0 0 0 0
44 406 2153 0 0 0 0 0 0 0 1 0 0 0
45 398 1895 0 0 0 0 0 0 0 0 1 0 0
46 397 2475 0 0 0 0 0 0 0 0 0 1 0
47 385 1793 0 0 0 0 0 0 0 0 0 0 1
48 390 2308 0 0 0 0 0 0 0 0 0 0 0
49 413 2051 1 0 0 0 0 0 0 0 0 0 0
50 413 1898 0 1 0 0 0 0 0 0 0 0 0
51 401 2142 0 0 1 0 0 0 0 0 0 0 0
52 397 1874 0 0 0 1 0 0 0 0 0 0 0
53 397 1560 0 0 0 0 1 0 0 0 0 0 0
54 409 1808 0 0 0 0 0 1 0 0 0 0 0
55 419 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 428 1997 0 0 0 0 0 0 0 0 1 0 0
58 430 1753 0 0 0 0 0 0 0 0 0 1 0
59 424 1623 0 0 0 0 0 0 0 0 0 0 1
60 433 2251 0 0 0 0 0 0 0 0 0 0 0
61 456 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
320.08343 0.04641 33.88867 37.97792 30.59913 17.51051
M5 M6 M7 M8 M9 M10
27.44224 25.34764 19.26675 25.05593 12.97150 13.90874
M11
11.70255
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-51.854 -20.338 5.893 15.900 38.724
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 320.08343 30.44337 10.514 4.78e-14 ***
bvg 0.04641 0.01206 3.848 0.000351 ***
M1 33.88867 15.48827 2.188 0.033568 *
M2 37.97792 16.37375 2.319 0.024672 *
M3 30.59913 16.28703 1.879 0.066360 .
M4 17.51051 15.83195 1.106 0.274229
M5 27.44224 16.96005 1.618 0.112205
M6 25.34764 16.29766 1.555 0.126446
M7 19.26675 15.78791 1.220 0.228295
M8 25.05593 16.11818 1.555 0.126631
M9 12.97150 15.75452 0.823 0.414381
M10 13.90874 15.86363 0.877 0.384979
M11 11.70255 15.99518 0.732 0.467951
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 24.76 on 48 degrees of freedom
Multiple R-squared: 0.2832, Adjusted R-squared: 0.104
F-statistic: 1.581 on 12 and 48 DF, p-value: 0.1295
> 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.0552073133 0.110414627 0.9447927
[2,] 0.0368852401 0.073770480 0.9631148
[3,] 0.0285916080 0.057183216 0.9714084
[4,] 0.0160703810 0.032140762 0.9839296
[5,] 0.0063479418 0.012695884 0.9936521
[6,] 0.0024474358 0.004894872 0.9975526
[7,] 0.0014522018 0.002904404 0.9985478
[8,] 0.0008402498 0.001680500 0.9991598
[9,] 0.0006196575 0.001239315 0.9993803
[10,] 0.0009414132 0.001882826 0.9990586
[11,] 0.0016323119 0.003264624 0.9983677
[12,] 0.0032222877 0.006444575 0.9967777
[13,] 0.0048132120 0.009626424 0.9951868
[14,] 0.0084738702 0.016947740 0.9915261
[15,] 0.0182025034 0.036405007 0.9817975
[16,] 0.0265348164 0.053069633 0.9734652
[17,] 0.0521583721 0.104316744 0.9478416
[18,] 0.1059460030 0.211892006 0.8940540
[19,] 0.1348990349 0.269798070 0.8651010
[20,] 0.1592408509 0.318481702 0.8407591
[21,] 0.3401738208 0.680347642 0.6598262
[22,] 0.3550003876 0.710000775 0.6449996
[23,] 0.3334067879 0.666813576 0.6665932
[24,] 0.3092819524 0.618563905 0.6907180
[25,] 0.4006517555 0.801303511 0.5993482
[26,] 0.4032648669 0.806529734 0.5967351
[27,] 0.3312506874 0.662501375 0.6687493
[28,] 0.3052945796 0.610589159 0.6947054
[29,] 0.2747037944 0.549407589 0.7252962
[30,] 0.2934527918 0.586905584 0.7065472
> postscript(file="/var/www/html/rcomp/tmp/1w0al1260819535.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/2iki01260819535.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/34g971260819535.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/4qxod1260819535.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/523ry1260819535.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
-2.0307499 2.2512728 14.5328228 19.3827576 23.8455253 29.1091766
7 8 9 10 11 12
29.3425506 15.9003240 17.1460681 5.8928010 2.8757774 15.3462869
13 14 15 16 17 18
22.7504593 11.2357891 16.4543527 17.9275045 11.6488240 4.7533471
19 20 21 22 23 24
-14.2431354 14.5633021 10.9062818 34.4895107 14.9343460 38.7242058
25 26 27 28 29 30
22.2023885 26.2346952 26.9427610 17.2976817 13.2377919 5.7467413
31 32 33 34 35 36
8.5834308 0.5069634 -7.3191147 -3.1812092 -4.7054067 -25.3255270
37 38 39 40 41 42
-21.0804827 -6.5763454 -8.8395422 -27.0437582 -25.8086170 -19.2709303
43 44 45 46 47 48
-30.2386755 -39.0576850 -22.9997644 -51.8541591 -29.9971137 -37.1951361
49 50 51 52 53 54
-36.1567229 -33.1454117 -49.0903942 -27.5641856 -22.9235242 -20.3383347
55 56 57 58 59 60
6.5558295 8.0870955 2.2665292 14.6530567 16.8923970 8.4501704
61
14.3151078
> postscript(file="/var/www/html/rcomp/tmp/6o8vg1260819535.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 -2.0307499 NA
1 2.2512728 -2.0307499
2 14.5328228 2.2512728
3 19.3827576 14.5328228
4 23.8455253 19.3827576
5 29.1091766 23.8455253
6 29.3425506 29.1091766
7 15.9003240 29.3425506
8 17.1460681 15.9003240
9 5.8928010 17.1460681
10 2.8757774 5.8928010
11 15.3462869 2.8757774
12 22.7504593 15.3462869
13 11.2357891 22.7504593
14 16.4543527 11.2357891
15 17.9275045 16.4543527
16 11.6488240 17.9275045
17 4.7533471 11.6488240
18 -14.2431354 4.7533471
19 14.5633021 -14.2431354
20 10.9062818 14.5633021
21 34.4895107 10.9062818
22 14.9343460 34.4895107
23 38.7242058 14.9343460
24 22.2023885 38.7242058
25 26.2346952 22.2023885
26 26.9427610 26.2346952
27 17.2976817 26.9427610
28 13.2377919 17.2976817
29 5.7467413 13.2377919
30 8.5834308 5.7467413
31 0.5069634 8.5834308
32 -7.3191147 0.5069634
33 -3.1812092 -7.3191147
34 -4.7054067 -3.1812092
35 -25.3255270 -4.7054067
36 -21.0804827 -25.3255270
37 -6.5763454 -21.0804827
38 -8.8395422 -6.5763454
39 -27.0437582 -8.8395422
40 -25.8086170 -27.0437582
41 -19.2709303 -25.8086170
42 -30.2386755 -19.2709303
43 -39.0576850 -30.2386755
44 -22.9997644 -39.0576850
45 -51.8541591 -22.9997644
46 -29.9971137 -51.8541591
47 -37.1951361 -29.9971137
48 -36.1567229 -37.1951361
49 -33.1454117 -36.1567229
50 -49.0903942 -33.1454117
51 -27.5641856 -49.0903942
52 -22.9235242 -27.5641856
53 -20.3383347 -22.9235242
54 6.5558295 -20.3383347
55 8.0870955 6.5558295
56 2.2665292 8.0870955
57 14.6530567 2.2665292
58 16.8923970 14.6530567
59 8.4501704 16.8923970
60 14.3151078 8.4501704
61 NA 14.3151078
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2.2512728 -2.0307499
[2,] 14.5328228 2.2512728
[3,] 19.3827576 14.5328228
[4,] 23.8455253 19.3827576
[5,] 29.1091766 23.8455253
[6,] 29.3425506 29.1091766
[7,] 15.9003240 29.3425506
[8,] 17.1460681 15.9003240
[9,] 5.8928010 17.1460681
[10,] 2.8757774 5.8928010
[11,] 15.3462869 2.8757774
[12,] 22.7504593 15.3462869
[13,] 11.2357891 22.7504593
[14,] 16.4543527 11.2357891
[15,] 17.9275045 16.4543527
[16,] 11.6488240 17.9275045
[17,] 4.7533471 11.6488240
[18,] -14.2431354 4.7533471
[19,] 14.5633021 -14.2431354
[20,] 10.9062818 14.5633021
[21,] 34.4895107 10.9062818
[22,] 14.9343460 34.4895107
[23,] 38.7242058 14.9343460
[24,] 22.2023885 38.7242058
[25,] 26.2346952 22.2023885
[26,] 26.9427610 26.2346952
[27,] 17.2976817 26.9427610
[28,] 13.2377919 17.2976817
[29,] 5.7467413 13.2377919
[30,] 8.5834308 5.7467413
[31,] 0.5069634 8.5834308
[32,] -7.3191147 0.5069634
[33,] -3.1812092 -7.3191147
[34,] -4.7054067 -3.1812092
[35,] -25.3255270 -4.7054067
[36,] -21.0804827 -25.3255270
[37,] -6.5763454 -21.0804827
[38,] -8.8395422 -6.5763454
[39,] -27.0437582 -8.8395422
[40,] -25.8086170 -27.0437582
[41,] -19.2709303 -25.8086170
[42,] -30.2386755 -19.2709303
[43,] -39.0576850 -30.2386755
[44,] -22.9997644 -39.0576850
[45,] -51.8541591 -22.9997644
[46,] -29.9971137 -51.8541591
[47,] -37.1951361 -29.9971137
[48,] -36.1567229 -37.1951361
[49,] -33.1454117 -36.1567229
[50,] -49.0903942 -33.1454117
[51,] -27.5641856 -49.0903942
[52,] -22.9235242 -27.5641856
[53,] -20.3383347 -22.9235242
[54,] 6.5558295 -20.3383347
[55,] 8.0870955 6.5558295
[56,] 2.2665292 8.0870955
[57,] 14.6530567 2.2665292
[58,] 16.8923970 14.6530567
[59,] 8.4501704 16.8923970
[60,] 14.3151078 8.4501704
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2.2512728 -2.0307499
2 14.5328228 2.2512728
3 19.3827576 14.5328228
4 23.8455253 19.3827576
5 29.1091766 23.8455253
6 29.3425506 29.1091766
7 15.9003240 29.3425506
8 17.1460681 15.9003240
9 5.8928010 17.1460681
10 2.8757774 5.8928010
11 15.3462869 2.8757774
12 22.7504593 15.3462869
13 11.2357891 22.7504593
14 16.4543527 11.2357891
15 17.9275045 16.4543527
16 11.6488240 17.9275045
17 4.7533471 11.6488240
18 -14.2431354 4.7533471
19 14.5633021 -14.2431354
20 10.9062818 14.5633021
21 34.4895107 10.9062818
22 14.9343460 34.4895107
23 38.7242058 14.9343460
24 22.2023885 38.7242058
25 26.2346952 22.2023885
26 26.9427610 26.2346952
27 17.2976817 26.9427610
28 13.2377919 17.2976817
29 5.7467413 13.2377919
30 8.5834308 5.7467413
31 0.5069634 8.5834308
32 -7.3191147 0.5069634
33 -3.1812092 -7.3191147
34 -4.7054067 -3.1812092
35 -25.3255270 -4.7054067
36 -21.0804827 -25.3255270
37 -6.5763454 -21.0804827
38 -8.8395422 -6.5763454
39 -27.0437582 -8.8395422
40 -25.8086170 -27.0437582
41 -19.2709303 -25.8086170
42 -30.2386755 -19.2709303
43 -39.0576850 -30.2386755
44 -22.9997644 -39.0576850
45 -51.8541591 -22.9997644
46 -29.9971137 -51.8541591
47 -37.1951361 -29.9971137
48 -36.1567229 -37.1951361
49 -33.1454117 -36.1567229
50 -49.0903942 -33.1454117
51 -27.5641856 -49.0903942
52 -22.9235242 -27.5641856
53 -20.3383347 -22.9235242
54 6.5558295 -20.3383347
55 8.0870955 6.5558295
56 2.2665292 8.0870955
57 14.6530567 2.2665292
58 16.8923970 14.6530567
59 8.4501704 16.8923970
60 14.3151078 8.4501704
> 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/7z2rb1260819535.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/8abfk1260819535.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/97pgt1260819535.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/107tfv1260819535.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/11y6fw1260819535.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/12fwdt1260819535.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/13x3a71260819535.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/14ppov1260819535.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/15cesb1260819535.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/16228n1260819535.tab")
+ }
>
> try(system("convert tmp/1w0al1260819535.ps tmp/1w0al1260819535.png",intern=TRUE))
character(0)
> try(system("convert tmp/2iki01260819535.ps tmp/2iki01260819535.png",intern=TRUE))
character(0)
> try(system("convert tmp/34g971260819535.ps tmp/34g971260819535.png",intern=TRUE))
character(0)
> try(system("convert tmp/4qxod1260819535.ps tmp/4qxod1260819535.png",intern=TRUE))
character(0)
> try(system("convert tmp/523ry1260819535.ps tmp/523ry1260819535.png",intern=TRUE))
character(0)
> try(system("convert tmp/6o8vg1260819535.ps tmp/6o8vg1260819535.png",intern=TRUE))
character(0)
> try(system("convert tmp/7z2rb1260819535.ps tmp/7z2rb1260819535.png",intern=TRUE))
character(0)
> try(system("convert tmp/8abfk1260819535.ps tmp/8abfk1260819535.png",intern=TRUE))
character(0)
> try(system("convert tmp/97pgt1260819535.ps tmp/97pgt1260819535.png",intern=TRUE))
character(0)
> try(system("convert tmp/107tfv1260819535.ps tmp/107tfv1260819535.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.456 1.571 3.614