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(2120.88,0,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,1,4356.98,1,4591.27,1,4696.96,1,4621.4,1,4562.84,1,4202.52,1,4296.49,1,4435.23,1,4105.18,1,4116.68,1,3844.49,1,3720.98,1,3674.4,1,3857.62,1,3801.06,1,3504.37,1,3032.6,1,3047.03,1,2962.34,1,2197.82,1),dim=c(2,61),dimnames=list(c('Bel20','dummy'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Bel20','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
Bel20 dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 2120.88 0 1 0 0 0 0 0 0 0 0 0 0 1
2 2174.56 0 0 1 0 0 0 0 0 0 0 0 0 2
3 2196.72 0 0 0 1 0 0 0 0 0 0 0 0 3
4 2350.44 0 0 0 0 1 0 0 0 0 0 0 0 4
5 2440.25 0 0 0 0 0 1 0 0 0 0 0 0 5
6 2408.64 0 0 0 0 0 0 1 0 0 0 0 0 6
7 2472.81 0 0 0 0 0 0 0 1 0 0 0 0 7
8 2407.60 0 0 0 0 0 0 0 0 1 0 0 0 8
9 2454.62 0 0 0 0 0 0 0 0 0 1 0 0 9
10 2448.05 0 0 0 0 0 0 0 0 0 0 1 0 10
11 2497.84 0 0 0 0 0 0 0 0 0 0 0 1 11
12 2645.64 0 0 0 0 0 0 0 0 0 0 0 0 12
13 2756.76 0 1 0 0 0 0 0 0 0 0 0 0 13
14 2849.27 0 0 1 0 0 0 0 0 0 0 0 0 14
15 2921.44 0 0 0 1 0 0 0 0 0 0 0 0 15
16 2981.85 0 0 0 0 1 0 0 0 0 0 0 0 16
17 3080.58 0 0 0 0 0 1 0 0 0 0 0 0 17
18 3106.22 0 0 0 0 0 0 1 0 0 0 0 0 18
19 3119.31 0 0 0 0 0 0 0 1 0 0 0 0 19
20 3061.26 0 0 0 0 0 0 0 0 1 0 0 0 20
21 3097.31 0 0 0 0 0 0 0 0 0 1 0 0 21
22 3161.69 0 0 0 0 0 0 0 0 0 0 1 0 22
23 3257.16 0 0 0 0 0 0 0 0 0 0 0 1 23
24 3277.01 0 0 0 0 0 0 0 0 0 0 0 0 24
25 3295.32 0 1 0 0 0 0 0 0 0 0 0 0 25
26 3363.99 0 0 1 0 0 0 0 0 0 0 0 0 26
27 3494.17 0 0 0 1 0 0 0 0 0 0 0 0 27
28 3667.03 0 0 0 0 1 0 0 0 0 0 0 0 28
29 3813.06 0 0 0 0 0 1 0 0 0 0 0 0 29
30 3917.96 0 0 0 0 0 0 1 0 0 0 0 0 30
31 3895.51 0 0 0 0 0 0 0 1 0 0 0 0 31
32 3801.06 0 0 0 0 0 0 0 0 1 0 0 0 32
33 3570.12 0 0 0 0 0 0 0 0 0 1 0 0 33
34 3701.61 0 0 0 0 0 0 0 0 0 0 1 0 34
35 3862.27 0 0 0 0 0 0 0 0 0 0 0 1 35
36 3970.10 0 0 0 0 0 0 0 0 0 0 0 0 36
37 4138.52 0 1 0 0 0 0 0 0 0 0 0 0 37
38 4199.75 0 0 1 0 0 0 0 0 0 0 0 0 38
39 4290.89 0 0 0 1 0 0 0 0 0 0 0 0 39
40 4443.91 0 0 0 0 1 0 0 0 0 0 0 0 40
41 4502.64 1 0 0 0 0 1 0 0 0 0 0 0 41
42 4356.98 1 0 0 0 0 0 1 0 0 0 0 0 42
43 4591.27 1 0 0 0 0 0 0 1 0 0 0 0 43
44 4696.96 1 0 0 0 0 0 0 0 1 0 0 0 44
45 4621.40 1 0 0 0 0 0 0 0 0 1 0 0 45
46 4562.84 1 0 0 0 0 0 0 0 0 0 1 0 46
47 4202.52 1 0 0 0 0 0 0 0 0 0 0 1 47
48 4296.49 1 0 0 0 0 0 0 0 0 0 0 0 48
49 4435.23 1 1 0 0 0 0 0 0 0 0 0 0 49
50 4105.18 1 0 1 0 0 0 0 0 0 0 0 0 50
51 4116.68 1 0 0 1 0 0 0 0 0 0 0 0 51
52 3844.49 1 0 0 0 1 0 0 0 0 0 0 0 52
53 3720.98 1 0 0 0 0 1 0 0 0 0 0 0 53
54 3674.40 1 0 0 0 0 0 1 0 0 0 0 0 54
55 3857.62 1 0 0 0 0 0 0 1 0 0 0 0 55
56 3801.06 1 0 0 0 0 0 0 0 1 0 0 0 56
57 3504.37 1 0 0 0 0 0 0 0 0 1 0 0 57
58 3032.60 1 0 0 0 0 0 0 0 0 0 1 0 58
59 3047.03 1 0 0 0 0 0 0 0 0 0 0 1 59
60 2962.34 1 0 0 0 0 0 0 0 0 0 0 0 60
61 2197.82 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
2222.11 -413.26 -109.68 207.11 234.39 249.80
M5 M6 M7 M8 M9 M10
348.26 291.44 347.75 275.88 133.71 27.35
M11 t
-18.80 38.15
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1828.69 -274.08 -34.31 259.82 971.60
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2222.108 326.719 6.801 1.64e-08 ***
dummy -413.262 284.966 -1.450 0.154
M1 -109.680 357.315 -0.307 0.760
M2 207.113 374.857 0.553 0.583
M3 234.390 374.325 0.626 0.534
M4 249.801 373.950 0.668 0.507
M5 348.258 376.139 0.926 0.359
M6 291.443 375.119 0.777 0.441
M7 347.754 374.253 0.929 0.358
M8 275.885 373.543 0.739 0.464
M9 133.707 372.990 0.358 0.722
M10 27.348 372.595 0.073 0.942
M11 -18.799 372.357 -0.050 0.960
t 38.153 7.679 4.968 9.36e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 588.6 on 47 degrees of freedom
Multiple R-squared: 0.5095, Adjusted R-squared: 0.3738
F-statistic: 3.755 on 13 and 47 DF, p-value: 0.0004043
> 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,] 2.009216e-04 4.018432e-04 0.9997991
[2,] 1.249641e-05 2.499282e-05 0.9999875
[3,] 6.826859e-07 1.365372e-06 0.9999993
[4,] 3.665048e-08 7.330096e-08 1.0000000
[5,] 2.119960e-09 4.239920e-09 1.0000000
[6,] 2.963138e-10 5.926276e-10 1.0000000
[7,] 1.878890e-10 3.757779e-10 1.0000000
[8,] 2.264060e-11 4.528120e-11 1.0000000
[9,] 1.328272e-10 2.656543e-10 1.0000000
[10,] 2.277129e-10 4.554259e-10 1.0000000
[11,] 1.214690e-10 2.429381e-10 1.0000000
[12,] 1.301822e-10 2.603644e-10 1.0000000
[13,] 6.054182e-11 1.210836e-10 1.0000000
[14,] 2.872747e-10 5.745494e-10 1.0000000
[15,] 1.974015e-10 3.948031e-10 1.0000000
[16,] 1.647798e-10 3.295596e-10 1.0000000
[17,] 6.833034e-09 1.366607e-08 1.0000000
[18,] 1.303483e-08 2.606965e-08 1.0000000
[19,] 6.982033e-09 1.396407e-08 1.0000000
[20,] 4.646749e-09 9.293499e-09 1.0000000
[21,] 1.821819e-09 3.643637e-09 1.0000000
[22,] 7.002956e-10 1.400591e-09 1.0000000
[23,] 3.441592e-10 6.883183e-10 1.0000000
[24,] 1.192903e-10 2.385806e-10 1.0000000
[25,] 5.327086e-11 1.065417e-10 1.0000000
[26,] 2.352234e-10 4.704468e-10 1.0000000
[27,] 1.729939e-09 3.459879e-09 1.0000000
[28,] 8.165316e-08 1.633063e-07 0.9999999
> postscript(file="/var/www/html/rcomp/tmp/1pzkb1227826172.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/2jv0d1227826172.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/35t8n1227826172.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/4z9kg1227826172.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/58wa11227826172.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
-29.70191 -330.96736 -374.23736 -274.08136 -320.88186 -333.82986
7 8 9 10 11 12
-364.12386 -395.61786 -244.57386 -182.93786 -125.15386 -34.30586
13 14 15 16 17 18
148.34052 -114.09493 -107.35493 -100.50893 -138.38943 -94.08743
19 20 21 22 23 24
-175.46143 -199.79543 -59.72143 72.86457 176.32857 139.22657
25 26 27 28 29 30
229.06295 -57.21250 7.53750 126.83350 136.25300 259.81500
31 32 33 34 35 36
142.90100 82.16700 -44.74900 154.94700 323.60100 374.47900
37 38 39 40 41 42
614.42538 320.70993 346.41993 445.87593 781.25793 654.25993
43 44 45 46 47 48
794.08593 933.49193 961.95593 971.60193 619.27593 656.29393
49 50 51 52 53 54
866.56031 181.56486 127.63486 -198.11914 -458.23964 -486.15764
55 56 57 58 59 60
-397.40164 -420.24564 -612.91164 -1016.47564 -994.05164 -1135.69364
61
-1828.68725
> postscript(file="/var/www/html/rcomp/tmp/6ymwc1227826172.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 -29.70191 NA
1 -330.96736 -29.70191
2 -374.23736 -330.96736
3 -274.08136 -374.23736
4 -320.88186 -274.08136
5 -333.82986 -320.88186
6 -364.12386 -333.82986
7 -395.61786 -364.12386
8 -244.57386 -395.61786
9 -182.93786 -244.57386
10 -125.15386 -182.93786
11 -34.30586 -125.15386
12 148.34052 -34.30586
13 -114.09493 148.34052
14 -107.35493 -114.09493
15 -100.50893 -107.35493
16 -138.38943 -100.50893
17 -94.08743 -138.38943
18 -175.46143 -94.08743
19 -199.79543 -175.46143
20 -59.72143 -199.79543
21 72.86457 -59.72143
22 176.32857 72.86457
23 139.22657 176.32857
24 229.06295 139.22657
25 -57.21250 229.06295
26 7.53750 -57.21250
27 126.83350 7.53750
28 136.25300 126.83350
29 259.81500 136.25300
30 142.90100 259.81500
31 82.16700 142.90100
32 -44.74900 82.16700
33 154.94700 -44.74900
34 323.60100 154.94700
35 374.47900 323.60100
36 614.42538 374.47900
37 320.70993 614.42538
38 346.41993 320.70993
39 445.87593 346.41993
40 781.25793 445.87593
41 654.25993 781.25793
42 794.08593 654.25993
43 933.49193 794.08593
44 961.95593 933.49193
45 971.60193 961.95593
46 619.27593 971.60193
47 656.29393 619.27593
48 866.56031 656.29393
49 181.56486 866.56031
50 127.63486 181.56486
51 -198.11914 127.63486
52 -458.23964 -198.11914
53 -486.15764 -458.23964
54 -397.40164 -486.15764
55 -420.24564 -397.40164
56 -612.91164 -420.24564
57 -1016.47564 -612.91164
58 -994.05164 -1016.47564
59 -1135.69364 -994.05164
60 -1828.68725 -1135.69364
61 NA -1828.68725
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -330.96736 -29.70191
[2,] -374.23736 -330.96736
[3,] -274.08136 -374.23736
[4,] -320.88186 -274.08136
[5,] -333.82986 -320.88186
[6,] -364.12386 -333.82986
[7,] -395.61786 -364.12386
[8,] -244.57386 -395.61786
[9,] -182.93786 -244.57386
[10,] -125.15386 -182.93786
[11,] -34.30586 -125.15386
[12,] 148.34052 -34.30586
[13,] -114.09493 148.34052
[14,] -107.35493 -114.09493
[15,] -100.50893 -107.35493
[16,] -138.38943 -100.50893
[17,] -94.08743 -138.38943
[18,] -175.46143 -94.08743
[19,] -199.79543 -175.46143
[20,] -59.72143 -199.79543
[21,] 72.86457 -59.72143
[22,] 176.32857 72.86457
[23,] 139.22657 176.32857
[24,] 229.06295 139.22657
[25,] -57.21250 229.06295
[26,] 7.53750 -57.21250
[27,] 126.83350 7.53750
[28,] 136.25300 126.83350
[29,] 259.81500 136.25300
[30,] 142.90100 259.81500
[31,] 82.16700 142.90100
[32,] -44.74900 82.16700
[33,] 154.94700 -44.74900
[34,] 323.60100 154.94700
[35,] 374.47900 323.60100
[36,] 614.42538 374.47900
[37,] 320.70993 614.42538
[38,] 346.41993 320.70993
[39,] 445.87593 346.41993
[40,] 781.25793 445.87593
[41,] 654.25993 781.25793
[42,] 794.08593 654.25993
[43,] 933.49193 794.08593
[44,] 961.95593 933.49193
[45,] 971.60193 961.95593
[46,] 619.27593 971.60193
[47,] 656.29393 619.27593
[48,] 866.56031 656.29393
[49,] 181.56486 866.56031
[50,] 127.63486 181.56486
[51,] -198.11914 127.63486
[52,] -458.23964 -198.11914
[53,] -486.15764 -458.23964
[54,] -397.40164 -486.15764
[55,] -420.24564 -397.40164
[56,] -612.91164 -420.24564
[57,] -1016.47564 -612.91164
[58,] -994.05164 -1016.47564
[59,] -1135.69364 -994.05164
[60,] -1828.68725 -1135.69364
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -330.96736 -29.70191
2 -374.23736 -330.96736
3 -274.08136 -374.23736
4 -320.88186 -274.08136
5 -333.82986 -320.88186
6 -364.12386 -333.82986
7 -395.61786 -364.12386
8 -244.57386 -395.61786
9 -182.93786 -244.57386
10 -125.15386 -182.93786
11 -34.30586 -125.15386
12 148.34052 -34.30586
13 -114.09493 148.34052
14 -107.35493 -114.09493
15 -100.50893 -107.35493
16 -138.38943 -100.50893
17 -94.08743 -138.38943
18 -175.46143 -94.08743
19 -199.79543 -175.46143
20 -59.72143 -199.79543
21 72.86457 -59.72143
22 176.32857 72.86457
23 139.22657 176.32857
24 229.06295 139.22657
25 -57.21250 229.06295
26 7.53750 -57.21250
27 126.83350 7.53750
28 136.25300 126.83350
29 259.81500 136.25300
30 142.90100 259.81500
31 82.16700 142.90100
32 -44.74900 82.16700
33 154.94700 -44.74900
34 323.60100 154.94700
35 374.47900 323.60100
36 614.42538 374.47900
37 320.70993 614.42538
38 346.41993 320.70993
39 445.87593 346.41993
40 781.25793 445.87593
41 654.25993 781.25793
42 794.08593 654.25993
43 933.49193 794.08593
44 961.95593 933.49193
45 971.60193 961.95593
46 619.27593 971.60193
47 656.29393 619.27593
48 866.56031 656.29393
49 181.56486 866.56031
50 127.63486 181.56486
51 -198.11914 127.63486
52 -458.23964 -198.11914
53 -486.15764 -458.23964
54 -397.40164 -486.15764
55 -420.24564 -397.40164
56 -612.91164 -420.24564
57 -1016.47564 -612.91164
58 -994.05164 -1016.47564
59 -1135.69364 -994.05164
60 -1828.68725 -1135.69364
> 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/73kjq1227826172.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/8mftu1227826172.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/9ac5z1227826172.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/10qmdj1227826172.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/11thms1227826172.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/12c5ue1227826172.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/132t9n1227826173.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/14z3nh1227826173.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/151y0v1227826173.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/1657q91227826173.tab")
+ }
>
> system("convert tmp/1pzkb1227826172.ps tmp/1pzkb1227826172.png")
> system("convert tmp/2jv0d1227826172.ps tmp/2jv0d1227826172.png")
> system("convert tmp/35t8n1227826172.ps tmp/35t8n1227826172.png")
> system("convert tmp/4z9kg1227826172.ps tmp/4z9kg1227826172.png")
> system("convert tmp/58wa11227826172.ps tmp/58wa11227826172.png")
> system("convert tmp/6ymwc1227826172.ps tmp/6ymwc1227826172.png")
> system("convert tmp/73kjq1227826172.ps tmp/73kjq1227826172.png")
> system("convert tmp/8mftu1227826172.ps tmp/8mftu1227826172.png")
> system("convert tmp/9ac5z1227826172.ps tmp/9ac5z1227826172.png")
> system("convert tmp/10qmdj1227826172.ps tmp/10qmdj1227826172.png")
>
>
> proc.time()
user system elapsed
2.354 1.542 2.796