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 = '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
Energiedragers Invoer M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 -1.2 23.6 1 0 0 0 0 0 0 0 0 0 0
2 -2.4 25.7 0 1 0 0 0 0 0 0 0 0 0
3 0.8 32.5 0 0 1 0 0 0 0 0 0 0 0
4 -0.1 33.5 0 0 0 1 0 0 0 0 0 0 0
5 -1.5 34.5 0 0 0 0 1 0 0 0 0 0 0
6 -4.4 27.9 0 0 0 0 0 1 0 0 0 0 0
7 -4.2 45.3 0 0 0 0 0 0 1 0 0 0 0
8 3.5 40.8 0 0 0 0 0 0 0 1 0 0 0
9 10.0 58.5 0 0 0 0 0 0 0 0 1 0 0
10 8.6 32.5 0 0 0 0 0 0 0 0 0 1 0
11 9.5 35.5 0 0 0 0 0 0 0 0 0 0 1
12 9.9 46.7 0 0 0 0 0 0 0 0 0 0 0
13 10.4 53.2 1 0 0 0 0 0 0 0 0 0 0
14 16.0 36.1 0 1 0 0 0 0 0 0 0 0 0
15 12.7 54.0 0 0 1 0 0 0 0 0 0 0 0
16 10.2 58.1 0 0 0 1 0 0 0 0 0 0 0
17 8.9 41.8 0 0 0 0 1 0 0 0 0 0 0
18 12.6 43.1 0 0 0 0 0 1 0 0 0 0 0
19 13.6 76.0 0 0 0 0 0 0 1 0 0 0 0
20 14.8 42.8 0 0 0 0 0 0 0 1 0 0 0
21 9.5 41.0 0 0 0 0 0 0 0 0 1 0 0
22 13.7 61.4 0 0 0 0 0 0 0 0 0 1 0
23 17.0 34.2 0 0 0 0 0 0 0 0 0 0 1
24 14.7 53.8 0 0 0 0 0 0 0 0 0 0 0
25 17.4 80.7 1 0 0 0 0 0 0 0 0 0 0
26 9.0 79.5 0 1 0 0 0 0 0 0 0 0 0
27 9.1 96.5 0 0 1 0 0 0 0 0 0 0 0
28 12.2 108.3 0 0 0 1 0 0 0 0 0 0 0
29 15.9 100.1 0 0 0 0 1 0 0 0 0 0 0
30 12.9 108.5 0 0 0 0 0 1 0 0 0 0 0
31 10.9 127.4 0 0 0 0 0 0 1 0 0 0 0
32 10.6 86.5 0 0 0 0 0 0 0 1 0 0 0
33 13.2 71.4 0 0 0 0 0 0 0 0 1 0 0
34 9.6 88.2 0 0 0 0 0 0 0 0 0 1 0
35 6.4 135.6 0 0 0 0 0 0 0 0 0 0 1
36 5.8 70.5 0 0 0 0 0 0 0 0 0 0 0
37 -1.0 87.5 1 0 0 0 0 0 0 0 0 0 0
38 -0.2 73.3 0 1 0 0 0 0 0 0 0 0 0
39 2.7 92.2 0 0 1 0 0 0 0 0 0 0 0
40 3.6 61.1 0 0 0 1 0 0 0 0 0 0 0
41 -0.9 45.7 0 0 0 0 1 0 0 0 0 0 0
42 0.3 30.5 0 0 0 0 0 1 0 0 0 0 0
43 -1.1 34.8 0 0 0 0 0 0 1 0 0 0 0
44 -2.5 29.2 0 0 0 0 0 0 0 1 0 0 0
45 -3.4 56.7 0 0 0 0 0 0 0 0 1 0 0
46 -3.5 67.1 0 0 0 0 0 0 0 0 0 1 0
47 -3.9 41.8 0 0 0 0 0 0 0 0 0 0 1
48 -4.6 46.8 0 0 0 0 0 0 0 0 0 0 0
49 -0.1 50.1 1 0 0 0 0 0 0 0 0 0 0
50 4.3 81.9 0 1 0 0 0 0 0 0 0 0 0
51 10.2 115.8 0 0 1 0 0 0 0 0 0 0 0
52 8.7 102.5 0 0 0 1 0 0 0 0 0 0 0
53 13.3 106.6 0 0 0 0 1 0 0 0 0 0 0
54 15.0 101.4 0 0 0 0 0 1 0 0 0 0 0
55 20.7 136.1 0 0 0 0 0 0 1 0 0 0 0
56 20.7 143.4 0 0 0 0 0 0 0 1 0 0 0
57 26.4 127.5 0 0 0 0 0 0 0 0 1 0 0
58 31.2 113.8 0 0 0 0 0 0 0 0 0 1 0
59 31.4 75.3 0 0 0 0 0 0 0 0 0 0 1
60 26.6 98.5 0 0 0 0 0 0 0 0 0 0 0
61 26.6 113.7 1 0 0 0 0 0 0 0 0 0 0
62 19.2 103.7 0 1 0 0 0 0 0 0 0 0 0
63 6.5 73.9 0 0 1 0 0 0 0 0 0 0 0
64 3.1 52.5 0 0 0 1 0 0 0 0 0 0 0
65 -0.2 63.9 0 0 0 0 1 0 0 0 0 0 0
66 -4.0 44.9 0 0 0 0 0 1 0 0 0 0 0
67 -12.6 31.3 0 0 0 0 0 0 1 0 0 0 0
68 -13.0 24.9 0 0 0 0 0 0 0 1 0 0 0
69 -17.6 22.8 0 0 0 0 0 0 0 0 1 0 0
70 -21.7 24.8 0 0 0 0 0 0 0 0 0 1 0
71 -23.2 22.8 0 0 0 0 0 0 0 0 0 0 1
72 -16.8 20.9 0 0 0 0 0 0 0 0 0 0 0
73 -19.8 21.5 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Invoer M1 M2 M3 M4
-7.36747 0.23667 -2.56663 -0.76836 -3.97044 -2.75825
M5 M6 M7 M8 M9 M10
-2.20191 -1.28673 -5.86821 -1.44912 -1.18874 -1.61258
M11
-0.04889
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-21.1797 -5.8684 -0.3466 6.5752 20.9952
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.36747 4.34689 -1.695 0.0953 .
Invoer 0.23667 0.03430 6.900 3.7e-09 ***
M1 -2.56663 5.31263 -0.483 0.6308
M2 -0.76836 5.52174 -0.139 0.8898
M3 -3.97044 5.55813 -0.714 0.4778
M4 -2.75825 5.52836 -0.499 0.6197
M5 -2.20191 5.51907 -0.399 0.6913
M6 -1.28673 5.51106 -0.233 0.8162
M7 -5.86821 5.54818 -1.058 0.2944
M8 -1.44912 5.51272 -0.263 0.7936
M9 -1.18874 5.51489 -0.216 0.8301
M10 -1.61258 5.51757 -0.292 0.7711
M11 -0.04889 5.51017 -0.009 0.9929
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.544 on 60 degrees of freedom
Multiple R-squared: 0.4455, Adjusted R-squared: 0.3346
F-statistic: 4.017 on 12 and 60 DF, p-value: 0.0001463
> 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.18710504 0.37421007 0.8128950
[2,] 0.11983942 0.23967884 0.8801606
[3,] 0.09769804 0.19539609 0.9023020
[4,] 0.04816594 0.09633189 0.9518341
[5,] 0.05704398 0.11408796 0.9429560
[6,] 0.05570783 0.11141567 0.9442922
[7,] 0.05540590 0.11081180 0.9445941
[8,] 0.07132332 0.14264663 0.9286767
[9,] 0.04931977 0.09863955 0.9506802
[10,] 0.04179964 0.08359928 0.9582004
[11,] 0.08994446 0.17988892 0.9100555
[12,] 0.10262171 0.20524341 0.8973783
[13,] 0.07902934 0.15805867 0.9209707
[14,] 0.05024889 0.10049778 0.9497511
[15,] 0.03654117 0.07308235 0.9634588
[16,] 0.03008472 0.06016944 0.9699153
[17,] 0.02112890 0.04225781 0.9788711
[18,] 0.01466693 0.02933387 0.9853331
[19,] 0.01158714 0.02317428 0.9884129
[20,] 0.10150574 0.20301148 0.8984943
[21,] 0.08526584 0.17053168 0.9147342
[22,] 0.12711893 0.25423786 0.8728811
[23,] 0.11798805 0.23597610 0.8820120
[24,] 0.09640685 0.19281369 0.9035932
[25,] 0.06927579 0.13855158 0.9307242
[26,] 0.06044576 0.12089151 0.9395542
[27,] 0.05414669 0.10829339 0.9458533
[28,] 0.05736865 0.11473730 0.9426313
[29,] 0.08860746 0.17721492 0.9113925
[30,] 0.09333165 0.18666330 0.9066683
[31,] 0.10722840 0.21445680 0.8927716
[32,] 0.10022381 0.20044762 0.8997762
[33,] 0.09394080 0.18788160 0.9060592
[34,] 0.06601954 0.13203908 0.9339805
[35,] 0.04453194 0.08906389 0.9554681
[36,] 0.03964999 0.07929999 0.9603500
[37,] 0.04127078 0.08254156 0.9587292
[38,] 0.02479512 0.04959023 0.9752049
[39,] 0.01466036 0.02932072 0.9853396
[40,] 0.01581868 0.03163735 0.9841813
[41,] 0.07890448 0.15780896 0.9210955
[42,] 0.11572750 0.23145499 0.8842725
> postscript(file="/var/www/html/rcomp/tmp/16eok1261311652.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/2zkj41261311652.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/3ydk21261311652.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/4pise1261311652.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/5xpq31261311652.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
3.14870871 -0.34656960 4.44616205 2.09730706 -0.09570454 -2.34886991
7 8 9 10 11 12
-1.68542933 2.66049290 4.71106617 9.88829830 8.51460772 6.21502249
13 14 15 16 17 18
7.74330530 15.59207244 11.25777782 6.57524882 8.57661151 11.05376077
19 20 21 22 23 24
8.84883132 13.48715483 8.35277427 8.14856321 16.32227746 9.33467235
25 26 27 28 29 30
8.23490686 -1.67936364 -2.40065613 -3.30553670 1.77880682 -4.12439406
31 32 33 34 35 36
-6.01595703 -1.05528196 4.85803563 -2.29416690 -18.27596260 -3.51770052
37 38 39 40 41 42
-11.77444257 -9.41201563 -7.78297929 -0.73475829 -2.14639773 1.73579060
43 44 45 46 47 48
3.89959553 -0.59414631 -8.26292957 -10.40045028 -6.37640720 -8.30864441
49 50 51 52 53 54
-2.02302070 -6.94736932 -5.86836849 -5.43285630 -2.35954191 -0.34404392
55 56 57 58 59 60
1.72502237 -4.42175000 4.78090281 13.24710582 20.99518016 10.65556652
61 62 63 64 65 66
9.62482873 2.79324574 0.34806404 0.80059541 -5.75377415 -5.97224349
67 68 69 70 71 72
-6.77206285 -10.07646946 -14.43984931 -18.58935014 -21.17969555 -14.37891643
73
-14.95428632
> postscript(file="/var/www/html/rcomp/tmp/6efkw1261311652.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 3.14870871 NA
1 -0.34656960 3.14870871
2 4.44616205 -0.34656960
3 2.09730706 4.44616205
4 -0.09570454 2.09730706
5 -2.34886991 -0.09570454
6 -1.68542933 -2.34886991
7 2.66049290 -1.68542933
8 4.71106617 2.66049290
9 9.88829830 4.71106617
10 8.51460772 9.88829830
11 6.21502249 8.51460772
12 7.74330530 6.21502249
13 15.59207244 7.74330530
14 11.25777782 15.59207244
15 6.57524882 11.25777782
16 8.57661151 6.57524882
17 11.05376077 8.57661151
18 8.84883132 11.05376077
19 13.48715483 8.84883132
20 8.35277427 13.48715483
21 8.14856321 8.35277427
22 16.32227746 8.14856321
23 9.33467235 16.32227746
24 8.23490686 9.33467235
25 -1.67936364 8.23490686
26 -2.40065613 -1.67936364
27 -3.30553670 -2.40065613
28 1.77880682 -3.30553670
29 -4.12439406 1.77880682
30 -6.01595703 -4.12439406
31 -1.05528196 -6.01595703
32 4.85803563 -1.05528196
33 -2.29416690 4.85803563
34 -18.27596260 -2.29416690
35 -3.51770052 -18.27596260
36 -11.77444257 -3.51770052
37 -9.41201563 -11.77444257
38 -7.78297929 -9.41201563
39 -0.73475829 -7.78297929
40 -2.14639773 -0.73475829
41 1.73579060 -2.14639773
42 3.89959553 1.73579060
43 -0.59414631 3.89959553
44 -8.26292957 -0.59414631
45 -10.40045028 -8.26292957
46 -6.37640720 -10.40045028
47 -8.30864441 -6.37640720
48 -2.02302070 -8.30864441
49 -6.94736932 -2.02302070
50 -5.86836849 -6.94736932
51 -5.43285630 -5.86836849
52 -2.35954191 -5.43285630
53 -0.34404392 -2.35954191
54 1.72502237 -0.34404392
55 -4.42175000 1.72502237
56 4.78090281 -4.42175000
57 13.24710582 4.78090281
58 20.99518016 13.24710582
59 10.65556652 20.99518016
60 9.62482873 10.65556652
61 2.79324574 9.62482873
62 0.34806404 2.79324574
63 0.80059541 0.34806404
64 -5.75377415 0.80059541
65 -5.97224349 -5.75377415
66 -6.77206285 -5.97224349
67 -10.07646946 -6.77206285
68 -14.43984931 -10.07646946
69 -18.58935014 -14.43984931
70 -21.17969555 -18.58935014
71 -14.37891643 -21.17969555
72 -14.95428632 -14.37891643
73 NA -14.95428632
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.34656960 3.14870871
[2,] 4.44616205 -0.34656960
[3,] 2.09730706 4.44616205
[4,] -0.09570454 2.09730706
[5,] -2.34886991 -0.09570454
[6,] -1.68542933 -2.34886991
[7,] 2.66049290 -1.68542933
[8,] 4.71106617 2.66049290
[9,] 9.88829830 4.71106617
[10,] 8.51460772 9.88829830
[11,] 6.21502249 8.51460772
[12,] 7.74330530 6.21502249
[13,] 15.59207244 7.74330530
[14,] 11.25777782 15.59207244
[15,] 6.57524882 11.25777782
[16,] 8.57661151 6.57524882
[17,] 11.05376077 8.57661151
[18,] 8.84883132 11.05376077
[19,] 13.48715483 8.84883132
[20,] 8.35277427 13.48715483
[21,] 8.14856321 8.35277427
[22,] 16.32227746 8.14856321
[23,] 9.33467235 16.32227746
[24,] 8.23490686 9.33467235
[25,] -1.67936364 8.23490686
[26,] -2.40065613 -1.67936364
[27,] -3.30553670 -2.40065613
[28,] 1.77880682 -3.30553670
[29,] -4.12439406 1.77880682
[30,] -6.01595703 -4.12439406
[31,] -1.05528196 -6.01595703
[32,] 4.85803563 -1.05528196
[33,] -2.29416690 4.85803563
[34,] -18.27596260 -2.29416690
[35,] -3.51770052 -18.27596260
[36,] -11.77444257 -3.51770052
[37,] -9.41201563 -11.77444257
[38,] -7.78297929 -9.41201563
[39,] -0.73475829 -7.78297929
[40,] -2.14639773 -0.73475829
[41,] 1.73579060 -2.14639773
[42,] 3.89959553 1.73579060
[43,] -0.59414631 3.89959553
[44,] -8.26292957 -0.59414631
[45,] -10.40045028 -8.26292957
[46,] -6.37640720 -10.40045028
[47,] -8.30864441 -6.37640720
[48,] -2.02302070 -8.30864441
[49,] -6.94736932 -2.02302070
[50,] -5.86836849 -6.94736932
[51,] -5.43285630 -5.86836849
[52,] -2.35954191 -5.43285630
[53,] -0.34404392 -2.35954191
[54,] 1.72502237 -0.34404392
[55,] -4.42175000 1.72502237
[56,] 4.78090281 -4.42175000
[57,] 13.24710582 4.78090281
[58,] 20.99518016 13.24710582
[59,] 10.65556652 20.99518016
[60,] 9.62482873 10.65556652
[61,] 2.79324574 9.62482873
[62,] 0.34806404 2.79324574
[63,] 0.80059541 0.34806404
[64,] -5.75377415 0.80059541
[65,] -5.97224349 -5.75377415
[66,] -6.77206285 -5.97224349
[67,] -10.07646946 -6.77206285
[68,] -14.43984931 -10.07646946
[69,] -18.58935014 -14.43984931
[70,] -21.17969555 -18.58935014
[71,] -14.37891643 -21.17969555
[72,] -14.95428632 -14.37891643
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.34656960 3.14870871
2 4.44616205 -0.34656960
3 2.09730706 4.44616205
4 -0.09570454 2.09730706
5 -2.34886991 -0.09570454
6 -1.68542933 -2.34886991
7 2.66049290 -1.68542933
8 4.71106617 2.66049290
9 9.88829830 4.71106617
10 8.51460772 9.88829830
11 6.21502249 8.51460772
12 7.74330530 6.21502249
13 15.59207244 7.74330530
14 11.25777782 15.59207244
15 6.57524882 11.25777782
16 8.57661151 6.57524882
17 11.05376077 8.57661151
18 8.84883132 11.05376077
19 13.48715483 8.84883132
20 8.35277427 13.48715483
21 8.14856321 8.35277427
22 16.32227746 8.14856321
23 9.33467235 16.32227746
24 8.23490686 9.33467235
25 -1.67936364 8.23490686
26 -2.40065613 -1.67936364
27 -3.30553670 -2.40065613
28 1.77880682 -3.30553670
29 -4.12439406 1.77880682
30 -6.01595703 -4.12439406
31 -1.05528196 -6.01595703
32 4.85803563 -1.05528196
33 -2.29416690 4.85803563
34 -18.27596260 -2.29416690
35 -3.51770052 -18.27596260
36 -11.77444257 -3.51770052
37 -9.41201563 -11.77444257
38 -7.78297929 -9.41201563
39 -0.73475829 -7.78297929
40 -2.14639773 -0.73475829
41 1.73579060 -2.14639773
42 3.89959553 1.73579060
43 -0.59414631 3.89959553
44 -8.26292957 -0.59414631
45 -10.40045028 -8.26292957
46 -6.37640720 -10.40045028
47 -8.30864441 -6.37640720
48 -2.02302070 -8.30864441
49 -6.94736932 -2.02302070
50 -5.86836849 -6.94736932
51 -5.43285630 -5.86836849
52 -2.35954191 -5.43285630
53 -0.34404392 -2.35954191
54 1.72502237 -0.34404392
55 -4.42175000 1.72502237
56 4.78090281 -4.42175000
57 13.24710582 4.78090281
58 20.99518016 13.24710582
59 10.65556652 20.99518016
60 9.62482873 10.65556652
61 2.79324574 9.62482873
62 0.34806404 2.79324574
63 0.80059541 0.34806404
64 -5.75377415 0.80059541
65 -5.97224349 -5.75377415
66 -6.77206285 -5.97224349
67 -10.07646946 -6.77206285
68 -14.43984931 -10.07646946
69 -18.58935014 -14.43984931
70 -21.17969555 -18.58935014
71 -14.37891643 -21.17969555
72 -14.95428632 -14.37891643
> 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/7zcc71261311652.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/8ogco1261311652.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/9fpot1261311652.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/10kvm51261311652.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/11yv1o1261311652.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/12woz61261311652.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/13qjqq1261311652.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/14gtvp1261311652.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/159kws1261311652.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/168tuy1261311652.tab")
+ }
>
> try(system("convert tmp/16eok1261311652.ps tmp/16eok1261311652.png",intern=TRUE))
character(0)
> try(system("convert tmp/2zkj41261311652.ps tmp/2zkj41261311652.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ydk21261311652.ps tmp/3ydk21261311652.png",intern=TRUE))
character(0)
> try(system("convert tmp/4pise1261311652.ps tmp/4pise1261311652.png",intern=TRUE))
character(0)
> try(system("convert tmp/5xpq31261311652.ps tmp/5xpq31261311652.png",intern=TRUE))
character(0)
> try(system("convert tmp/6efkw1261311652.ps tmp/6efkw1261311652.png",intern=TRUE))
character(0)
> try(system("convert tmp/7zcc71261311652.ps tmp/7zcc71261311652.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ogco1261311652.ps tmp/8ogco1261311652.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fpot1261311652.ps tmp/9fpot1261311652.png",intern=TRUE))
character(0)
> try(system("convert tmp/10kvm51261311652.ps tmp/10kvm51261311652.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.532 1.573 3.317