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(1
+ ,102.89
+ ,167.16
+ ,100.70
+ ,106.88
+ ,97.69
+ ,1
+ ,106.14
+ ,196.07
+ ,100.12
+ ,109.96
+ ,100.90
+ ,1
+ ,108.26
+ ,369.59
+ ,99.10
+ ,113.15
+ ,106.83
+ ,1
+ ,112.79
+ ,239.36
+ ,99.84
+ ,118.27
+ ,110.67
+ ,1
+ ,117.31
+ ,201.19
+ ,97.90
+ ,121.81
+ ,102.75
+ ,2
+ ,102.64
+ ,179.84
+ ,99.62
+ ,107.45
+ ,101.69
+ ,2
+ ,105.85
+ ,199.98
+ ,100.40
+ ,110.49
+ ,103.91
+ ,2
+ ,108.93
+ ,425.00
+ ,99.48
+ ,113.28
+ ,112.51
+ ,2
+ ,113.87
+ ,224.69
+ ,99.68
+ ,118.88
+ ,113.28
+ ,2
+ ,117.77
+ ,194.37
+ ,97.88
+ ,122.21
+ ,107.21
+ ,3
+ ,103.33
+ ,174.44
+ ,99.83
+ ,107.65
+ ,102.72
+ ,3
+ ,106.27
+ ,199.10
+ ,100.51
+ ,111.37
+ ,103.81
+ ,3
+ ,109.43
+ ,439.72
+ ,99.74
+ ,113.83
+ ,113.61
+ ,3
+ ,114.28
+ ,230.98
+ ,99.74
+ ,119.11
+ ,112.08
+ ,3
+ ,118.37
+ ,191.08
+ ,97.56
+ ,122.82
+ ,107.24
+ ,4
+ ,103.56
+ ,180.35
+ ,100.74
+ ,107.72
+ ,101.85
+ ,4
+ ,106.51
+ ,198.31
+ ,100.70
+ ,111.56
+ ,104.59
+ ,4
+ ,109.61
+ ,362.23
+ ,100.42
+ ,114.49
+ ,114.96
+ ,4
+ ,115.51
+ ,233.47
+ ,99.71
+ ,119.29
+ ,111.41
+ ,4
+ ,117.91
+ ,192.87
+ ,96.86
+ ,123.02
+ ,106.01
+ ,5
+ ,103.60
+ ,193.17
+ ,100.84
+ ,108.10
+ ,114.94
+ ,5
+ ,106.82
+ ,195.72
+ ,100.62
+ ,111.90
+ ,104.94
+ ,5
+ ,109.74
+ ,328.76
+ ,100.80
+ ,114.76
+ ,118.66
+ ,5
+ ,116.76
+ ,256.70
+ ,99.35
+ ,119.36
+ ,113.81
+ ,5
+ ,118.12
+ ,181.61
+ ,96.86
+ ,123.14
+ ,121.36
+ ,6
+ ,104.24
+ ,195.16
+ ,100.85
+ ,108.38
+ ,106.20
+ ,6
+ ,106.53
+ ,223.04
+ ,99.70
+ ,111.96
+ ,111.64
+ ,6
+ ,110.12
+ ,348.55
+ ,100.66
+ ,114.96
+ ,116.84
+ ,6
+ ,116.91
+ ,253.41
+ ,99.21
+ ,119.48
+ ,109.16
+ ,6
+ ,118.02
+ ,157.67
+ ,96.75
+ ,123.12
+ ,120.44
+ ,7
+ ,105.31
+ ,202.43
+ ,99.71
+ ,108.62
+ ,106.76
+ ,7
+ ,107.14
+ ,238.41
+ ,99.48
+ ,112.25
+ ,111.27
+ ,7
+ ,110.16
+ ,328.18
+ ,101.03
+ ,115.41
+ ,121.19
+ ,7
+ ,116.47
+ ,224.95
+ ,99.21
+ ,120.10
+ ,105.09
+ ,7
+ ,117.77
+ ,196.14
+ ,97.12
+ ,123.42
+ ,109.40
+ ,8
+ ,105.40
+ ,189.91
+ ,100.80
+ ,108.79
+ ,107.24
+ ,8
+ ,107.39
+ ,259.73
+ ,99.36
+ ,112.39
+ ,106.82
+ ,8
+ ,110.44
+ ,329.34
+ ,101.22
+ ,115.84
+ ,117.42
+ ,8
+ ,116.94
+ ,210.37
+ ,99.16
+ ,120.30
+ ,102.23
+ ,8
+ ,117.85
+ ,246.35
+ ,97.22
+ ,123.50
+ ,111.51
+ ,9
+ ,105.89
+ ,195.98
+ ,100.06
+ ,109.03
+ ,106.50
+ ,9
+ ,107.33
+ ,326.54
+ ,99.39
+ ,112.30
+ ,106.07
+ ,9
+ ,111.23
+ ,295.55
+ ,101.23
+ ,116.31
+ ,116.88
+ ,9
+ ,117.24
+ ,191.09
+ ,99.20
+ ,120.54
+ ,101.95
+ ,9
+ ,118.68
+ ,271.90
+ ,97.52
+ ,125.77
+ ,111.97
+ ,10
+ ,105.89
+ ,212.09
+ ,100.57
+ ,109.34
+ ,106.77
+ ,10
+ ,107.53
+ ,335.15
+ ,99.45
+ ,112.49
+ ,111.35
+ ,10
+ ,112.86
+ ,237.38
+ ,100.10
+ ,117.23
+ ,115.01
+ ,10
+ ,116.82
+ ,198.85
+ ,99.08
+ ,120.86
+ ,104.75
+ ,10
+ ,118.90
+ ,270.29
+ ,97.57
+ ,125.99
+ ,114.64
+ ,11
+ ,105.54
+ ,205.81
+ ,99.79
+ ,109.73
+ ,108.24
+ ,11
+ ,107.42
+ ,321.81
+ ,99.28
+ ,112.77
+ ,112.59
+ ,11
+ ,112.77
+ ,226.85
+ ,99.98
+ ,117.97
+ ,111.81
+ ,11
+ ,117.48
+ ,211.04
+ ,98.16
+ ,121.10
+ ,107.25
+ ,12
+ ,106.15
+ ,204.31
+ ,99.90
+ ,109.76
+ ,104.43
+ ,12
+ ,108.25
+ ,368.62
+ ,99.40
+ ,113.15
+ ,108.59
+ ,12
+ ,113.04
+ ,220.14
+ ,99.91
+ ,118.08
+ ,110.61
+ ,12
+ ,117.11
+ ,206.25
+ ,98.00
+ ,121.42
+ ,105.25)
+ ,dim=c(6
+ ,58)
+ ,dimnames=list(c('month'
+ ,'Bier'
+ ,'Tarwe'
+ ,'suiker'
+ ,'minerwater'
+ ,'fruit')
+ ,1:58))
> y <- array(NA,dim=c(6,58),dimnames=list(c('month','Bier','Tarwe','suiker','minerwater','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 = '4'
> #'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
suiker month Bier Tarwe minerwater fruit
1 100.70 1 102.89 167.16 106.88 97.69
2 100.12 1 106.14 196.07 109.96 100.90
3 99.10 1 108.26 369.59 113.15 106.83
4 99.84 1 112.79 239.36 118.27 110.67
5 97.90 1 117.31 201.19 121.81 102.75
6 99.62 2 102.64 179.84 107.45 101.69
7 100.40 2 105.85 199.98 110.49 103.91
8 99.48 2 108.93 425.00 113.28 112.51
9 99.68 2 113.87 224.69 118.88 113.28
10 97.88 2 117.77 194.37 122.21 107.21
11 99.83 3 103.33 174.44 107.65 102.72
12 100.51 3 106.27 199.10 111.37 103.81
13 99.74 3 109.43 439.72 113.83 113.61
14 99.74 3 114.28 230.98 119.11 112.08
15 97.56 3 118.37 191.08 122.82 107.24
16 100.74 4 103.56 180.35 107.72 101.85
17 100.70 4 106.51 198.31 111.56 104.59
18 100.42 4 109.61 362.23 114.49 114.96
19 99.71 4 115.51 233.47 119.29 111.41
20 96.86 4 117.91 192.87 123.02 106.01
21 100.84 5 103.60 193.17 108.10 114.94
22 100.62 5 106.82 195.72 111.90 104.94
23 100.80 5 109.74 328.76 114.76 118.66
24 99.35 5 116.76 256.70 119.36 113.81
25 96.86 5 118.12 181.61 123.14 121.36
26 100.85 6 104.24 195.16 108.38 106.20
27 99.70 6 106.53 223.04 111.96 111.64
28 100.66 6 110.12 348.55 114.96 116.84
29 99.21 6 116.91 253.41 119.48 109.16
30 96.75 6 118.02 157.67 123.12 120.44
31 99.71 7 105.31 202.43 108.62 106.76
32 99.48 7 107.14 238.41 112.25 111.27
33 101.03 7 110.16 328.18 115.41 121.19
34 99.21 7 116.47 224.95 120.10 105.09
35 97.12 7 117.77 196.14 123.42 109.40
36 100.80 8 105.40 189.91 108.79 107.24
37 99.36 8 107.39 259.73 112.39 106.82
38 101.22 8 110.44 329.34 115.84 117.42
39 99.16 8 116.94 210.37 120.30 102.23
40 97.22 8 117.85 246.35 123.50 111.51
41 100.06 9 105.89 195.98 109.03 106.50
42 99.39 9 107.33 326.54 112.30 106.07
43 101.23 9 111.23 295.55 116.31 116.88
44 99.20 9 117.24 191.09 120.54 101.95
45 97.52 9 118.68 271.90 125.77 111.97
46 100.57 10 105.89 212.09 109.34 106.77
47 99.45 10 107.53 335.15 112.49 111.35
48 100.10 10 112.86 237.38 117.23 115.01
49 99.08 10 116.82 198.85 120.86 104.75
50 97.57 10 118.90 270.29 125.99 114.64
51 99.79 11 105.54 205.81 109.73 108.24
52 99.28 11 107.42 321.81 112.77 112.59
53 99.98 11 112.77 226.85 117.97 111.81
54 98.16 11 117.48 211.04 121.10 107.25
55 99.90 12 106.15 204.31 109.76 104.43
56 99.40 12 108.25 368.62 113.15 108.59
57 99.91 12 113.04 220.14 118.08 110.61
58 98.00 12 117.11 206.25 121.42 105.25
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) month Bier Tarwe minerwater fruit
114.597308 0.001065 0.166951 0.001865 -0.338281 0.045539
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.68678 -0.60293 0.01033 0.67601 1.52496
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 114.597308 2.694921 42.523 <2e-16 ***
month 0.001065 0.030571 0.035 0.9723
Bier 0.166951 0.128996 1.294 0.2013
Tarwe 0.001865 0.001782 1.046 0.3004
minerwater -0.338281 0.126346 -2.677 0.0099 **
fruit 0.045539 0.023818 1.912 0.0614 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7738 on 52 degrees of freedom
Multiple R-squared: 0.6208, Adjusted R-squared: 0.5844
F-statistic: 17.03 on 5 and 52 DF, p-value: 6.063e-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.21749151 0.4349830 0.7825085
[2,] 0.11298021 0.2259604 0.8870198
[3,] 0.05181876 0.1036375 0.9481812
[4,] 0.24171052 0.4834210 0.7582895
[5,] 0.20942707 0.4188541 0.7905729
[6,] 0.16899859 0.3379972 0.8310014
[7,] 0.12796130 0.2559226 0.8720387
[8,] 0.11309032 0.2261806 0.8869097
[9,] 0.11739129 0.2347826 0.8826087
[10,] 0.07844484 0.1568897 0.9215552
[11,] 0.06810148 0.1362030 0.9318985
[12,] 0.09454663 0.1890933 0.9054534
[13,] 0.17162971 0.3432594 0.8283703
[14,] 0.20515040 0.4103008 0.7948496
[15,] 0.16256962 0.3251392 0.8374304
[16,] 0.11994960 0.2398992 0.8800504
[17,] 0.52280686 0.9543863 0.4771931
[18,] 0.49279377 0.9855875 0.5072062
[19,] 0.44696334 0.8939267 0.5530367
[20,] 0.39774496 0.7954899 0.6022550
[21,] 0.35211240 0.7042248 0.6478876
[22,] 0.81937288 0.3612542 0.1806271
[23,] 0.85180946 0.2963811 0.1481905
[24,] 0.81339375 0.3732125 0.1866062
[25,] 0.80409924 0.3918015 0.1959008
[26,] 0.75813643 0.4837271 0.2418636
[27,] 0.78895819 0.4220836 0.2110418
[28,] 0.71876447 0.5624711 0.2812355
[29,] 0.65151721 0.6969656 0.3484828
[30,] 0.72652639 0.5469472 0.2734736
[31,] 0.69280329 0.6143934 0.3071967
[32,] 0.82421351 0.3515730 0.1757865
[33,] 0.83241698 0.3351660 0.1675830
[34,] 0.77267920 0.4546416 0.2273208
[35,] 0.86872629 0.2625474 0.1312737
[36,] 0.85633884 0.2873223 0.1436612
[37,] 0.77325339 0.4534932 0.2267466
[38,] 0.68112791 0.6377442 0.3188721
[39,] 0.56675379 0.8664924 0.4332462
[40,] 0.43844730 0.8768946 0.5615527
[41,] 0.79570543 0.4085891 0.2042946
> postscript(file="/var/www/rcomp/tmp/1x2rn1290602748.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/2psqq1290602748.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/3psqq1290602748.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/4psqq1290602748.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/50kpt1290602748.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.319149902 0.038381505 -0.850017867 0.933640854 -0.131618689 -0.733157490
7 8 9 10 11 12
0.400656848 -0.900940869 0.707115066 -0.284558745 -0.608599570 0.743354592
13 14 15 16 17 18
-0.616965895 0.818324371 -0.494674649 0.314216333 0.922446869 0.338186933
19 20 21 22 23 24
0.668669424 -0.998610267 -0.084995150 0.893532906 0.680663546 -0.030011899
25 26 27 28 29 30
-1.672177993 0.306114800 -0.314873096 0.589796362 0.062366992 -1.686780743
31 32 33 34 35 36
-0.971457182 -0.551484673 0.944163220 0.582904415 -0.743596134 0.161444752
37 38 39 40 41 42
-0.504029012 1.401333823 0.678456393 -0.820660871 -0.557857147 -0.585936599
43 44 45 46 47 48
1.524963399 0.797192579 0.039016526 0.013611712 -0.752620523 0.626605624
49 50 51 52 53 54
0.712517330 0.007055519 -0.632324974 -0.842198656 0.936253849 -0.374124992
55 56 57 58
-0.438779309 -0.638406284 0.924481284 -0.265158452
> postscript(file="/var/www/rcomp/tmp/60kpt1290602748.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.319149902 NA
1 0.038381505 0.319149902
2 -0.850017867 0.038381505
3 0.933640854 -0.850017867
4 -0.131618689 0.933640854
5 -0.733157490 -0.131618689
6 0.400656848 -0.733157490
7 -0.900940869 0.400656848
8 0.707115066 -0.900940869
9 -0.284558745 0.707115066
10 -0.608599570 -0.284558745
11 0.743354592 -0.608599570
12 -0.616965895 0.743354592
13 0.818324371 -0.616965895
14 -0.494674649 0.818324371
15 0.314216333 -0.494674649
16 0.922446869 0.314216333
17 0.338186933 0.922446869
18 0.668669424 0.338186933
19 -0.998610267 0.668669424
20 -0.084995150 -0.998610267
21 0.893532906 -0.084995150
22 0.680663546 0.893532906
23 -0.030011899 0.680663546
24 -1.672177993 -0.030011899
25 0.306114800 -1.672177993
26 -0.314873096 0.306114800
27 0.589796362 -0.314873096
28 0.062366992 0.589796362
29 -1.686780743 0.062366992
30 -0.971457182 -1.686780743
31 -0.551484673 -0.971457182
32 0.944163220 -0.551484673
33 0.582904415 0.944163220
34 -0.743596134 0.582904415
35 0.161444752 -0.743596134
36 -0.504029012 0.161444752
37 1.401333823 -0.504029012
38 0.678456393 1.401333823
39 -0.820660871 0.678456393
40 -0.557857147 -0.820660871
41 -0.585936599 -0.557857147
42 1.524963399 -0.585936599
43 0.797192579 1.524963399
44 0.039016526 0.797192579
45 0.013611712 0.039016526
46 -0.752620523 0.013611712
47 0.626605624 -0.752620523
48 0.712517330 0.626605624
49 0.007055519 0.712517330
50 -0.632324974 0.007055519
51 -0.842198656 -0.632324974
52 0.936253849 -0.842198656
53 -0.374124992 0.936253849
54 -0.438779309 -0.374124992
55 -0.638406284 -0.438779309
56 0.924481284 -0.638406284
57 -0.265158452 0.924481284
58 NA -0.265158452
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.038381505 0.319149902
[2,] -0.850017867 0.038381505
[3,] 0.933640854 -0.850017867
[4,] -0.131618689 0.933640854
[5,] -0.733157490 -0.131618689
[6,] 0.400656848 -0.733157490
[7,] -0.900940869 0.400656848
[8,] 0.707115066 -0.900940869
[9,] -0.284558745 0.707115066
[10,] -0.608599570 -0.284558745
[11,] 0.743354592 -0.608599570
[12,] -0.616965895 0.743354592
[13,] 0.818324371 -0.616965895
[14,] -0.494674649 0.818324371
[15,] 0.314216333 -0.494674649
[16,] 0.922446869 0.314216333
[17,] 0.338186933 0.922446869
[18,] 0.668669424 0.338186933
[19,] -0.998610267 0.668669424
[20,] -0.084995150 -0.998610267
[21,] 0.893532906 -0.084995150
[22,] 0.680663546 0.893532906
[23,] -0.030011899 0.680663546
[24,] -1.672177993 -0.030011899
[25,] 0.306114800 -1.672177993
[26,] -0.314873096 0.306114800
[27,] 0.589796362 -0.314873096
[28,] 0.062366992 0.589796362
[29,] -1.686780743 0.062366992
[30,] -0.971457182 -1.686780743
[31,] -0.551484673 -0.971457182
[32,] 0.944163220 -0.551484673
[33,] 0.582904415 0.944163220
[34,] -0.743596134 0.582904415
[35,] 0.161444752 -0.743596134
[36,] -0.504029012 0.161444752
[37,] 1.401333823 -0.504029012
[38,] 0.678456393 1.401333823
[39,] -0.820660871 0.678456393
[40,] -0.557857147 -0.820660871
[41,] -0.585936599 -0.557857147
[42,] 1.524963399 -0.585936599
[43,] 0.797192579 1.524963399
[44,] 0.039016526 0.797192579
[45,] 0.013611712 0.039016526
[46,] -0.752620523 0.013611712
[47,] 0.626605624 -0.752620523
[48,] 0.712517330 0.626605624
[49,] 0.007055519 0.712517330
[50,] -0.632324974 0.007055519
[51,] -0.842198656 -0.632324974
[52,] 0.936253849 -0.842198656
[53,] -0.374124992 0.936253849
[54,] -0.438779309 -0.374124992
[55,] -0.638406284 -0.438779309
[56,] 0.924481284 -0.638406284
[57,] -0.265158452 0.924481284
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.038381505 0.319149902
2 -0.850017867 0.038381505
3 0.933640854 -0.850017867
4 -0.131618689 0.933640854
5 -0.733157490 -0.131618689
6 0.400656848 -0.733157490
7 -0.900940869 0.400656848
8 0.707115066 -0.900940869
9 -0.284558745 0.707115066
10 -0.608599570 -0.284558745
11 0.743354592 -0.608599570
12 -0.616965895 0.743354592
13 0.818324371 -0.616965895
14 -0.494674649 0.818324371
15 0.314216333 -0.494674649
16 0.922446869 0.314216333
17 0.338186933 0.922446869
18 0.668669424 0.338186933
19 -0.998610267 0.668669424
20 -0.084995150 -0.998610267
21 0.893532906 -0.084995150
22 0.680663546 0.893532906
23 -0.030011899 0.680663546
24 -1.672177993 -0.030011899
25 0.306114800 -1.672177993
26 -0.314873096 0.306114800
27 0.589796362 -0.314873096
28 0.062366992 0.589796362
29 -1.686780743 0.062366992
30 -0.971457182 -1.686780743
31 -0.551484673 -0.971457182
32 0.944163220 -0.551484673
33 0.582904415 0.944163220
34 -0.743596134 0.582904415
35 0.161444752 -0.743596134
36 -0.504029012 0.161444752
37 1.401333823 -0.504029012
38 0.678456393 1.401333823
39 -0.820660871 0.678456393
40 -0.557857147 -0.820660871
41 -0.585936599 -0.557857147
42 1.524963399 -0.585936599
43 0.797192579 1.524963399
44 0.039016526 0.797192579
45 0.013611712 0.039016526
46 -0.752620523 0.013611712
47 0.626605624 -0.752620523
48 0.712517330 0.626605624
49 0.007055519 0.712517330
50 -0.632324974 0.007055519
51 -0.842198656 -0.632324974
52 0.936253849 -0.842198656
53 -0.374124992 0.936253849
54 -0.438779309 -0.374124992
55 -0.638406284 -0.438779309
56 0.924481284 -0.638406284
57 -0.265158452 0.924481284
> 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/7tbow1290602748.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/8tbow1290602748.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/9325h1290602748.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/10325h1290602748.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/117l4n1290602748.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/12sm3t1290602748.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/13h5i41290602748.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/14aehp1290602748.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/15vwxv1290602748.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/16rov41290602748.tab")
+ }
>
> try(system("convert tmp/1x2rn1290602748.ps tmp/1x2rn1290602748.png",intern=TRUE))
character(0)
> try(system("convert tmp/2psqq1290602748.ps tmp/2psqq1290602748.png",intern=TRUE))
character(0)
> try(system("convert tmp/3psqq1290602748.ps tmp/3psqq1290602748.png",intern=TRUE))
character(0)
> try(system("convert tmp/4psqq1290602748.ps tmp/4psqq1290602748.png",intern=TRUE))
character(0)
> try(system("convert tmp/50kpt1290602748.ps tmp/50kpt1290602748.png",intern=TRUE))
character(0)
> try(system("convert tmp/60kpt1290602748.ps tmp/60kpt1290602748.png",intern=TRUE))
character(0)
> try(system("convert tmp/7tbow1290602748.ps tmp/7tbow1290602748.png",intern=TRUE))
character(0)
> try(system("convert tmp/8tbow1290602748.ps tmp/8tbow1290602748.png",intern=TRUE))
character(0)
> try(system("convert tmp/9325h1290602748.ps tmp/9325h1290602748.png",intern=TRUE))
character(0)
> try(system("convert tmp/10325h1290602748.ps tmp/10325h1290602748.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.780 0.960 4.738