R version 2.11.1 (2010-05-31)
Copyright (C) 2010 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(43071,990633,45552,1047696,36329,835567,37703,867169,50519,1161937,36798,846354,37056,852288,44927,1033321,37635,865605,62924,1447252,8170,187910,27438,631074,27429,630867,33666,774318,27733,637859,33228,764244,25699,591077,303936,6990528,30169,693887,35117,807691,34870,802010,56676,1303548,7054,162242,29722,683606,41629,957467,41117,945691,39341,904843,39486,908178,48138,1107174,45633,1049559,41756,960388,47221,1086083,50530,1162190,68184,1568232,8771,201733,37898,871654,41888,963424,40439,930097,40898,940654,38401,883223,52073,1197679,41547,955581,38529,886167,51321,1180383,41519,954937,69116,1589668,12657,291111,34801,800423,37967,873241,39401,906223,33425,768775,36222,833106,48428,1113844,40891,940493,36432,837936,50669,1165387,39556,909788,68906,1584838),dim=c(2,58),dimnames=list(c('Verkoopcijfers','Totaleuitstootkm/u'),1:58))
> y <- array(NA,dim=c(2,58),dimnames=list(c('Verkoopcijfers','Totaleuitstootkm/u'),1:58))
> 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
Verkoopcijfers Totaleuitstootkm/u
1 43071 990633
2 45552 1047696
3 36329 835567
4 37703 867169
5 50519 1161937
6 36798 846354
7 37056 852288
8 44927 1033321
9 37635 865605
10 62924 1447252
11 8170 187910
12 27438 631074
13 27429 630867
14 33666 774318
15 27733 637859
16 33228 764244
17 25699 591077
18 303936 6990528
19 30169 693887
20 35117 807691
21 34870 802010
22 56676 1303548
23 7054 162242
24 29722 683606
25 41629 957467
26 41117 945691
27 39341 904843
28 39486 908178
29 48138 1107174
30 45633 1049559
31 41756 960388
32 47221 1086083
33 50530 1162190
34 68184 1568232
35 8771 201733
36 37898 871654
37 41888 963424
38 40439 930097
39 40898 940654
40 38401 883223
41 52073 1197679
42 41547 955581
43 38529 886167
44 51321 1180383
45 41519 954937
46 69116 1589668
47 12657 291111
48 34801 800423
49 37967 873241
50 39401 906223
51 33425 768775
52 36222 833106
53 48428 1113844
54 40891 940493
55 36432 837936
56 50669 1165387
57 39556 909788
58 68906 1584838
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `Totaleuitstootkm/u`
-2.161e-12 4.348e-02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.551e-12 -1.258e-12 -5.474e-13 -2.253e-13 3.047e-11
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.161e-12 1.086e-12 -1.989e+00 0.0516 .
`Totaleuitstootkm/u` 4.348e-02 8.236e-19 5.279e+16 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.286e-12 on 56 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 2.787e+33 on 1 and 56 DF, p-value: < 2.2e-16
> 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,] 6.436248e-01 7.127504e-01 3.563752e-01
[2,] 3.389929e-02 6.779859e-02 9.661007e-01
[3,] 1.753361e-05 3.506722e-05 9.999825e-01
[4,] 1.884253e-03 3.768506e-03 9.981157e-01
[5,] 1.093500e-09 2.187001e-09 1.000000e+00
[6,] 1.089473e-01 2.178945e-01 8.910527e-01
[7,] 1.319261e-02 2.638522e-02 9.868074e-01
[8,] 6.402307e-01 7.195386e-01 3.597693e-01
[9,] 1.656231e-02 3.312462e-02 9.834377e-01
[10,] 1.557030e-04 3.114060e-04 9.998443e-01
[11,] 1.562193e-08 3.124385e-08 1.000000e+00
[12,] 3.597121e-02 7.194242e-02 9.640288e-01
[13,] 4.076187e-07 8.152373e-07 9.999996e-01
[14,] 1.000000e+00 1.498157e-13 7.490783e-14
[15,] 1.000000e+00 2.701707e-09 1.350853e-09
[16,] 9.999998e-01 4.628501e-07 2.314250e-07
[17,] 9.998403e-01 3.194470e-04 1.597235e-04
[18,] 9.999288e-01 1.423513e-04 7.117564e-05
[19,] 9.999228e-01 1.544861e-04 7.724303e-05
[20,] 9.999779e-01 4.420771e-05 2.210385e-05
[21,] 9.988776e-01 2.244755e-03 1.122378e-03
[22,] 1.000000e+00 9.455790e-10 4.727895e-10
[23,] 1.373536e-02 2.747072e-02 9.862646e-01
[24,] 9.825921e-01 3.481581e-02 1.740790e-02
[25,] 9.999400e-01 1.200002e-04 6.000008e-05
[26,] 1.000000e+00 5.799741e-09 2.899871e-09
[27,] 1.000000e+00 1.044363e-10 5.221813e-11
[28,] 1.000000e+00 3.348098e-10 1.674049e-10
[29,] 9.999981e-01 3.714643e-06 1.857322e-06
[30,] 9.218618e-06 1.843724e-05 9.999908e-01
[31,] 9.999988e-01 2.405895e-06 1.202947e-06
[32,] 1.000000e+00 1.731789e-14 8.658947e-15
[33,] 9.999976e-01 4.827371e-06 2.413686e-06
[34,] 9.999996e-01 8.964540e-07 4.482270e-07
[35,] 1.000000e+00 7.869571e-16 3.934786e-16
[36,] 9.999838e-01 3.249696e-05 1.624848e-05
[37,] 4.492982e-01 8.985965e-01 5.507018e-01
[38,] 9.984214e-01 3.157149e-03 1.578574e-03
[39,] 3.503825e-03 7.007651e-03 9.964962e-01
[40,] 1.940355e-04 3.880711e-04 9.998060e-01
[41,] 2.189865e-01 4.379729e-01 7.810135e-01
[42,] 9.999596e-01 8.089711e-05 4.044855e-05
[43,] 9.999963e-01 7.315298e-06 3.657649e-06
[44,] 9.999497e-01 1.006929e-04 5.034644e-05
[45,] 9.999900e-01 1.999466e-05 9.997329e-06
[46,] 9.925685e-01 1.486308e-02 7.431540e-03
[47,] 9.970320e-01 5.936030e-03 2.968015e-03
[48,] 9.855696e-01 2.886089e-02 1.443044e-02
[49,] 6.562836e-02 1.312567e-01 9.343716e-01
> postscript(file="/var/www/rcomp/tmp/1j2yj1290764213.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/rcomp/tmp/2utf41290764213.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/rcomp/tmp/3utf41290764213.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/rcomp/tmp/4utf41290764213.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/rcomp/tmp/5utf41290764213.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 = 58
Frequency = 1
1 2 3 4 5
2.223660e-11 3.046721e-11 4.264234e-14 -2.133852e-13 4.864978e-13
6 7 8 9 10
-3.006591e-13 -5.388082e-13 -3.454553e-13 -5.558745e-13 -3.584014e-12
11 12 13 14 15
-3.354957e-13 -4.437174e-12 -3.256791e-12 -1.437033e-12 -2.754957e-12
16 17 18 19 20
-8.526191e-13 -3.370098e-12 -1.840402e-12 -1.763246e-12 3.742773e-13
21 22 23 24 25
-6.130548e-13 -5.412007e-13 -1.284844e-12 -1.936455e-12 -2.608717e-13
26 27 28 29 30
-8.459932e-13 -8.486301e-13 -3.750319e-13 -3.416514e-13 -2.001313e-13
31 32 33 34 35
-3.332936e-13 -1.474550e-13 2.286869e-13 -2.951145e-12 -1.835410e-12
36 37 38 39 40
5.646005e-14 -2.806516e-13 -9.289156e-13 -5.535642e-13 -1.083711e-12
41 42 43 44 45
8.523198e-13 -2.752380e-13 -1.179348e-12 1.051590e-12 -5.871067e-13
46 47 48 49 50
-1.560006e-12 -5.550830e-12 3.556371e-14 -5.926025e-13 -3.190862e-13
51 52 53 54 55
-6.287603e-13 -3.539048e-13 1.219320e-13 -8.735600e-13 -4.351756e-13
56 57 58
8.433534e-13 -1.047532e-12 -2.445962e-12
> postscript(file="/var/www/rcomp/tmp/6nlwp1290764213.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 = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 2.223660e-11 NA
1 3.046721e-11 2.223660e-11
2 4.264234e-14 3.046721e-11
3 -2.133852e-13 4.264234e-14
4 4.864978e-13 -2.133852e-13
5 -3.006591e-13 4.864978e-13
6 -5.388082e-13 -3.006591e-13
7 -3.454553e-13 -5.388082e-13
8 -5.558745e-13 -3.454553e-13
9 -3.584014e-12 -5.558745e-13
10 -3.354957e-13 -3.584014e-12
11 -4.437174e-12 -3.354957e-13
12 -3.256791e-12 -4.437174e-12
13 -1.437033e-12 -3.256791e-12
14 -2.754957e-12 -1.437033e-12
15 -8.526191e-13 -2.754957e-12
16 -3.370098e-12 -8.526191e-13
17 -1.840402e-12 -3.370098e-12
18 -1.763246e-12 -1.840402e-12
19 3.742773e-13 -1.763246e-12
20 -6.130548e-13 3.742773e-13
21 -5.412007e-13 -6.130548e-13
22 -1.284844e-12 -5.412007e-13
23 -1.936455e-12 -1.284844e-12
24 -2.608717e-13 -1.936455e-12
25 -8.459932e-13 -2.608717e-13
26 -8.486301e-13 -8.459932e-13
27 -3.750319e-13 -8.486301e-13
28 -3.416514e-13 -3.750319e-13
29 -2.001313e-13 -3.416514e-13
30 -3.332936e-13 -2.001313e-13
31 -1.474550e-13 -3.332936e-13
32 2.286869e-13 -1.474550e-13
33 -2.951145e-12 2.286869e-13
34 -1.835410e-12 -2.951145e-12
35 5.646005e-14 -1.835410e-12
36 -2.806516e-13 5.646005e-14
37 -9.289156e-13 -2.806516e-13
38 -5.535642e-13 -9.289156e-13
39 -1.083711e-12 -5.535642e-13
40 8.523198e-13 -1.083711e-12
41 -2.752380e-13 8.523198e-13
42 -1.179348e-12 -2.752380e-13
43 1.051590e-12 -1.179348e-12
44 -5.871067e-13 1.051590e-12
45 -1.560006e-12 -5.871067e-13
46 -5.550830e-12 -1.560006e-12
47 3.556371e-14 -5.550830e-12
48 -5.926025e-13 3.556371e-14
49 -3.190862e-13 -5.926025e-13
50 -6.287603e-13 -3.190862e-13
51 -3.539048e-13 -6.287603e-13
52 1.219320e-13 -3.539048e-13
53 -8.735600e-13 1.219320e-13
54 -4.351756e-13 -8.735600e-13
55 8.433534e-13 -4.351756e-13
56 -1.047532e-12 8.433534e-13
57 -2.445962e-12 -1.047532e-12
58 NA -2.445962e-12
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.046721e-11 2.223660e-11
[2,] 4.264234e-14 3.046721e-11
[3,] -2.133852e-13 4.264234e-14
[4,] 4.864978e-13 -2.133852e-13
[5,] -3.006591e-13 4.864978e-13
[6,] -5.388082e-13 -3.006591e-13
[7,] -3.454553e-13 -5.388082e-13
[8,] -5.558745e-13 -3.454553e-13
[9,] -3.584014e-12 -5.558745e-13
[10,] -3.354957e-13 -3.584014e-12
[11,] -4.437174e-12 -3.354957e-13
[12,] -3.256791e-12 -4.437174e-12
[13,] -1.437033e-12 -3.256791e-12
[14,] -2.754957e-12 -1.437033e-12
[15,] -8.526191e-13 -2.754957e-12
[16,] -3.370098e-12 -8.526191e-13
[17,] -1.840402e-12 -3.370098e-12
[18,] -1.763246e-12 -1.840402e-12
[19,] 3.742773e-13 -1.763246e-12
[20,] -6.130548e-13 3.742773e-13
[21,] -5.412007e-13 -6.130548e-13
[22,] -1.284844e-12 -5.412007e-13
[23,] -1.936455e-12 -1.284844e-12
[24,] -2.608717e-13 -1.936455e-12
[25,] -8.459932e-13 -2.608717e-13
[26,] -8.486301e-13 -8.459932e-13
[27,] -3.750319e-13 -8.486301e-13
[28,] -3.416514e-13 -3.750319e-13
[29,] -2.001313e-13 -3.416514e-13
[30,] -3.332936e-13 -2.001313e-13
[31,] -1.474550e-13 -3.332936e-13
[32,] 2.286869e-13 -1.474550e-13
[33,] -2.951145e-12 2.286869e-13
[34,] -1.835410e-12 -2.951145e-12
[35,] 5.646005e-14 -1.835410e-12
[36,] -2.806516e-13 5.646005e-14
[37,] -9.289156e-13 -2.806516e-13
[38,] -5.535642e-13 -9.289156e-13
[39,] -1.083711e-12 -5.535642e-13
[40,] 8.523198e-13 -1.083711e-12
[41,] -2.752380e-13 8.523198e-13
[42,] -1.179348e-12 -2.752380e-13
[43,] 1.051590e-12 -1.179348e-12
[44,] -5.871067e-13 1.051590e-12
[45,] -1.560006e-12 -5.871067e-13
[46,] -5.550830e-12 -1.560006e-12
[47,] 3.556371e-14 -5.550830e-12
[48,] -5.926025e-13 3.556371e-14
[49,] -3.190862e-13 -5.926025e-13
[50,] -6.287603e-13 -3.190862e-13
[51,] -3.539048e-13 -6.287603e-13
[52,] 1.219320e-13 -3.539048e-13
[53,] -8.735600e-13 1.219320e-13
[54,] -4.351756e-13 -8.735600e-13
[55,] 8.433534e-13 -4.351756e-13
[56,] -1.047532e-12 8.433534e-13
[57,] -2.445962e-12 -1.047532e-12
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.046721e-11 2.223660e-11
2 4.264234e-14 3.046721e-11
3 -2.133852e-13 4.264234e-14
4 4.864978e-13 -2.133852e-13
5 -3.006591e-13 4.864978e-13
6 -5.388082e-13 -3.006591e-13
7 -3.454553e-13 -5.388082e-13
8 -5.558745e-13 -3.454553e-13
9 -3.584014e-12 -5.558745e-13
10 -3.354957e-13 -3.584014e-12
11 -4.437174e-12 -3.354957e-13
12 -3.256791e-12 -4.437174e-12
13 -1.437033e-12 -3.256791e-12
14 -2.754957e-12 -1.437033e-12
15 -8.526191e-13 -2.754957e-12
16 -3.370098e-12 -8.526191e-13
17 -1.840402e-12 -3.370098e-12
18 -1.763246e-12 -1.840402e-12
19 3.742773e-13 -1.763246e-12
20 -6.130548e-13 3.742773e-13
21 -5.412007e-13 -6.130548e-13
22 -1.284844e-12 -5.412007e-13
23 -1.936455e-12 -1.284844e-12
24 -2.608717e-13 -1.936455e-12
25 -8.459932e-13 -2.608717e-13
26 -8.486301e-13 -8.459932e-13
27 -3.750319e-13 -8.486301e-13
28 -3.416514e-13 -3.750319e-13
29 -2.001313e-13 -3.416514e-13
30 -3.332936e-13 -2.001313e-13
31 -1.474550e-13 -3.332936e-13
32 2.286869e-13 -1.474550e-13
33 -2.951145e-12 2.286869e-13
34 -1.835410e-12 -2.951145e-12
35 5.646005e-14 -1.835410e-12
36 -2.806516e-13 5.646005e-14
37 -9.289156e-13 -2.806516e-13
38 -5.535642e-13 -9.289156e-13
39 -1.083711e-12 -5.535642e-13
40 8.523198e-13 -1.083711e-12
41 -2.752380e-13 8.523198e-13
42 -1.179348e-12 -2.752380e-13
43 1.051590e-12 -1.179348e-12
44 -5.871067e-13 1.051590e-12
45 -1.560006e-12 -5.871067e-13
46 -5.550830e-12 -1.560006e-12
47 3.556371e-14 -5.550830e-12
48 -5.926025e-13 3.556371e-14
49 -3.190862e-13 -5.926025e-13
50 -6.287603e-13 -3.190862e-13
51 -3.539048e-13 -6.287603e-13
52 1.219320e-13 -3.539048e-13
53 -8.735600e-13 1.219320e-13
54 -4.351756e-13 -8.735600e-13
55 8.433534e-13 -4.351756e-13
56 -1.047532e-12 8.433534e-13
57 -2.445962e-12 -1.047532e-12
> 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/rcomp/tmp/7fcwa1290764213.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/rcomp/tmp/8fcwa1290764213.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/rcomp/tmp/9fcwa1290764213.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/rcomp/tmp/10qlvv1290764213.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11t4bj1290764213.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/rcomp/tmp/12fmap1290764213.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/rcomp/tmp/13exav1290764214.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/rcomp/tmp/14po9y1290764214.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/rcomp/tmp/15a7841290764214.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/rcomp/tmp/16e86a1290764214.tab")
+ }
>
> try(system("convert tmp/1j2yj1290764213.ps tmp/1j2yj1290764213.png",intern=TRUE))
character(0)
> try(system("convert tmp/2utf41290764213.ps tmp/2utf41290764213.png",intern=TRUE))
character(0)
> try(system("convert tmp/3utf41290764213.ps tmp/3utf41290764213.png",intern=TRUE))
character(0)
> try(system("convert tmp/4utf41290764213.ps tmp/4utf41290764213.png",intern=TRUE))
character(0)
> try(system("convert tmp/5utf41290764213.ps tmp/5utf41290764213.png",intern=TRUE))
character(0)
> try(system("convert tmp/6nlwp1290764213.ps tmp/6nlwp1290764213.png",intern=TRUE))
character(0)
> try(system("convert tmp/7fcwa1290764213.ps tmp/7fcwa1290764213.png",intern=TRUE))
character(0)
> try(system("convert tmp/8fcwa1290764213.ps tmp/8fcwa1290764213.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fcwa1290764213.ps tmp/9fcwa1290764213.png",intern=TRUE))
character(0)
> try(system("convert tmp/10qlvv1290764213.ps tmp/10qlvv1290764213.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.590 1.780 5.338