R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(1.071,630,776,1.762,388,222,414,864,1.068,1.234,449,547,638,919,668,808,673,724,503,837,853,527,845,730,861,626,612,667,871,558,1.077,424,312,298,744,710,255,804,625,940,338,250,518,635,303,-1.284,1.069,1.507,403,487,388,752,360,107,-132,481,868,1.052,68,62,549,262,346,666,251,222,887,120,76,-459,548,991,103,408,565,315,409,338,703,182,172,424,278,351,338,271,430,775,107,121,586,215,200,624,196,174,615,236,141,654,181,148,567,206,175,601,162,158,642,166,100,712,91,80,566,145,145,530,170,150,296,309,245,-109,404,552,427,199,205,-1.859,965,1.715,251,255,300,421,194,176,195,278,311,-1.019,681,1.120,504,153,111,448,197,119,438,132,167,467,141,111,190,215,303,696,NA,NA,458,105,94,8,224,412,-39,314,357,-42,469,195,355,160,99,382,138,76,271,157,159,66,232,288,435,84,55,-453,451,566,326,132,101,295,116,143,258,160,135,259,168,125,-126,256,419,92,182,274),dim=c(3,70),dimnames=list(c('weekdagen','zaterdag','zondag'),1:70))
> y <- array(NA,dim=c(3,70),dimnames=list(c('weekdagen','zaterdag','zondag'),1:70))
> 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
> 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
weekdagen zaterdag zondag
1 1.071 630.000 776.000
2 1.762 388.000 222.000
3 414.000 864.000 1.068
4 1.234 449.000 547.000
5 638.000 919.000 668.000
6 808.000 673.000 724.000
7 503.000 837.000 853.000
8 527.000 845.000 730.000
9 861.000 626.000 612.000
10 667.000 871.000 558.000
11 1.077 424.000 312.000
12 298.000 744.000 710.000
13 255.000 804.000 625.000
14 940.000 338.000 250.000
15 518.000 635.000 303.000
16 -1.284 1.069 1.507
17 403.000 487.000 388.000
18 752.000 360.000 107.000
19 -132.000 481.000 868.000
20 1.052 68.000 62.000
21 549.000 262.000 346.000
22 666.000 251.000 222.000
23 887.000 120.000 76.000
24 -459.000 548.000 991.000
25 103.000 408.000 565.000
26 315.000 409.000 338.000
27 703.000 182.000 172.000
28 424.000 278.000 351.000
29 338.000 271.000 430.000
30 775.000 107.000 121.000
31 586.000 215.000 200.000
32 624.000 196.000 174.000
33 615.000 236.000 141.000
34 654.000 181.000 148.000
35 567.000 206.000 175.000
36 601.000 162.000 158.000
37 642.000 166.000 100.000
38 712.000 91.000 80.000
39 566.000 145.000 145.000
40 530.000 170.000 150.000
41 296.000 309.000 245.000
42 -109.000 404.000 552.000
43 427.000 199.000 205.000
44 -1.859 965.000 1.715
45 251.000 255.000 300.000
46 421.000 194.000 176.000
47 195.000 278.000 311.000
48 -1.019 681.000 1.120
49 504.000 153.000 111.000
50 448.000 197.000 119.000
51 438.000 132.000 167.000
52 467.000 141.000 111.000
53 190.000 215.000 303.000
54 696.000 NA NA
55 458.000 105.000 94.000
56 8.000 224.000 412.000
57 -39.000 314.000 357.000
58 -42.000 469.000 195.000
59 355.000 160.000 99.000
60 382.000 138.000 76.000
61 271.000 157.000 159.000
62 66.000 232.000 288.000
63 435.000 84.000 55.000
64 -453.000 451.000 566.000
65 326.000 132.000 101.000
66 295.000 116.000 143.000
67 258.000 160.000 135.000
68 259.000 168.000 125.000
69 -126.000 256.000 419.000
70 92.000 182.000 274.000
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) zaterdag zondag
432.2499 0.1113 -0.4154
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-700.29 -178.89 40.28 216.36 613.35
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 432.2499 64.9328 6.657 6.57e-09 ***
zaterdag 0.1113 0.1875 0.593 0.5550
zondag -0.4154 0.1913 -2.172 0.0335 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 300.4 on 66 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.07751, Adjusted R-squared: 0.04956
F-statistic: 2.773 on 2 and 66 DF, p-value: 0.06978
> 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.6831436 6.337127e-01 3.168564e-01
[2,] 0.5510058 8.979884e-01 4.489942e-01
[3,] 0.4235374 8.470749e-01 5.764626e-01
[4,] 0.7555480 4.889041e-01 2.444520e-01
[5,] 0.7231690 5.536621e-01 2.768310e-01
[6,] 0.6430914 7.138172e-01 3.569086e-01
[7,] 0.6180503 7.638994e-01 3.819497e-01
[8,] 0.6290307 7.419386e-01 3.709693e-01
[9,] 0.9682516 6.349672e-02 3.174836e-02
[10,] 0.9652457 6.950854e-02 3.475427e-02
[11,] 0.9787679 4.246424e-02 2.123212e-02
[12,] 0.9753344 4.933127e-02 2.466563e-02
[13,] 0.9897431 2.051373e-02 1.025687e-02
[14,] 0.9916443 1.671142e-02 8.355712e-03
[15,] 0.9971763 5.647439e-03 2.823719e-03
[16,] 0.9983469 3.306275e-03 1.653138e-03
[17,] 0.9992419 1.516273e-03 7.581367e-04
[18,] 0.9998515 2.969034e-04 1.484517e-04
[19,] 0.9999415 1.169144e-04 5.845721e-05
[20,] 0.9999340 1.320760e-04 6.603799e-05
[21,] 0.9999291 1.418528e-04 7.092639e-05
[22,] 0.9999602 7.965828e-05 3.982914e-05
[23,] 0.9999685 6.298428e-05 3.149214e-05
[24,] 0.9999800 4.008792e-05 2.004396e-05
[25,] 0.9999892 2.159062e-05 1.079531e-05
[26,] 0.9999918 1.648259e-05 8.241295e-06
[27,] 0.9999942 1.162739e-05 5.813697e-06
[28,] 0.9999957 8.615801e-06 4.307901e-06
[29,] 0.9999975 4.922653e-06 2.461326e-06
[30,] 0.9999982 3.603352e-06 1.801676e-06
[31,] 0.9999987 2.695274e-06 1.347637e-06
[32,] 0.9999989 2.131336e-06 1.065668e-06
[33,] 0.9999993 1.375919e-06 6.879596e-07
[34,] 0.9999994 1.192050e-06 5.960249e-07
[35,] 0.9999995 1.000987e-06 5.004935e-07
[36,] 0.9999995 9.844031e-07 4.922015e-07
[37,] 0.9999997 5.540085e-07 2.770043e-07
[38,] 0.9999998 3.932918e-07 1.966459e-07
[39,] 1.0000000 8.145706e-08 4.072853e-08
[40,] 1.0000000 5.536199e-08 2.768100e-08
[41,] 1.0000000 4.595899e-08 2.297949e-08
[42,] 1.0000000 2.259202e-08 1.129601e-08
[43,] 1.0000000 5.237265e-08 2.618633e-08
[44,] 1.0000000 4.389938e-08 2.194969e-08
[45,] 1.0000000 2.289073e-08 1.144536e-08
[46,] 1.0000000 1.251937e-08 6.259683e-09
[47,] 1.0000000 5.576512e-09 2.788256e-09
[48,] 1.0000000 9.788477e-10 4.894238e-10
[49,] 1.0000000 8.029461e-10 4.014730e-10
[50,] 1.0000000 2.024627e-10 1.012314e-10
[51,] 1.0000000 9.082579e-11 4.541289e-11
[52,] 1.0000000 1.018681e-09 5.093407e-10
[53,] 1.0000000 3.371759e-09 1.685879e-09
[54,] 1.0000000 2.210670e-08 1.105335e-08
[55,] 0.9999999 1.208617e-07 6.043084e-08
[56,] 0.9999999 2.006245e-07 1.003123e-07
[57,] 0.9999999 1.775112e-07 8.875558e-08
[58,] 0.9999982 3.609050e-06 1.804525e-06
[59,] 0.9999496 1.007438e-04 5.037188e-05
> postscript(file="/var/wessaorg/rcomp/tmp/14tfp1322147393.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)
Warning message:
In x[, 1] - mysum$resid :
longer object length is not a multiple of shorter object length
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/2j6pj1322147393.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/wessaorg/rcomp/tmp/3p1am1322147393.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/wessaorg/rcomp/tmp/4jfzd1322147393.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/wessaorg/rcomp/tmp/5n4wi1322147393.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 = 69
Frequency = 1
1 2 3 4 5 6
-178.894159 -381.432612 -113.946030 -253.728062 381.008703 601.646872
7 8 9 10 11 12
331.990739 304.000602 613.346663 369.650660 -348.733229 77.930236
13 14 15 16 17 18
-7.059117 574.001546 140.972154 -433.026757 77.753528 324.144652
19 20 21 22 23 24
-257.164384 -413.006746 231.341181 298.049784 472.971314 -540.519720
25 26 27 28 29 30
-139.921829 -22.339519 321.955280 106.638048 54.237253 381.112968
31 32 33 34 35 36
212.915782 242.228352 215.067690 263.095831 184.531070 216.364485
37 38 39 40 41 42
232.823480 302.860017 177.855317 141.150727 -68.848772 -356.877545
43 44 45 46 47 48
57.773383 -540.774806 -84.990459 40.281792 -138.979823 -508.580474
49 50 51 52 53 55
100.839943 43.267507 60.441695 65.175219 -140.293200 53.118451
56 57 58 59 60 61
-278.010959 -357.875098 -445.424785 -53.924329 -34.031599 -112.663703
62 63 64 65 66 67
-272.416542 16.252759 -700.291120 -78.977792 -90.748660 -135.968245
68 69 70
-140.012896 -412.663567 -246.669149
> postscript(file="/var/wessaorg/rcomp/tmp/67ssr1322147393.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 -178.894159 NA
1 -381.432612 -178.894159
2 -113.946030 -381.432612
3 -253.728062 -113.946030
4 381.008703 -253.728062
5 601.646872 381.008703
6 331.990739 601.646872
7 304.000602 331.990739
8 613.346663 304.000602
9 369.650660 613.346663
10 -348.733229 369.650660
11 77.930236 -348.733229
12 -7.059117 77.930236
13 574.001546 -7.059117
14 140.972154 574.001546
15 -433.026757 140.972154
16 77.753528 -433.026757
17 324.144652 77.753528
18 -257.164384 324.144652
19 -413.006746 -257.164384
20 231.341181 -413.006746
21 298.049784 231.341181
22 472.971314 298.049784
23 -540.519720 472.971314
24 -139.921829 -540.519720
25 -22.339519 -139.921829
26 321.955280 -22.339519
27 106.638048 321.955280
28 54.237253 106.638048
29 381.112968 54.237253
30 212.915782 381.112968
31 242.228352 212.915782
32 215.067690 242.228352
33 263.095831 215.067690
34 184.531070 263.095831
35 216.364485 184.531070
36 232.823480 216.364485
37 302.860017 232.823480
38 177.855317 302.860017
39 141.150727 177.855317
40 -68.848772 141.150727
41 -356.877545 -68.848772
42 57.773383 -356.877545
43 -540.774806 57.773383
44 -84.990459 -540.774806
45 40.281792 -84.990459
46 -138.979823 40.281792
47 -508.580474 -138.979823
48 100.839943 -508.580474
49 43.267507 100.839943
50 60.441695 43.267507
51 65.175219 60.441695
52 -140.293200 65.175219
53 53.118451 -140.293200
54 -278.010959 53.118451
55 -357.875098 -278.010959
56 -445.424785 -357.875098
57 -53.924329 -445.424785
58 -34.031599 -53.924329
59 -112.663703 -34.031599
60 -272.416542 -112.663703
61 16.252759 -272.416542
62 -700.291120 16.252759
63 -78.977792 -700.291120
64 -90.748660 -78.977792
65 -135.968245 -90.748660
66 -140.012896 -135.968245
67 -412.663567 -140.012896
68 -246.669149 -412.663567
69 NA -246.669149
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -381.432612 -178.894159
[2,] -113.946030 -381.432612
[3,] -253.728062 -113.946030
[4,] 381.008703 -253.728062
[5,] 601.646872 381.008703
[6,] 331.990739 601.646872
[7,] 304.000602 331.990739
[8,] 613.346663 304.000602
[9,] 369.650660 613.346663
[10,] -348.733229 369.650660
[11,] 77.930236 -348.733229
[12,] -7.059117 77.930236
[13,] 574.001546 -7.059117
[14,] 140.972154 574.001546
[15,] -433.026757 140.972154
[16,] 77.753528 -433.026757
[17,] 324.144652 77.753528
[18,] -257.164384 324.144652
[19,] -413.006746 -257.164384
[20,] 231.341181 -413.006746
[21,] 298.049784 231.341181
[22,] 472.971314 298.049784
[23,] -540.519720 472.971314
[24,] -139.921829 -540.519720
[25,] -22.339519 -139.921829
[26,] 321.955280 -22.339519
[27,] 106.638048 321.955280
[28,] 54.237253 106.638048
[29,] 381.112968 54.237253
[30,] 212.915782 381.112968
[31,] 242.228352 212.915782
[32,] 215.067690 242.228352
[33,] 263.095831 215.067690
[34,] 184.531070 263.095831
[35,] 216.364485 184.531070
[36,] 232.823480 216.364485
[37,] 302.860017 232.823480
[38,] 177.855317 302.860017
[39,] 141.150727 177.855317
[40,] -68.848772 141.150727
[41,] -356.877545 -68.848772
[42,] 57.773383 -356.877545
[43,] -540.774806 57.773383
[44,] -84.990459 -540.774806
[45,] 40.281792 -84.990459
[46,] -138.979823 40.281792
[47,] -508.580474 -138.979823
[48,] 100.839943 -508.580474
[49,] 43.267507 100.839943
[50,] 60.441695 43.267507
[51,] 65.175219 60.441695
[52,] -140.293200 65.175219
[53,] 53.118451 -140.293200
[54,] -278.010959 53.118451
[55,] -357.875098 -278.010959
[56,] -445.424785 -357.875098
[57,] -53.924329 -445.424785
[58,] -34.031599 -53.924329
[59,] -112.663703 -34.031599
[60,] -272.416542 -112.663703
[61,] 16.252759 -272.416542
[62,] -700.291120 16.252759
[63,] -78.977792 -700.291120
[64,] -90.748660 -78.977792
[65,] -135.968245 -90.748660
[66,] -140.012896 -135.968245
[67,] -412.663567 -140.012896
[68,] -246.669149 -412.663567
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -381.432612 -178.894159
2 -113.946030 -381.432612
3 -253.728062 -113.946030
4 381.008703 -253.728062
5 601.646872 381.008703
6 331.990739 601.646872
7 304.000602 331.990739
8 613.346663 304.000602
9 369.650660 613.346663
10 -348.733229 369.650660
11 77.930236 -348.733229
12 -7.059117 77.930236
13 574.001546 -7.059117
14 140.972154 574.001546
15 -433.026757 140.972154
16 77.753528 -433.026757
17 324.144652 77.753528
18 -257.164384 324.144652
19 -413.006746 -257.164384
20 231.341181 -413.006746
21 298.049784 231.341181
22 472.971314 298.049784
23 -540.519720 472.971314
24 -139.921829 -540.519720
25 -22.339519 -139.921829
26 321.955280 -22.339519
27 106.638048 321.955280
28 54.237253 106.638048
29 381.112968 54.237253
30 212.915782 381.112968
31 242.228352 212.915782
32 215.067690 242.228352
33 263.095831 215.067690
34 184.531070 263.095831
35 216.364485 184.531070
36 232.823480 216.364485
37 302.860017 232.823480
38 177.855317 302.860017
39 141.150727 177.855317
40 -68.848772 141.150727
41 -356.877545 -68.848772
42 57.773383 -356.877545
43 -540.774806 57.773383
44 -84.990459 -540.774806
45 40.281792 -84.990459
46 -138.979823 40.281792
47 -508.580474 -138.979823
48 100.839943 -508.580474
49 43.267507 100.839943
50 60.441695 43.267507
51 65.175219 60.441695
52 -140.293200 65.175219
53 53.118451 -140.293200
54 -278.010959 53.118451
55 -357.875098 -278.010959
56 -445.424785 -357.875098
57 -53.924329 -445.424785
58 -34.031599 -53.924329
59 -112.663703 -34.031599
60 -272.416542 -112.663703
61 16.252759 -272.416542
62 -700.291120 16.252759
63 -78.977792 -700.291120
64 -90.748660 -78.977792
65 -135.968245 -90.748660
66 -140.012896 -135.968245
67 -412.663567 -140.012896
68 -246.669149 -412.663567
> 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/wessaorg/rcomp/tmp/76ekl1322147393.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/wessaorg/rcomp/tmp/8l3xq1322147393.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/wessaorg/rcomp/tmp/9xcfh1322147393.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/wessaorg/rcomp/tmp/10j5081322147393.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/116avq1322147393.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/wessaorg/rcomp/tmp/12qu3e1322147393.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/wessaorg/rcomp/tmp/13zzx61322147393.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/wessaorg/rcomp/tmp/14spxy1322147393.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/wessaorg/rcomp/tmp/153ikp1322147393.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/wessaorg/rcomp/tmp/16lhfj1322147393.tab")
+ }
>
> try(system("convert tmp/14tfp1322147393.ps tmp/14tfp1322147393.png",intern=TRUE))
character(0)
> try(system("convert tmp/2j6pj1322147393.ps tmp/2j6pj1322147393.png",intern=TRUE))
character(0)
> try(system("convert tmp/3p1am1322147393.ps tmp/3p1am1322147393.png",intern=TRUE))
character(0)
> try(system("convert tmp/4jfzd1322147393.ps tmp/4jfzd1322147393.png",intern=TRUE))
character(0)
> try(system("convert tmp/5n4wi1322147393.ps tmp/5n4wi1322147393.png",intern=TRUE))
character(0)
> try(system("convert tmp/67ssr1322147393.ps tmp/67ssr1322147393.png",intern=TRUE))
character(0)
> try(system("convert tmp/76ekl1322147393.ps tmp/76ekl1322147393.png",intern=TRUE))
character(0)
> try(system("convert tmp/8l3xq1322147393.ps tmp/8l3xq1322147393.png",intern=TRUE))
character(0)
> try(system("convert tmp/9xcfh1322147393.ps tmp/9xcfh1322147393.png",intern=TRUE))
character(0)
> try(system("convert tmp/10j5081322147393.ps tmp/10j5081322147393.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.214 0.446 3.707