R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(-1.2,23.6,-2.4,25.7,0.8,32.5,-0.1,33.5,-1.5,34.5,-4.4,27.9,-4.2,45.3,3.5,40.8,10,58.5,8.6,32.5,9.5,35.5,9.9,46.7,10.4,53.2,16,36.1,12.7,54,10.2,58.1,8.9,41.8,12.6,43.1,13.6,76,14.8,42.8,9.5,41,13.7,61.4,17,34.2,14.7,53.8,17.4,80.7,9,79.5,9.1,96.5,12.2,108.3,15.9,100.1,12.9,108.5,10.9,127.4,10.6,86.5,13.2,71.4,9.6,88.2,6.4,135.6,5.8,70.5,-1,87.5,-0.2,73.3,2.7,92.2,3.6,61.1,-0.9,45.7,0.3,30.5,-1.1,34.8,-2.5,29.2,-3.4,56.7,-3.5,67.1,-3.9,41.8,-4.6,46.8,-0.1,50.1,4.3,81.9,10.2,115.8,8.7,102.5,13.3,106.6,15,101.4,20.7,136.1,20.7,143.4,26.4,127.5,31.2,113.8,31.4,75.3,26.6,98.5,26.6,113.7,19.2,103.7,6.5,73.9,3.1,52.5,-0.2,63.9,-4,44.9,-12.6,31.3,-13,24.9,-17.6,22.8,-21.7,24.8,-23.2,22.8,-16.8,20.9,-19.8,21.5),dim=c(2,73),dimnames=list(c('Energiedragers','Invoer'),1:73))
> y <- array(NA,dim=c(2,73),dimnames=list(c('Energiedragers','Invoer'),1:73))
> 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
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
Energiedragers Invoer
1 -1.2 23.6
2 -2.4 25.7
3 0.8 32.5
4 -0.1 33.5
5 -1.5 34.5
6 -4.4 27.9
7 -4.2 45.3
8 3.5 40.8
9 10.0 58.5
10 8.6 32.5
11 9.5 35.5
12 9.9 46.7
13 10.4 53.2
14 16.0 36.1
15 12.7 54.0
16 10.2 58.1
17 8.9 41.8
18 12.6 43.1
19 13.6 76.0
20 14.8 42.8
21 9.5 41.0
22 13.7 61.4
23 17.0 34.2
24 14.7 53.8
25 17.4 80.7
26 9.0 79.5
27 9.1 96.5
28 12.2 108.3
29 15.9 100.1
30 12.9 108.5
31 10.9 127.4
32 10.6 86.5
33 13.2 71.4
34 9.6 88.2
35 6.4 135.6
36 5.8 70.5
37 -1.0 87.5
38 -0.2 73.3
39 2.7 92.2
40 3.6 61.1
41 -0.9 45.7
42 0.3 30.5
43 -1.1 34.8
44 -2.5 29.2
45 -3.4 56.7
46 -3.5 67.1
47 -3.9 41.8
48 -4.6 46.8
49 -0.1 50.1
50 4.3 81.9
51 10.2 115.8
52 8.7 102.5
53 13.3 106.6
54 15.0 101.4
55 20.7 136.1
56 20.7 143.4
57 26.4 127.5
58 31.2 113.8
59 31.4 75.3
60 26.6 98.5
61 26.6 113.7
62 19.2 103.7
63 6.5 73.9
64 3.1 52.5
65 -0.2 63.9
66 -4.0 44.9
67 -12.6 31.3
68 -13.0 24.9
69 -17.6 22.8
70 -21.7 24.8
71 -23.2 22.8
72 -16.8 20.9
73 -19.8 21.5
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Invoer
-8.855 0.229
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-19.5660 -5.7183 -0.3526 6.0588 23.0120
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.85490 2.28930 -3.868 0.000241 ***
Invoer 0.22899 0.03147 7.277 3.64e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.916 on 71 degrees of freedom
Multiple R-squared: 0.4272, Adjusted R-squared: 0.4191
F-statistic: 52.95 on 1 and 71 DF, p-value: 3.637e-10
> 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.003287994 0.006575987 0.996712006
[2,] 0.002791959 0.005583918 0.997208041
[3,] 0.001610749 0.003221499 0.998389251
[4,] 0.002673817 0.005347633 0.997326183
[5,] 0.003862024 0.007724048 0.996137976
[6,] 0.014602927 0.029205854 0.985397073
[7,] 0.023127226 0.046254453 0.976872774
[8,] 0.016769408 0.033538816 0.983230592
[9,] 0.009269545 0.018539090 0.990730455
[10,] 0.049106218 0.098212436 0.950893782
[11,] 0.034349637 0.068699274 0.965650363
[12,] 0.020865378 0.041730757 0.979134622
[13,] 0.014606188 0.029212375 0.985393812
[14,] 0.015521770 0.031043539 0.984478230
[15,] 0.010626760 0.021253519 0.989373240
[16,] 0.018195761 0.036391521 0.981804239
[17,] 0.015084029 0.030168057 0.984915971
[18,] 0.010979710 0.021959419 0.989020290
[19,] 0.059130078 0.118260156 0.940869922
[20,] 0.067382723 0.134765445 0.932617277
[21,] 0.056822835 0.113645670 0.943177165
[22,] 0.059820383 0.119640766 0.940179617
[23,] 0.069010444 0.138020888 0.930989556
[24,] 0.057948291 0.115896582 0.942051709
[25,] 0.040516149 0.081032297 0.959483851
[26,] 0.029688753 0.059377506 0.970311247
[27,] 0.033875065 0.067750129 0.966124935
[28,] 0.022945240 0.045890481 0.977054760
[29,] 0.019679028 0.039358056 0.980320972
[30,] 0.013277808 0.026555617 0.986722192
[31,] 0.042844998 0.085689996 0.957155002
[32,] 0.032071771 0.064143542 0.967928229
[33,] 0.062203982 0.124407965 0.937796018
[34,] 0.070595878 0.141191757 0.929404122
[35,] 0.081270719 0.162541438 0.918729281
[36,] 0.064114724 0.128229448 0.935885276
[37,] 0.060078331 0.120156662 0.939921669
[38,] 0.066929191 0.133858382 0.933070809
[39,] 0.069729602 0.139459204 0.930270398
[40,] 0.080841995 0.161683991 0.919158005
[41,] 0.082811892 0.165623785 0.917188108
[42,] 0.095912423 0.191824845 0.904087577
[43,] 0.090694118 0.181388235 0.909305882
[44,] 0.085310665 0.170621330 0.914689335
[45,] 0.070003446 0.140006891 0.929996554
[46,] 0.053420112 0.106840224 0.946579888
[47,] 0.061047469 0.122094938 0.938952531
[48,] 0.056385085 0.112770171 0.943614915
[49,] 0.046206377 0.092412753 0.953793623
[50,] 0.034075512 0.068151024 0.965924488
[51,] 0.054519944 0.109039888 0.945480056
[52,] 0.251726543 0.503453085 0.748273457
[53,] 0.382971446 0.765942891 0.617028554
[54,] 0.412904118 0.825808237 0.587095882
[55,] 0.990600377 0.018799246 0.009399623
[56,] 0.995700623 0.008598754 0.004299377
[57,] 0.991691162 0.016617677 0.008308838
[58,] 0.986104198 0.027791605 0.013895802
[59,] 0.974993854 0.050012291 0.025006146
[60,] 0.978951939 0.042096122 0.021048061
[61,] 0.974316891 0.051366218 0.025683109
[62,] 0.941603171 0.116793657 0.058396829
[63,] 0.888579470 0.222841061 0.111420530
[64,] 0.957738068 0.084523863 0.042261932
> postscript(file="/var/www/html/rcomp/tmp/1rqwm1261301633.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/2o2p21261301633.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/3hhp01261301633.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/4ptpu1261301633.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/5do9k1261301633.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 = 73
Frequency = 1
1 2 3 4 5 6
2.25076246 0.56988569 2.21276094 1.08377200 -0.54521693 -1.93388996
7 8 9 10 11 12
-5.71829743 3.01215278 5.45904863 10.01276094 10.22579413 8.06111806
13 14 15 16 17 18
7.07268998 16.58840077 9.18949884 5.75064420 8.18316384 11.58547823
19 20 21 22 23 24
5.05174227 13.85417491 8.96635499 8.49498072 18.02347975 11.23529662
25 26 27 28 29 30
7.77549427 -0.34971900 -4.14253090 -3.74460033 1.83310894 -3.09039812
31 32 33 34 35 36
-9.41828899 -0.35264155 5.70509137 -1.74192274 -15.79599825 -1.48881859
37 38 39 40 41 42
-12.18163048 -8.12998761 -9.55787848 -1.53632260 -2.50989300 2.17073881
43 44 45 46 47 48
-0.21391361 -0.33157558 -7.52877129 -10.01025621 -4.61683616 -6.46178083
49 50 51 52 53 54
-2.71744432 -5.59929245 -7.46201734 -5.91646451 -2.25531914 0.63542332
55 56 57 58 59 60
-1.61049272 -3.28211195 6.05881212 13.99596053 23.01203452 12.89949123
61 62 63 64 65 66
9.41885942 4.30874877 -1.56738097 -0.06701776 -5.97749162 -5.42670186
67 68 69 70 71 72
-10.91245234 -9.84692316 -13.96604639 -18.52402426 -19.56604639 -12.73096742
73
-15.86836078
> postscript(file="/var/www/html/rcomp/tmp/6l3hs1261301633.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 2.25076246 NA
1 0.56988569 2.25076246
2 2.21276094 0.56988569
3 1.08377200 2.21276094
4 -0.54521693 1.08377200
5 -1.93388996 -0.54521693
6 -5.71829743 -1.93388996
7 3.01215278 -5.71829743
8 5.45904863 3.01215278
9 10.01276094 5.45904863
10 10.22579413 10.01276094
11 8.06111806 10.22579413
12 7.07268998 8.06111806
13 16.58840077 7.07268998
14 9.18949884 16.58840077
15 5.75064420 9.18949884
16 8.18316384 5.75064420
17 11.58547823 8.18316384
18 5.05174227 11.58547823
19 13.85417491 5.05174227
20 8.96635499 13.85417491
21 8.49498072 8.96635499
22 18.02347975 8.49498072
23 11.23529662 18.02347975
24 7.77549427 11.23529662
25 -0.34971900 7.77549427
26 -4.14253090 -0.34971900
27 -3.74460033 -4.14253090
28 1.83310894 -3.74460033
29 -3.09039812 1.83310894
30 -9.41828899 -3.09039812
31 -0.35264155 -9.41828899
32 5.70509137 -0.35264155
33 -1.74192274 5.70509137
34 -15.79599825 -1.74192274
35 -1.48881859 -15.79599825
36 -12.18163048 -1.48881859
37 -8.12998761 -12.18163048
38 -9.55787848 -8.12998761
39 -1.53632260 -9.55787848
40 -2.50989300 -1.53632260
41 2.17073881 -2.50989300
42 -0.21391361 2.17073881
43 -0.33157558 -0.21391361
44 -7.52877129 -0.33157558
45 -10.01025621 -7.52877129
46 -4.61683616 -10.01025621
47 -6.46178083 -4.61683616
48 -2.71744432 -6.46178083
49 -5.59929245 -2.71744432
50 -7.46201734 -5.59929245
51 -5.91646451 -7.46201734
52 -2.25531914 -5.91646451
53 0.63542332 -2.25531914
54 -1.61049272 0.63542332
55 -3.28211195 -1.61049272
56 6.05881212 -3.28211195
57 13.99596053 6.05881212
58 23.01203452 13.99596053
59 12.89949123 23.01203452
60 9.41885942 12.89949123
61 4.30874877 9.41885942
62 -1.56738097 4.30874877
63 -0.06701776 -1.56738097
64 -5.97749162 -0.06701776
65 -5.42670186 -5.97749162
66 -10.91245234 -5.42670186
67 -9.84692316 -10.91245234
68 -13.96604639 -9.84692316
69 -18.52402426 -13.96604639
70 -19.56604639 -18.52402426
71 -12.73096742 -19.56604639
72 -15.86836078 -12.73096742
73 NA -15.86836078
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.56988569 2.25076246
[2,] 2.21276094 0.56988569
[3,] 1.08377200 2.21276094
[4,] -0.54521693 1.08377200
[5,] -1.93388996 -0.54521693
[6,] -5.71829743 -1.93388996
[7,] 3.01215278 -5.71829743
[8,] 5.45904863 3.01215278
[9,] 10.01276094 5.45904863
[10,] 10.22579413 10.01276094
[11,] 8.06111806 10.22579413
[12,] 7.07268998 8.06111806
[13,] 16.58840077 7.07268998
[14,] 9.18949884 16.58840077
[15,] 5.75064420 9.18949884
[16,] 8.18316384 5.75064420
[17,] 11.58547823 8.18316384
[18,] 5.05174227 11.58547823
[19,] 13.85417491 5.05174227
[20,] 8.96635499 13.85417491
[21,] 8.49498072 8.96635499
[22,] 18.02347975 8.49498072
[23,] 11.23529662 18.02347975
[24,] 7.77549427 11.23529662
[25,] -0.34971900 7.77549427
[26,] -4.14253090 -0.34971900
[27,] -3.74460033 -4.14253090
[28,] 1.83310894 -3.74460033
[29,] -3.09039812 1.83310894
[30,] -9.41828899 -3.09039812
[31,] -0.35264155 -9.41828899
[32,] 5.70509137 -0.35264155
[33,] -1.74192274 5.70509137
[34,] -15.79599825 -1.74192274
[35,] -1.48881859 -15.79599825
[36,] -12.18163048 -1.48881859
[37,] -8.12998761 -12.18163048
[38,] -9.55787848 -8.12998761
[39,] -1.53632260 -9.55787848
[40,] -2.50989300 -1.53632260
[41,] 2.17073881 -2.50989300
[42,] -0.21391361 2.17073881
[43,] -0.33157558 -0.21391361
[44,] -7.52877129 -0.33157558
[45,] -10.01025621 -7.52877129
[46,] -4.61683616 -10.01025621
[47,] -6.46178083 -4.61683616
[48,] -2.71744432 -6.46178083
[49,] -5.59929245 -2.71744432
[50,] -7.46201734 -5.59929245
[51,] -5.91646451 -7.46201734
[52,] -2.25531914 -5.91646451
[53,] 0.63542332 -2.25531914
[54,] -1.61049272 0.63542332
[55,] -3.28211195 -1.61049272
[56,] 6.05881212 -3.28211195
[57,] 13.99596053 6.05881212
[58,] 23.01203452 13.99596053
[59,] 12.89949123 23.01203452
[60,] 9.41885942 12.89949123
[61,] 4.30874877 9.41885942
[62,] -1.56738097 4.30874877
[63,] -0.06701776 -1.56738097
[64,] -5.97749162 -0.06701776
[65,] -5.42670186 -5.97749162
[66,] -10.91245234 -5.42670186
[67,] -9.84692316 -10.91245234
[68,] -13.96604639 -9.84692316
[69,] -18.52402426 -13.96604639
[70,] -19.56604639 -18.52402426
[71,] -12.73096742 -19.56604639
[72,] -15.86836078 -12.73096742
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.56988569 2.25076246
2 2.21276094 0.56988569
3 1.08377200 2.21276094
4 -0.54521693 1.08377200
5 -1.93388996 -0.54521693
6 -5.71829743 -1.93388996
7 3.01215278 -5.71829743
8 5.45904863 3.01215278
9 10.01276094 5.45904863
10 10.22579413 10.01276094
11 8.06111806 10.22579413
12 7.07268998 8.06111806
13 16.58840077 7.07268998
14 9.18949884 16.58840077
15 5.75064420 9.18949884
16 8.18316384 5.75064420
17 11.58547823 8.18316384
18 5.05174227 11.58547823
19 13.85417491 5.05174227
20 8.96635499 13.85417491
21 8.49498072 8.96635499
22 18.02347975 8.49498072
23 11.23529662 18.02347975
24 7.77549427 11.23529662
25 -0.34971900 7.77549427
26 -4.14253090 -0.34971900
27 -3.74460033 -4.14253090
28 1.83310894 -3.74460033
29 -3.09039812 1.83310894
30 -9.41828899 -3.09039812
31 -0.35264155 -9.41828899
32 5.70509137 -0.35264155
33 -1.74192274 5.70509137
34 -15.79599825 -1.74192274
35 -1.48881859 -15.79599825
36 -12.18163048 -1.48881859
37 -8.12998761 -12.18163048
38 -9.55787848 -8.12998761
39 -1.53632260 -9.55787848
40 -2.50989300 -1.53632260
41 2.17073881 -2.50989300
42 -0.21391361 2.17073881
43 -0.33157558 -0.21391361
44 -7.52877129 -0.33157558
45 -10.01025621 -7.52877129
46 -4.61683616 -10.01025621
47 -6.46178083 -4.61683616
48 -2.71744432 -6.46178083
49 -5.59929245 -2.71744432
50 -7.46201734 -5.59929245
51 -5.91646451 -7.46201734
52 -2.25531914 -5.91646451
53 0.63542332 -2.25531914
54 -1.61049272 0.63542332
55 -3.28211195 -1.61049272
56 6.05881212 -3.28211195
57 13.99596053 6.05881212
58 23.01203452 13.99596053
59 12.89949123 23.01203452
60 9.41885942 12.89949123
61 4.30874877 9.41885942
62 -1.56738097 4.30874877
63 -0.06701776 -1.56738097
64 -5.97749162 -0.06701776
65 -5.42670186 -5.97749162
66 -10.91245234 -5.42670186
67 -9.84692316 -10.91245234
68 -13.96604639 -9.84692316
69 -18.52402426 -13.96604639
70 -19.56604639 -18.52402426
71 -12.73096742 -19.56604639
72 -15.86836078 -12.73096742
> 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/79d2t1261301633.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/8rszz1261301633.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/9rc8c1261301633.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/102q1a1261301633.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/11nvre1261301633.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/12r3ug1261301633.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/1365oh1261301633.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/14d3yv1261301633.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/156h0z1261301633.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/16vn6j1261301633.tab")
+ }
>
> try(system("convert tmp/1rqwm1261301633.ps tmp/1rqwm1261301633.png",intern=TRUE))
character(0)
> try(system("convert tmp/2o2p21261301633.ps tmp/2o2p21261301633.png",intern=TRUE))
character(0)
> try(system("convert tmp/3hhp01261301633.ps tmp/3hhp01261301633.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ptpu1261301633.ps tmp/4ptpu1261301633.png",intern=TRUE))
character(0)
> try(system("convert tmp/5do9k1261301633.ps tmp/5do9k1261301633.png",intern=TRUE))
character(0)
> try(system("convert tmp/6l3hs1261301633.ps tmp/6l3hs1261301633.png",intern=TRUE))
character(0)
> try(system("convert tmp/79d2t1261301633.ps tmp/79d2t1261301633.png",intern=TRUE))
character(0)
> try(system("convert tmp/8rszz1261301633.ps tmp/8rszz1261301633.png",intern=TRUE))
character(0)
> try(system("convert tmp/9rc8c1261301633.ps tmp/9rc8c1261301633.png",intern=TRUE))
character(0)
> try(system("convert tmp/102q1a1261301633.ps tmp/102q1a1261301633.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.401 1.467 3.480