R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
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(2174.56,0,2196.72,0,2350.44,0,2440.25,0,2408.64,0,2472.81,0,2407.6,0,2454.62,0,2448.05,0,2497.84,0,2645.64,0,2756.76,0,2849.27,0,2921.44,0,2981.85,0,3080.58,0,3106.22,0,3119.31,0,3061.26,0,3097.31,0,3161.69,0,3257.16,0,3277.01,0,3295.32,0,3363.99,0,3494.17,0,3667.03,0,3813.06,0,3917.96,0,3895.51,0,3801.06,0,3570.12,0,3701.61,0,3862.27,0,3970.1,0,4138.52,0,4199.75,0,4290.89,0,4443.91,0,4502.64,0,4356.98,0,4591.27,0,4696.96,0,4621.4,0,4562.84,0,4202.52,0,4296.49,0,4435.23,0,4105.18,0,4116.68,0,3844.49,0,3720.98,0,3674.4,0,3857.62,0,3801.06,0,3504.37,1,3032.6,1,3047.03,1,2962.34,1,2197.82,1,2014.45,1),dim=c(2,61),dimnames=list(c('Bel_20','dummy'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Bel_20','dummy'),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 = '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
Bel_20 dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 2174.56 0 1 0 0 0 0 0 0 0 0 0 0 1
2 2196.72 0 0 1 0 0 0 0 0 0 0 0 0 2
3 2350.44 0 0 0 1 0 0 0 0 0 0 0 0 3
4 2440.25 0 0 0 0 1 0 0 0 0 0 0 0 4
5 2408.64 0 0 0 0 0 1 0 0 0 0 0 0 5
6 2472.81 0 0 0 0 0 0 1 0 0 0 0 0 6
7 2407.60 0 0 0 0 0 0 0 1 0 0 0 0 7
8 2454.62 0 0 0 0 0 0 0 0 1 0 0 0 8
9 2448.05 0 0 0 0 0 0 0 0 0 1 0 0 9
10 2497.84 0 0 0 0 0 0 0 0 0 0 1 0 10
11 2645.64 0 0 0 0 0 0 0 0 0 0 0 1 11
12 2756.76 0 0 0 0 0 0 0 0 0 0 0 0 12
13 2849.27 0 1 0 0 0 0 0 0 0 0 0 0 13
14 2921.44 0 0 1 0 0 0 0 0 0 0 0 0 14
15 2981.85 0 0 0 1 0 0 0 0 0 0 0 0 15
16 3080.58 0 0 0 0 1 0 0 0 0 0 0 0 16
17 3106.22 0 0 0 0 0 1 0 0 0 0 0 0 17
18 3119.31 0 0 0 0 0 0 1 0 0 0 0 0 18
19 3061.26 0 0 0 0 0 0 0 1 0 0 0 0 19
20 3097.31 0 0 0 0 0 0 0 0 1 0 0 0 20
21 3161.69 0 0 0 0 0 0 0 0 0 1 0 0 21
22 3257.16 0 0 0 0 0 0 0 0 0 0 1 0 22
23 3277.01 0 0 0 0 0 0 0 0 0 0 0 1 23
24 3295.32 0 0 0 0 0 0 0 0 0 0 0 0 24
25 3363.99 0 1 0 0 0 0 0 0 0 0 0 0 25
26 3494.17 0 0 1 0 0 0 0 0 0 0 0 0 26
27 3667.03 0 0 0 1 0 0 0 0 0 0 0 0 27
28 3813.06 0 0 0 0 1 0 0 0 0 0 0 0 28
29 3917.96 0 0 0 0 0 1 0 0 0 0 0 0 29
30 3895.51 0 0 0 0 0 0 1 0 0 0 0 0 30
31 3801.06 0 0 0 0 0 0 0 1 0 0 0 0 31
32 3570.12 0 0 0 0 0 0 0 0 1 0 0 0 32
33 3701.61 0 0 0 0 0 0 0 0 0 1 0 0 33
34 3862.27 0 0 0 0 0 0 0 0 0 0 1 0 34
35 3970.10 0 0 0 0 0 0 0 0 0 0 0 1 35
36 4138.52 0 0 0 0 0 0 0 0 0 0 0 0 36
37 4199.75 0 1 0 0 0 0 0 0 0 0 0 0 37
38 4290.89 0 0 1 0 0 0 0 0 0 0 0 0 38
39 4443.91 0 0 0 1 0 0 0 0 0 0 0 0 39
40 4502.64 0 0 0 0 1 0 0 0 0 0 0 0 40
41 4356.98 0 0 0 0 0 1 0 0 0 0 0 0 41
42 4591.27 0 0 0 0 0 0 1 0 0 0 0 0 42
43 4696.96 0 0 0 0 0 0 0 1 0 0 0 0 43
44 4621.40 0 0 0 0 0 0 0 0 1 0 0 0 44
45 4562.84 0 0 0 0 0 0 0 0 0 1 0 0 45
46 4202.52 0 0 0 0 0 0 0 0 0 0 1 0 46
47 4296.49 0 0 0 0 0 0 0 0 0 0 0 1 47
48 4435.23 0 0 0 0 0 0 0 0 0 0 0 0 48
49 4105.18 0 1 0 0 0 0 0 0 0 0 0 0 49
50 4116.68 0 0 1 0 0 0 0 0 0 0 0 0 50
51 3844.49 0 0 0 1 0 0 0 0 0 0 0 0 51
52 3720.98 0 0 0 0 1 0 0 0 0 0 0 0 52
53 3674.40 0 0 0 0 0 1 0 0 0 0 0 0 53
54 3857.62 0 0 0 0 0 0 1 0 0 0 0 0 54
55 3801.06 0 0 0 0 0 0 0 1 0 0 0 0 55
56 3504.37 1 0 0 0 0 0 0 0 1 0 0 0 56
57 3032.60 1 0 0 0 0 0 0 0 0 1 0 0 57
58 3047.03 1 0 0 0 0 0 0 0 0 0 1 0 58
59 2962.34 1 0 0 0 0 0 0 0 0 0 0 1 59
60 2197.82 1 0 0 0 0 0 0 0 0 0 0 0 60
61 2014.45 1 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) dummy M1 M2 M3 M4
2288.3279 -1960.8993 -108.2569 55.0096 67.7796 80.9437
M5 M6 M7 M8 M9 M10
21.4877 75.1578 0.6478 248.0098 139.0098 90.2219
M11 t
106.3799 40.7939
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-797.495 -166.933 6.987 247.472 653.845
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2288.3279 208.0743 10.998 1.36e-14 ***
dummy -1960.8993 209.6422 -9.354 2.65e-12 ***
M1 -108.2569 241.5752 -0.448 0.656
M2 55.0096 254.7634 0.216 0.830
M3 67.7796 254.6138 0.266 0.791
M4 80.9437 254.5096 0.318 0.752
M5 21.4877 254.4508 0.084 0.933
M6 75.1578 254.4375 0.295 0.769
M7 0.6478 254.4696 0.003 0.998
M8 248.0098 252.2118 0.983 0.330
M9 139.0098 252.0512 0.552 0.584
M10 90.2219 251.9364 0.358 0.722
M11 106.3799 251.8675 0.422 0.675
t 40.7939 3.4015 11.993 6.62e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 398.2 on 47 degrees of freedom
Multiple R-squared: 0.7774, Adjusted R-squared: 0.7159
F-statistic: 12.63 on 13 and 47 DF, p-value: 3.014e-11
> 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,] 5.001429e-04 1.000286e-03 0.9994999
[2,] 3.500445e-05 7.000890e-05 0.9999650
[3,] 2.104210e-06 4.208420e-06 0.9999979
[4,] 1.629705e-07 3.259410e-07 0.9999998
[5,] 2.383981e-08 4.767962e-08 1.0000000
[6,] 1.386600e-08 2.773201e-08 1.0000000
[7,] 1.948392e-09 3.896783e-09 1.0000000
[8,] 3.733205e-09 7.466410e-09 1.0000000
[9,] 4.634161e-09 9.268321e-09 1.0000000
[10,] 7.150600e-10 1.430120e-09 1.0000000
[11,] 1.019518e-10 2.039036e-10 1.0000000
[12,] 2.668990e-11 5.337979e-11 1.0000000
[13,] 7.513105e-11 1.502621e-10 1.0000000
[14,] 2.322268e-11 4.644535e-11 1.0000000
[15,] 5.306756e-12 1.061351e-11 1.0000000
[16,] 3.725288e-10 7.450576e-10 1.0000000
[17,] 1.799112e-09 3.598225e-09 1.0000000
[18,] 4.120032e-09 8.240063e-09 1.0000000
[19,] 4.019451e-08 8.038902e-08 1.0000000
[20,] 1.273883e-07 2.547766e-07 0.9999999
[21,] 7.369412e-08 1.473882e-07 0.9999999
[22,] 1.772684e-07 3.545369e-07 0.9999998
[23,] 8.963119e-08 1.792624e-07 0.9999999
[24,] 2.183182e-08 4.366363e-08 1.0000000
[25,] 1.613571e-08 3.227141e-08 1.0000000
[26,] 5.103399e-09 1.020680e-08 1.0000000
[27,] 1.855196e-08 3.710393e-08 1.0000000
[28,] 2.176731e-07 4.353462e-07 0.9999998
> postscript(file="/var/www/html/rcomp/tmp/1pohj1229721805.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/2je411229721805.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/37eyc1229721805.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/4qbek1229721805.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/5ed1i1229721805.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
-46.304957 -228.205391 -128.049391 -92.197391 -105.145391 -135.439391
7 8 9 10 11 12
-166.933391 -408.069252 -346.433252 -288.649252 -197.801252 -21.095252
13 14 15 16 17 18
138.877739 6.987304 13.833304 58.605304 102.907304 21.533304
19 20 21 22 23 24
-2.800696 -254.906557 -122.320557 -18.856557 -55.958557 27.937443
25 26 27 28 29 30
164.070435 90.190000 209.486000 301.558000 425.120000 308.206000
31 32 33 34 35 36
247.472000 -271.623861 -71.927861 96.726139 147.604139 381.610139
37 38 39 40 41 42
510.303130 397.382696 496.838696 501.610696 374.612696 514.438696
43 44 45 46 47 48
653.844696 290.128835 299.774835 -52.551165 -15.533165 188.792835
49 50 51 52 53 54
-73.794174 -266.354609 -592.108609 -769.576609 -797.494609 -708.738609
55 56 57 58 59 60
-731.582609 644.470835 240.906835 263.330835 121.688835 -577.245165
61
-693.152174
> postscript(file="/var/www/html/rcomp/tmp/675mo1229721805.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 -46.304957 NA
1 -228.205391 -46.304957
2 -128.049391 -228.205391
3 -92.197391 -128.049391
4 -105.145391 -92.197391
5 -135.439391 -105.145391
6 -166.933391 -135.439391
7 -408.069252 -166.933391
8 -346.433252 -408.069252
9 -288.649252 -346.433252
10 -197.801252 -288.649252
11 -21.095252 -197.801252
12 138.877739 -21.095252
13 6.987304 138.877739
14 13.833304 6.987304
15 58.605304 13.833304
16 102.907304 58.605304
17 21.533304 102.907304
18 -2.800696 21.533304
19 -254.906557 -2.800696
20 -122.320557 -254.906557
21 -18.856557 -122.320557
22 -55.958557 -18.856557
23 27.937443 -55.958557
24 164.070435 27.937443
25 90.190000 164.070435
26 209.486000 90.190000
27 301.558000 209.486000
28 425.120000 301.558000
29 308.206000 425.120000
30 247.472000 308.206000
31 -271.623861 247.472000
32 -71.927861 -271.623861
33 96.726139 -71.927861
34 147.604139 96.726139
35 381.610139 147.604139
36 510.303130 381.610139
37 397.382696 510.303130
38 496.838696 397.382696
39 501.610696 496.838696
40 374.612696 501.610696
41 514.438696 374.612696
42 653.844696 514.438696
43 290.128835 653.844696
44 299.774835 290.128835
45 -52.551165 299.774835
46 -15.533165 -52.551165
47 188.792835 -15.533165
48 -73.794174 188.792835
49 -266.354609 -73.794174
50 -592.108609 -266.354609
51 -769.576609 -592.108609
52 -797.494609 -769.576609
53 -708.738609 -797.494609
54 -731.582609 -708.738609
55 644.470835 -731.582609
56 240.906835 644.470835
57 263.330835 240.906835
58 121.688835 263.330835
59 -577.245165 121.688835
60 -693.152174 -577.245165
61 NA -693.152174
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -228.205391 -46.304957
[2,] -128.049391 -228.205391
[3,] -92.197391 -128.049391
[4,] -105.145391 -92.197391
[5,] -135.439391 -105.145391
[6,] -166.933391 -135.439391
[7,] -408.069252 -166.933391
[8,] -346.433252 -408.069252
[9,] -288.649252 -346.433252
[10,] -197.801252 -288.649252
[11,] -21.095252 -197.801252
[12,] 138.877739 -21.095252
[13,] 6.987304 138.877739
[14,] 13.833304 6.987304
[15,] 58.605304 13.833304
[16,] 102.907304 58.605304
[17,] 21.533304 102.907304
[18,] -2.800696 21.533304
[19,] -254.906557 -2.800696
[20,] -122.320557 -254.906557
[21,] -18.856557 -122.320557
[22,] -55.958557 -18.856557
[23,] 27.937443 -55.958557
[24,] 164.070435 27.937443
[25,] 90.190000 164.070435
[26,] 209.486000 90.190000
[27,] 301.558000 209.486000
[28,] 425.120000 301.558000
[29,] 308.206000 425.120000
[30,] 247.472000 308.206000
[31,] -271.623861 247.472000
[32,] -71.927861 -271.623861
[33,] 96.726139 -71.927861
[34,] 147.604139 96.726139
[35,] 381.610139 147.604139
[36,] 510.303130 381.610139
[37,] 397.382696 510.303130
[38,] 496.838696 397.382696
[39,] 501.610696 496.838696
[40,] 374.612696 501.610696
[41,] 514.438696 374.612696
[42,] 653.844696 514.438696
[43,] 290.128835 653.844696
[44,] 299.774835 290.128835
[45,] -52.551165 299.774835
[46,] -15.533165 -52.551165
[47,] 188.792835 -15.533165
[48,] -73.794174 188.792835
[49,] -266.354609 -73.794174
[50,] -592.108609 -266.354609
[51,] -769.576609 -592.108609
[52,] -797.494609 -769.576609
[53,] -708.738609 -797.494609
[54,] -731.582609 -708.738609
[55,] 644.470835 -731.582609
[56,] 240.906835 644.470835
[57,] 263.330835 240.906835
[58,] 121.688835 263.330835
[59,] -577.245165 121.688835
[60,] -693.152174 -577.245165
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -228.205391 -46.304957
2 -128.049391 -228.205391
3 -92.197391 -128.049391
4 -105.145391 -92.197391
5 -135.439391 -105.145391
6 -166.933391 -135.439391
7 -408.069252 -166.933391
8 -346.433252 -408.069252
9 -288.649252 -346.433252
10 -197.801252 -288.649252
11 -21.095252 -197.801252
12 138.877739 -21.095252
13 6.987304 138.877739
14 13.833304 6.987304
15 58.605304 13.833304
16 102.907304 58.605304
17 21.533304 102.907304
18 -2.800696 21.533304
19 -254.906557 -2.800696
20 -122.320557 -254.906557
21 -18.856557 -122.320557
22 -55.958557 -18.856557
23 27.937443 -55.958557
24 164.070435 27.937443
25 90.190000 164.070435
26 209.486000 90.190000
27 301.558000 209.486000
28 425.120000 301.558000
29 308.206000 425.120000
30 247.472000 308.206000
31 -271.623861 247.472000
32 -71.927861 -271.623861
33 96.726139 -71.927861
34 147.604139 96.726139
35 381.610139 147.604139
36 510.303130 381.610139
37 397.382696 510.303130
38 496.838696 397.382696
39 501.610696 496.838696
40 374.612696 501.610696
41 514.438696 374.612696
42 653.844696 514.438696
43 290.128835 653.844696
44 299.774835 290.128835
45 -52.551165 299.774835
46 -15.533165 -52.551165
47 188.792835 -15.533165
48 -73.794174 188.792835
49 -266.354609 -73.794174
50 -592.108609 -266.354609
51 -769.576609 -592.108609
52 -797.494609 -769.576609
53 -708.738609 -797.494609
54 -731.582609 -708.738609
55 644.470835 -731.582609
56 240.906835 644.470835
57 263.330835 240.906835
58 121.688835 263.330835
59 -577.245165 121.688835
60 -693.152174 -577.245165
> 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/791i81229721805.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/8tskm1229721805.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/9xjut1229721805.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/10p5ku1229721805.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/11t65v1229721805.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/12x1wp1229721805.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/13iv3b1229721805.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/14wh2c1229721805.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/15c0oy1229721805.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/16k9tf1229721806.tab")
+ }
>
> system("convert tmp/1pohj1229721805.ps tmp/1pohj1229721805.png")
> system("convert tmp/2je411229721805.ps tmp/2je411229721805.png")
> system("convert tmp/37eyc1229721805.ps tmp/37eyc1229721805.png")
> system("convert tmp/4qbek1229721805.ps tmp/4qbek1229721805.png")
> system("convert tmp/5ed1i1229721805.ps tmp/5ed1i1229721805.png")
> system("convert tmp/675mo1229721805.ps tmp/675mo1229721805.png")
> system("convert tmp/791i81229721805.ps tmp/791i81229721805.png")
> system("convert tmp/8tskm1229721805.ps tmp/8tskm1229721805.png")
> system("convert tmp/9xjut1229721805.ps tmp/9xjut1229721805.png")
> system("convert tmp/10p5ku1229721805.ps tmp/10p5ku1229721805.png")
>
>
> proc.time()
user system elapsed
2.401 1.544 2.847