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(102.89
+ ,167.16
+ ,100.70
+ ,106.88
+ ,97.69
+ ,102.64
+ ,179.84
+ ,99.62
+ ,107.45
+ ,101.69
+ ,103.33
+ ,174.44
+ ,99.83
+ ,107.65
+ ,102.72
+ ,103.56
+ ,180.35
+ ,100.74
+ ,107.72
+ ,101.85
+ ,103.60
+ ,193.17
+ ,100.84
+ ,108.10
+ ,114.94
+ ,104.24
+ ,195.16
+ ,100.85
+ ,108.38
+ ,106.20
+ ,105.31
+ ,202.43
+ ,99.71
+ ,108.62
+ ,106.76
+ ,105.40
+ ,189.91
+ ,100.80
+ ,108.79
+ ,107.24
+ ,105.89
+ ,195.98
+ ,100.06
+ ,109.03
+ ,106.50
+ ,105.89
+ ,212.09
+ ,100.57
+ ,109.34
+ ,106.77
+ ,105.54
+ ,205.81
+ ,99.79
+ ,109.73
+ ,108.24
+ ,106.15
+ ,204.31
+ ,99.90
+ ,109.76
+ ,104.43
+ ,106.14
+ ,196.07
+ ,100.12
+ ,109.96
+ ,100.90
+ ,105.85
+ ,199.98
+ ,100.40
+ ,110.49
+ ,103.91
+ ,106.27
+ ,199.10
+ ,100.51
+ ,111.37
+ ,103.81
+ ,106.51
+ ,198.31
+ ,100.70
+ ,111.56
+ ,104.59
+ ,106.82
+ ,195.72
+ ,100.62
+ ,111.90
+ ,104.94
+ ,106.53
+ ,223.04
+ ,99.70
+ ,111.96
+ ,111.64
+ ,107.14
+ ,238.41
+ ,99.48
+ ,112.25
+ ,111.27
+ ,107.39
+ ,259.73
+ ,99.36
+ ,112.39
+ ,106.82
+ ,107.33
+ ,326.54
+ ,99.39
+ ,112.30
+ ,106.07
+ ,107.53
+ ,335.15
+ ,99.45
+ ,112.49
+ ,111.35
+ ,107.42
+ ,321.81
+ ,99.28
+ ,112.77
+ ,112.59
+ ,108.25
+ ,368.62
+ ,99.40
+ ,113.15
+ ,108.59
+ ,108.26
+ ,369.59
+ ,99.10
+ ,113.15
+ ,106.83
+ ,108.93
+ ,425.00
+ ,99.48
+ ,113.28
+ ,112.51
+ ,109.43
+ ,439.72
+ ,99.74
+ ,113.83
+ ,113.61
+ ,109.61
+ ,362.23
+ ,100.42
+ ,114.49
+ ,114.96
+ ,109.74
+ ,328.76
+ ,100.80
+ ,114.76
+ ,118.66
+ ,110.12
+ ,348.55
+ ,100.66
+ ,114.96
+ ,116.84
+ ,110.16
+ ,328.18
+ ,101.03
+ ,115.41
+ ,121.19
+ ,110.44
+ ,329.34
+ ,101.22
+ ,115.84
+ ,117.42
+ ,111.23
+ ,295.55
+ ,101.23
+ ,116.31
+ ,116.88
+ ,112.86
+ ,237.38
+ ,100.10
+ ,117.23
+ ,115.01
+ ,112.77
+ ,226.85
+ ,99.98
+ ,117.97
+ ,111.81
+ ,113.04
+ ,220.14
+ ,99.91
+ ,118.08
+ ,110.61
+ ,112.79
+ ,239.36
+ ,99.84
+ ,118.27
+ ,110.67
+ ,113.87
+ ,224.69
+ ,99.68
+ ,118.88
+ ,113.28
+ ,114.28
+ ,230.98
+ ,99.74
+ ,119.11
+ ,112.08
+ ,115.51
+ ,233.47
+ ,99.71
+ ,119.29
+ ,111.41
+ ,116.76
+ ,256.70
+ ,99.35
+ ,119.36
+ ,113.81
+ ,116.91
+ ,253.41
+ ,99.21
+ ,119.48
+ ,109.16
+ ,116.47
+ ,224.95
+ ,99.21
+ ,120.10
+ ,105.09
+ ,116.94
+ ,210.37
+ ,99.16
+ ,120.30
+ ,102.23
+ ,117.24
+ ,191.09
+ ,99.20
+ ,120.54
+ ,101.95
+ ,116.82
+ ,198.85
+ ,99.08
+ ,120.86
+ ,104.75
+ ,117.48
+ ,211.04
+ ,98.16
+ ,121.10
+ ,107.25
+ ,117.11
+ ,206.25
+ ,98.00
+ ,121.42
+ ,105.25
+ ,117.31
+ ,201.19
+ ,97.90
+ ,121.81
+ ,102.75
+ ,117.77
+ ,194.37
+ ,97.88
+ ,122.21
+ ,107.21
+ ,118.37
+ ,191.08
+ ,97.56
+ ,122.82
+ ,107.24
+ ,117.91
+ ,192.87
+ ,96.86
+ ,123.02
+ ,106.01
+ ,118.12
+ ,181.61
+ ,96.86
+ ,123.14
+ ,121.36
+ ,118.02
+ ,157.67
+ ,96.75
+ ,123.12
+ ,120.44
+ ,117.77
+ ,196.14
+ ,97.12
+ ,123.42
+ ,109.40
+ ,117.85
+ ,246.35
+ ,97.22
+ ,123.50
+ ,111.51
+ ,118.68
+ ,271.90
+ ,97.52
+ ,125.77
+ ,111.97
+ ,118.90
+ ,270.29
+ ,97.57
+ ,125.99
+ ,114.64)
+ ,dim=c(5
+ ,58)
+ ,dimnames=list(c('Bier'
+ ,'Tarwe'
+ ,'Suiker'
+ ,'Mineraalwater'
+ ,'Fruit')
+ ,1:58))
> y <- array(NA,dim=c(5,58),dimnames=list(c('Bier','Tarwe','Suiker','Mineraalwater','Fruit'),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
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
Bier Tarwe Suiker Mineraalwater Fruit
1 102.89 167.16 100.70 106.88 97.69
2 102.64 179.84 99.62 107.45 101.69
3 103.33 174.44 99.83 107.65 102.72
4 103.56 180.35 100.74 107.72 101.85
5 103.60 193.17 100.84 108.10 114.94
6 104.24 195.16 100.85 108.38 106.20
7 105.31 202.43 99.71 108.62 106.76
8 105.40 189.91 100.80 108.79 107.24
9 105.89 195.98 100.06 109.03 106.50
10 105.89 212.09 100.57 109.34 106.77
11 105.54 205.81 99.79 109.73 108.24
12 106.15 204.31 99.90 109.76 104.43
13 106.14 196.07 100.12 109.96 100.90
14 105.85 199.98 100.40 110.49 103.91
15 106.27 199.10 100.51 111.37 103.81
16 106.51 198.31 100.70 111.56 104.59
17 106.82 195.72 100.62 111.90 104.94
18 106.53 223.04 99.70 111.96 111.64
19 107.14 238.41 99.48 112.25 111.27
20 107.39 259.73 99.36 112.39 106.82
21 107.33 326.54 99.39 112.30 106.07
22 107.53 335.15 99.45 112.49 111.35
23 107.42 321.81 99.28 112.77 112.59
24 108.25 368.62 99.40 113.15 108.59
25 108.26 369.59 99.10 113.15 106.83
26 108.93 425.00 99.48 113.28 112.51
27 109.43 439.72 99.74 113.83 113.61
28 109.61 362.23 100.42 114.49 114.96
29 109.74 328.76 100.80 114.76 118.66
30 110.12 348.55 100.66 114.96 116.84
31 110.16 328.18 101.03 115.41 121.19
32 110.44 329.34 101.22 115.84 117.42
33 111.23 295.55 101.23 116.31 116.88
34 112.86 237.38 100.10 117.23 115.01
35 112.77 226.85 99.98 117.97 111.81
36 113.04 220.14 99.91 118.08 110.61
37 112.79 239.36 99.84 118.27 110.67
38 113.87 224.69 99.68 118.88 113.28
39 114.28 230.98 99.74 119.11 112.08
40 115.51 233.47 99.71 119.29 111.41
41 116.76 256.70 99.35 119.36 113.81
42 116.91 253.41 99.21 119.48 109.16
43 116.47 224.95 99.21 120.10 105.09
44 116.94 210.37 99.16 120.30 102.23
45 117.24 191.09 99.20 120.54 101.95
46 116.82 198.85 99.08 120.86 104.75
47 117.48 211.04 98.16 121.10 107.25
48 117.11 206.25 98.00 121.42 105.25
49 117.31 201.19 97.90 121.81 102.75
50 117.77 194.37 97.88 122.21 107.21
51 118.37 191.08 97.56 122.82 107.24
52 117.91 192.87 96.86 123.02 106.01
53 118.12 181.61 96.86 123.14 121.36
54 118.02 157.67 96.75 123.12 120.44
55 117.77 196.14 97.12 123.42 109.40
56 117.85 246.35 97.22 123.50 111.51
57 118.68 271.90 97.52 125.77 111.97
58 118.90 270.29 97.57 125.99 114.64
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Tarwe Suiker Mineraalwater Fruit
-17.210137 -0.001912 0.187982 1.000666 -0.052018
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.9514 -0.4750 -0.0747 0.5220 2.2655
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.210137 16.756988 -1.027 0.309
Tarwe -0.001912 0.001873 -1.021 0.312
Suiker 0.187982 0.143322 1.312 0.195
Mineraalwater 1.000666 0.033099 30.233 <2e-16 ***
Fruit -0.052018 0.024704 -2.106 0.040 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8126 on 53 degrees of freedom
Multiple R-squared: 0.9779, Adjusted R-squared: 0.9762
F-statistic: 586.4 on 4 and 53 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,] 0.212088462 0.4241769248 0.7879115376
[2,] 0.108952444 0.2179048885 0.8910475557
[3,] 0.064027024 0.1280540484 0.9359729758
[4,] 0.128407432 0.2568148644 0.8715925678
[5,] 0.080016837 0.1600336748 0.9199831626
[6,] 0.055648807 0.1112976131 0.9443511934
[7,] 0.077686550 0.1553730991 0.9223134504
[8,] 0.069480825 0.1389616490 0.9305191755
[9,] 0.045571408 0.0911428168 0.9544285916
[10,] 0.032098882 0.0641977638 0.9679011181
[11,] 0.059997163 0.1199943269 0.9400028366
[12,] 0.070033796 0.1400675929 0.9299662035
[13,] 0.096566765 0.1931335307 0.9034332346
[14,] 0.101016726 0.2020334519 0.8989832740
[15,] 0.073800027 0.1476000535 0.9261999733
[16,] 0.070739907 0.1414798134 0.9292600933
[17,] 0.054711106 0.1094222121 0.9452888939
[18,] 0.057299290 0.1145985806 0.9427007097
[19,] 0.041951611 0.0839032224 0.9580483888
[20,] 0.028633494 0.0572669873 0.9713665063
[21,] 0.017898938 0.0357978761 0.9821010620
[22,] 0.010857376 0.0217147518 0.9891426241
[23,] 0.006515311 0.0130306211 0.9934846894
[24,] 0.003614974 0.0072299472 0.9963850264
[25,] 0.002284027 0.0045680541 0.9977159730
[26,] 0.001639750 0.0032794994 0.9983602503
[27,] 0.005835618 0.0116712361 0.9941643820
[28,] 0.007552895 0.0151057903 0.9924471049
[29,] 0.013434885 0.0268697708 0.9865651146
[30,] 0.101645935 0.2032918691 0.8983540655
[31,] 0.326870521 0.6537410428 0.6731294786
[32,] 0.920083605 0.1598327906 0.0799163953
[33,] 0.997170357 0.0056592855 0.0028296428
[34,] 0.999031927 0.0019361453 0.0009680727
[35,] 0.999669876 0.0006602477 0.0003301239
[36,] 0.999365174 0.0012696526 0.0006348263
[37,] 0.998226810 0.0035463793 0.0017731896
[38,] 0.996103017 0.0077939651 0.0038969825
[39,] 0.996230099 0.0075398013 0.0037699006
[40,] 0.993150537 0.0136989265 0.0068494633
[41,] 0.983774543 0.0324509136 0.0162254568
[42,] 0.973350311 0.0532993790 0.0266496895
[43,] 0.956691453 0.0866170942 0.0433085471
> postscript(file="/var/www/html/rcomp/tmp/1kk4b1290543700.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/2dtle1290543700.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/3dtle1290543700.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/4dtle1290543700.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/552kh1290543700.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 6
-0.37967567 -0.76472326 -0.27107663 -0.31614537 0.03022688 -0.06267372
7 8 9 10 11 12
1.02449346 0.74051531 1.10257217 0.74133523 0.21216346 0.57040904
13 14 15 16 17 18
0.11954436 -0.58939489 -1.07754312 -1.02432232 -1.02625488 -0.80260543
19 20 21 22 23 24
-0.43130821 -0.48956929 -0.37644960 -0.08674073 -0.40596831 -0.09737061
25 26 27 28 29 30
-0.12067359 0.74919022 0.73530690 0.04913534 -0.03399137 0.11535015
31 32 33 34 35 36
-0.17716314 -0.55705705 -0.33193208 0.38140439 -0.61311759 -0.51528056
37 38 39 40 41 42
-0.90238668 -0.29499148 -0.17682151 0.84860576 2.26548224 2.07354648
43 44 45 46 47 48
0.74701615 0.84963936 0.85053983 0.29336914 1.03950037 0.26617166
49 50 51 52 53 54
-0.04500780 0.23744920 0.28246861 -0.30663762 0.56023609 0.40730752
55 56 57 58
-0.71318715 -0.52629992 -1.95143740 -1.82517234
> postscript(file="/var/www/html/rcomp/tmp/652kh1290543700.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 -0.37967567 NA
1 -0.76472326 -0.37967567
2 -0.27107663 -0.76472326
3 -0.31614537 -0.27107663
4 0.03022688 -0.31614537
5 -0.06267372 0.03022688
6 1.02449346 -0.06267372
7 0.74051531 1.02449346
8 1.10257217 0.74051531
9 0.74133523 1.10257217
10 0.21216346 0.74133523
11 0.57040904 0.21216346
12 0.11954436 0.57040904
13 -0.58939489 0.11954436
14 -1.07754312 -0.58939489
15 -1.02432232 -1.07754312
16 -1.02625488 -1.02432232
17 -0.80260543 -1.02625488
18 -0.43130821 -0.80260543
19 -0.48956929 -0.43130821
20 -0.37644960 -0.48956929
21 -0.08674073 -0.37644960
22 -0.40596831 -0.08674073
23 -0.09737061 -0.40596831
24 -0.12067359 -0.09737061
25 0.74919022 -0.12067359
26 0.73530690 0.74919022
27 0.04913534 0.73530690
28 -0.03399137 0.04913534
29 0.11535015 -0.03399137
30 -0.17716314 0.11535015
31 -0.55705705 -0.17716314
32 -0.33193208 -0.55705705
33 0.38140439 -0.33193208
34 -0.61311759 0.38140439
35 -0.51528056 -0.61311759
36 -0.90238668 -0.51528056
37 -0.29499148 -0.90238668
38 -0.17682151 -0.29499148
39 0.84860576 -0.17682151
40 2.26548224 0.84860576
41 2.07354648 2.26548224
42 0.74701615 2.07354648
43 0.84963936 0.74701615
44 0.85053983 0.84963936
45 0.29336914 0.85053983
46 1.03950037 0.29336914
47 0.26617166 1.03950037
48 -0.04500780 0.26617166
49 0.23744920 -0.04500780
50 0.28246861 0.23744920
51 -0.30663762 0.28246861
52 0.56023609 -0.30663762
53 0.40730752 0.56023609
54 -0.71318715 0.40730752
55 -0.52629992 -0.71318715
56 -1.95143740 -0.52629992
57 -1.82517234 -1.95143740
58 NA -1.82517234
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.76472326 -0.37967567
[2,] -0.27107663 -0.76472326
[3,] -0.31614537 -0.27107663
[4,] 0.03022688 -0.31614537
[5,] -0.06267372 0.03022688
[6,] 1.02449346 -0.06267372
[7,] 0.74051531 1.02449346
[8,] 1.10257217 0.74051531
[9,] 0.74133523 1.10257217
[10,] 0.21216346 0.74133523
[11,] 0.57040904 0.21216346
[12,] 0.11954436 0.57040904
[13,] -0.58939489 0.11954436
[14,] -1.07754312 -0.58939489
[15,] -1.02432232 -1.07754312
[16,] -1.02625488 -1.02432232
[17,] -0.80260543 -1.02625488
[18,] -0.43130821 -0.80260543
[19,] -0.48956929 -0.43130821
[20,] -0.37644960 -0.48956929
[21,] -0.08674073 -0.37644960
[22,] -0.40596831 -0.08674073
[23,] -0.09737061 -0.40596831
[24,] -0.12067359 -0.09737061
[25,] 0.74919022 -0.12067359
[26,] 0.73530690 0.74919022
[27,] 0.04913534 0.73530690
[28,] -0.03399137 0.04913534
[29,] 0.11535015 -0.03399137
[30,] -0.17716314 0.11535015
[31,] -0.55705705 -0.17716314
[32,] -0.33193208 -0.55705705
[33,] 0.38140439 -0.33193208
[34,] -0.61311759 0.38140439
[35,] -0.51528056 -0.61311759
[36,] -0.90238668 -0.51528056
[37,] -0.29499148 -0.90238668
[38,] -0.17682151 -0.29499148
[39,] 0.84860576 -0.17682151
[40,] 2.26548224 0.84860576
[41,] 2.07354648 2.26548224
[42,] 0.74701615 2.07354648
[43,] 0.84963936 0.74701615
[44,] 0.85053983 0.84963936
[45,] 0.29336914 0.85053983
[46,] 1.03950037 0.29336914
[47,] 0.26617166 1.03950037
[48,] -0.04500780 0.26617166
[49,] 0.23744920 -0.04500780
[50,] 0.28246861 0.23744920
[51,] -0.30663762 0.28246861
[52,] 0.56023609 -0.30663762
[53,] 0.40730752 0.56023609
[54,] -0.71318715 0.40730752
[55,] -0.52629992 -0.71318715
[56,] -1.95143740 -0.52629992
[57,] -1.82517234 -1.95143740
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.76472326 -0.37967567
2 -0.27107663 -0.76472326
3 -0.31614537 -0.27107663
4 0.03022688 -0.31614537
5 -0.06267372 0.03022688
6 1.02449346 -0.06267372
7 0.74051531 1.02449346
8 1.10257217 0.74051531
9 0.74133523 1.10257217
10 0.21216346 0.74133523
11 0.57040904 0.21216346
12 0.11954436 0.57040904
13 -0.58939489 0.11954436
14 -1.07754312 -0.58939489
15 -1.02432232 -1.07754312
16 -1.02625488 -1.02432232
17 -0.80260543 -1.02625488
18 -0.43130821 -0.80260543
19 -0.48956929 -0.43130821
20 -0.37644960 -0.48956929
21 -0.08674073 -0.37644960
22 -0.40596831 -0.08674073
23 -0.09737061 -0.40596831
24 -0.12067359 -0.09737061
25 0.74919022 -0.12067359
26 0.73530690 0.74919022
27 0.04913534 0.73530690
28 -0.03399137 0.04913534
29 0.11535015 -0.03399137
30 -0.17716314 0.11535015
31 -0.55705705 -0.17716314
32 -0.33193208 -0.55705705
33 0.38140439 -0.33193208
34 -0.61311759 0.38140439
35 -0.51528056 -0.61311759
36 -0.90238668 -0.51528056
37 -0.29499148 -0.90238668
38 -0.17682151 -0.29499148
39 0.84860576 -0.17682151
40 2.26548224 0.84860576
41 2.07354648 2.26548224
42 0.74701615 2.07354648
43 0.84963936 0.74701615
44 0.85053983 0.84963936
45 0.29336914 0.85053983
46 1.03950037 0.29336914
47 0.26617166 1.03950037
48 -0.04500780 0.26617166
49 0.23744920 -0.04500780
50 0.28246861 0.23744920
51 -0.30663762 0.28246861
52 0.56023609 -0.30663762
53 0.40730752 0.56023609
54 -0.71318715 0.40730752
55 -0.52629992 -0.71318715
56 -1.95143740 -0.52629992
57 -1.82517234 -1.95143740
> 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/7yck21290543700.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/8yck21290543700.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/9eoqt1290543700.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/10eoqt1290543700.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/11ndhd1290543700.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/12qdfj1290543700.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/13m5da1290543700.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/148ncy1290543700.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/15toa41290543700.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/16e79a1290543700.tab")
+ }
>
> try(system("convert tmp/1kk4b1290543700.ps tmp/1kk4b1290543700.png",intern=TRUE))
character(0)
> try(system("convert tmp/2dtle1290543700.ps tmp/2dtle1290543700.png",intern=TRUE))
character(0)
> try(system("convert tmp/3dtle1290543700.ps tmp/3dtle1290543700.png",intern=TRUE))
character(0)
> try(system("convert tmp/4dtle1290543700.ps tmp/4dtle1290543700.png",intern=TRUE))
character(0)
> try(system("convert tmp/552kh1290543700.ps tmp/552kh1290543700.png",intern=TRUE))
character(0)
> try(system("convert tmp/652kh1290543700.ps tmp/652kh1290543700.png",intern=TRUE))
character(0)
> try(system("convert tmp/7yck21290543700.ps tmp/7yck21290543700.png",intern=TRUE))
character(0)
> try(system("convert tmp/8yck21290543700.ps tmp/8yck21290543700.png",intern=TRUE))
character(0)
> try(system("convert tmp/9eoqt1290543700.ps tmp/9eoqt1290543700.png",intern=TRUE))
character(0)
> try(system("convert tmp/10eoqt1290543700.ps tmp/10eoqt1290543700.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.436 1.571 5.907