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'
+ ,'minerwater'
+ ,'fruit
')
+ ,1:58))
> y <- array(NA,dim=c(5,58),dimnames=list(c('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 = '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 minerwater fruit\r t
1 102.89 167.16 100.70 106.88 97.69 1
2 102.64 179.84 99.62 107.45 101.69 2
3 103.33 174.44 99.83 107.65 102.72 3
4 103.56 180.35 100.74 107.72 101.85 4
5 103.60 193.17 100.84 108.10 114.94 5
6 104.24 195.16 100.85 108.38 106.20 6
7 105.31 202.43 99.71 108.62 106.76 7
8 105.40 189.91 100.80 108.79 107.24 8
9 105.89 195.98 100.06 109.03 106.50 9
10 105.89 212.09 100.57 109.34 106.77 10
11 105.54 205.81 99.79 109.73 108.24 11
12 106.15 204.31 99.90 109.76 104.43 12
13 106.14 196.07 100.12 109.96 100.90 13
14 105.85 199.98 100.40 110.49 103.91 14
15 106.27 199.10 100.51 111.37 103.81 15
16 106.51 198.31 100.70 111.56 104.59 16
17 106.82 195.72 100.62 111.90 104.94 17
18 106.53 223.04 99.70 111.96 111.64 18
19 107.14 238.41 99.48 112.25 111.27 19
20 107.39 259.73 99.36 112.39 106.82 20
21 107.33 326.54 99.39 112.30 106.07 21
22 107.53 335.15 99.45 112.49 111.35 22
23 107.42 321.81 99.28 112.77 112.59 23
24 108.25 368.62 99.40 113.15 108.59 24
25 108.26 369.59 99.10 113.15 106.83 25
26 108.93 425.00 99.48 113.28 112.51 26
27 109.43 439.72 99.74 113.83 113.61 27
28 109.61 362.23 100.42 114.49 114.96 28
29 109.74 328.76 100.80 114.76 118.66 29
30 110.12 348.55 100.66 114.96 116.84 30
31 110.16 328.18 101.03 115.41 121.19 31
32 110.44 329.34 101.22 115.84 117.42 32
33 111.23 295.55 101.23 116.31 116.88 33
34 112.86 237.38 100.10 117.23 115.01 34
35 112.77 226.85 99.98 117.97 111.81 35
36 113.04 220.14 99.91 118.08 110.61 36
37 112.79 239.36 99.84 118.27 110.67 37
38 113.87 224.69 99.68 118.88 113.28 38
39 114.28 230.98 99.74 119.11 112.08 39
40 115.51 233.47 99.71 119.29 111.41 40
41 116.76 256.70 99.35 119.36 113.81 41
42 116.91 253.41 99.21 119.48 109.16 42
43 116.47 224.95 99.21 120.10 105.09 43
44 116.94 210.37 99.16 120.30 102.23 44
45 117.24 191.09 99.20 120.54 101.95 45
46 116.82 198.85 99.08 120.86 104.75 46
47 117.48 211.04 98.16 121.10 107.25 47
48 117.11 206.25 98.00 121.42 105.25 48
49 117.31 201.19 97.90 121.81 102.75 49
50 117.77 194.37 97.88 122.21 107.21 50
51 118.37 191.08 97.56 122.82 107.24 51
52 117.91 192.87 96.86 123.02 106.01 52
53 118.12 181.61 96.86 123.14 121.36 53
54 118.02 157.67 96.75 123.12 120.44 54
55 117.77 196.14 97.12 123.42 109.40 55
56 117.85 246.35 97.22 123.50 111.51 56
57 118.68 271.90 97.52 125.77 111.97 57
58 118.90 270.29 97.57 125.99 114.64 58
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Tarwe suiker minerwater `fruit\r` t
69.863374 -0.006834 0.297906 0.089566 -0.055912 0.301836
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.2239 -0.4158 -0.1634 0.3889 2.3515
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 69.863374 28.616018 2.441 0.01807 *
Tarwe -0.006834 0.002180 -3.136 0.00282 **
suiker 0.297906 0.133131 2.238 0.02955 *
minerwater 0.089566 0.255794 0.350 0.72764
`fruit\r` -0.055912 0.022358 -2.501 0.01558 *
t 0.301836 0.084159 3.586 0.00074 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7345 on 52 degrees of freedom
Multiple R-squared: 0.9823, Adjusted R-squared: 0.9806
F-statistic: 576.7 on 5 and 52 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.0496990253 0.0993980506 0.9503009747
[2,] 0.0222756821 0.0445513641 0.9777243179
[3,] 0.0564725396 0.1129450792 0.9435274604
[4,] 0.0941182560 0.1882365120 0.9058817440
[5,] 0.0766250365 0.1532500730 0.9233749635
[6,] 0.0399444133 0.0798888265 0.9600555867
[7,] 0.0575262015 0.1150524029 0.9424737985
[8,] 0.0344679144 0.0689358288 0.9655320856
[9,] 0.0228079670 0.0456159340 0.9771920330
[10,] 0.0382155980 0.0764311961 0.9617844020
[11,] 0.0248468396 0.0496936792 0.9751531604
[12,] 0.0166701843 0.0333403686 0.9833298157
[13,] 0.0114028351 0.0228056703 0.9885971649
[14,] 0.0066673800 0.0133347600 0.9933326200
[15,] 0.0078581376 0.0157162753 0.9921418624
[16,] 0.0070421820 0.0140843639 0.9929578180
[17,] 0.0057048391 0.0114096781 0.9942951609
[18,] 0.0040994076 0.0081988152 0.9959005924
[19,] 0.0050964825 0.0101929649 0.9949035175
[20,] 0.0029587394 0.0059174789 0.9970412606
[21,] 0.0015395483 0.0030790967 0.9984604517
[22,] 0.0007502976 0.0015005952 0.9992497024
[23,] 0.0003830091 0.0007660182 0.9996169909
[24,] 0.0002752704 0.0005505408 0.9997247296
[25,] 0.0006316975 0.0012633951 0.9993683025
[26,] 0.0119296006 0.0238592012 0.9880703994
[27,] 0.0134119723 0.0268239445 0.9865880277
[28,] 0.0176538813 0.0353077626 0.9823461187
[29,] 0.1060317323 0.2120634646 0.8939682677
[30,] 0.2755605423 0.5511210845 0.7244394577
[31,] 0.8871378242 0.2257243516 0.1128621758
[32,] 0.9969383899 0.0061232202 0.0030616101
[33,] 0.9986782662 0.0026434676 0.0013217338
[34,] 0.9992330835 0.0015338329 0.0007669165
[35,] 0.9990368861 0.0019262278 0.0009631139
[36,] 0.9970986214 0.0058027572 0.0029013786
[37,] 0.9946794212 0.0106411577 0.0053205788
[38,] 0.9917974679 0.0164050642 0.0082025321
[39,] 0.9819974776 0.0360050448 0.0180025224
[40,] 0.9673475075 0.0653049851 0.0326524925
[41,] 0.9603722621 0.0792554758 0.0396277379
> postscript(file="/var/www/html/rcomp/tmp/1mvpe1290545353.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/2mvpe1290545353.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/3fmoz1290545353.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/4fmoz1290545353.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/5fmoz1290545353.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.24268364 -0.21352657 0.11484774 -0.24260447 0.25123533 0.08627380
7 8 9 10 11 12
1.25355161 0.64304338 1.03027223 0.67393657 0.25880747 0.30823971
13 14 15 16 17 18
-0.34073216 -0.86843481 -0.87346346 -0.97070692 -0.96729444 -0.72910721
19 20 21 22 23 24
-0.29702098 -0.42874601 -0.37678852 -0.15945837 -0.56756829 -0.01291849
25 26 27 28 29 30
-0.30715782 0.63242892 0.86598106 0.02834222 -0.30275317 -0.16730270
31 32 33 34 35 36
-0.47566811 -0.79547875 -0.61351475 0.46677351 -0.20647549 -0.34026274
37 38 39 40 41 42
-0.75355155 0.06331221 0.10889593 1.00943216 2.35152326 1.94817177
43 44 45 46 47 48
0.72873866 0.63433250 0.45166297 -0.05349762 0.78033466 -0.01705724
49 50 51 52 53 54
-0.29839425 0.03265786 0.35070930 -0.27704328 0.40166340 -0.18066476
55 56 57 58
-1.22394348 -1.02160831 -0.58579267 -0.56394721
> postscript(file="/var/www/html/rcomp/tmp/68wo21290545353.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.24268364 NA
1 -0.21352657 -0.24268364
2 0.11484774 -0.21352657
3 -0.24260447 0.11484774
4 0.25123533 -0.24260447
5 0.08627380 0.25123533
6 1.25355161 0.08627380
7 0.64304338 1.25355161
8 1.03027223 0.64304338
9 0.67393657 1.03027223
10 0.25880747 0.67393657
11 0.30823971 0.25880747
12 -0.34073216 0.30823971
13 -0.86843481 -0.34073216
14 -0.87346346 -0.86843481
15 -0.97070692 -0.87346346
16 -0.96729444 -0.97070692
17 -0.72910721 -0.96729444
18 -0.29702098 -0.72910721
19 -0.42874601 -0.29702098
20 -0.37678852 -0.42874601
21 -0.15945837 -0.37678852
22 -0.56756829 -0.15945837
23 -0.01291849 -0.56756829
24 -0.30715782 -0.01291849
25 0.63242892 -0.30715782
26 0.86598106 0.63242892
27 0.02834222 0.86598106
28 -0.30275317 0.02834222
29 -0.16730270 -0.30275317
30 -0.47566811 -0.16730270
31 -0.79547875 -0.47566811
32 -0.61351475 -0.79547875
33 0.46677351 -0.61351475
34 -0.20647549 0.46677351
35 -0.34026274 -0.20647549
36 -0.75355155 -0.34026274
37 0.06331221 -0.75355155
38 0.10889593 0.06331221
39 1.00943216 0.10889593
40 2.35152326 1.00943216
41 1.94817177 2.35152326
42 0.72873866 1.94817177
43 0.63433250 0.72873866
44 0.45166297 0.63433250
45 -0.05349762 0.45166297
46 0.78033466 -0.05349762
47 -0.01705724 0.78033466
48 -0.29839425 -0.01705724
49 0.03265786 -0.29839425
50 0.35070930 0.03265786
51 -0.27704328 0.35070930
52 0.40166340 -0.27704328
53 -0.18066476 0.40166340
54 -1.22394348 -0.18066476
55 -1.02160831 -1.22394348
56 -0.58579267 -1.02160831
57 -0.56394721 -0.58579267
58 NA -0.56394721
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.21352657 -0.24268364
[2,] 0.11484774 -0.21352657
[3,] -0.24260447 0.11484774
[4,] 0.25123533 -0.24260447
[5,] 0.08627380 0.25123533
[6,] 1.25355161 0.08627380
[7,] 0.64304338 1.25355161
[8,] 1.03027223 0.64304338
[9,] 0.67393657 1.03027223
[10,] 0.25880747 0.67393657
[11,] 0.30823971 0.25880747
[12,] -0.34073216 0.30823971
[13,] -0.86843481 -0.34073216
[14,] -0.87346346 -0.86843481
[15,] -0.97070692 -0.87346346
[16,] -0.96729444 -0.97070692
[17,] -0.72910721 -0.96729444
[18,] -0.29702098 -0.72910721
[19,] -0.42874601 -0.29702098
[20,] -0.37678852 -0.42874601
[21,] -0.15945837 -0.37678852
[22,] -0.56756829 -0.15945837
[23,] -0.01291849 -0.56756829
[24,] -0.30715782 -0.01291849
[25,] 0.63242892 -0.30715782
[26,] 0.86598106 0.63242892
[27,] 0.02834222 0.86598106
[28,] -0.30275317 0.02834222
[29,] -0.16730270 -0.30275317
[30,] -0.47566811 -0.16730270
[31,] -0.79547875 -0.47566811
[32,] -0.61351475 -0.79547875
[33,] 0.46677351 -0.61351475
[34,] -0.20647549 0.46677351
[35,] -0.34026274 -0.20647549
[36,] -0.75355155 -0.34026274
[37,] 0.06331221 -0.75355155
[38,] 0.10889593 0.06331221
[39,] 1.00943216 0.10889593
[40,] 2.35152326 1.00943216
[41,] 1.94817177 2.35152326
[42,] 0.72873866 1.94817177
[43,] 0.63433250 0.72873866
[44,] 0.45166297 0.63433250
[45,] -0.05349762 0.45166297
[46,] 0.78033466 -0.05349762
[47,] -0.01705724 0.78033466
[48,] -0.29839425 -0.01705724
[49,] 0.03265786 -0.29839425
[50,] 0.35070930 0.03265786
[51,] -0.27704328 0.35070930
[52,] 0.40166340 -0.27704328
[53,] -0.18066476 0.40166340
[54,] -1.22394348 -0.18066476
[55,] -1.02160831 -1.22394348
[56,] -0.58579267 -1.02160831
[57,] -0.56394721 -0.58579267
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.21352657 -0.24268364
2 0.11484774 -0.21352657
3 -0.24260447 0.11484774
4 0.25123533 -0.24260447
5 0.08627380 0.25123533
6 1.25355161 0.08627380
7 0.64304338 1.25355161
8 1.03027223 0.64304338
9 0.67393657 1.03027223
10 0.25880747 0.67393657
11 0.30823971 0.25880747
12 -0.34073216 0.30823971
13 -0.86843481 -0.34073216
14 -0.87346346 -0.86843481
15 -0.97070692 -0.87346346
16 -0.96729444 -0.97070692
17 -0.72910721 -0.96729444
18 -0.29702098 -0.72910721
19 -0.42874601 -0.29702098
20 -0.37678852 -0.42874601
21 -0.15945837 -0.37678852
22 -0.56756829 -0.15945837
23 -0.01291849 -0.56756829
24 -0.30715782 -0.01291849
25 0.63242892 -0.30715782
26 0.86598106 0.63242892
27 0.02834222 0.86598106
28 -0.30275317 0.02834222
29 -0.16730270 -0.30275317
30 -0.47566811 -0.16730270
31 -0.79547875 -0.47566811
32 -0.61351475 -0.79547875
33 0.46677351 -0.61351475
34 -0.20647549 0.46677351
35 -0.34026274 -0.20647549
36 -0.75355155 -0.34026274
37 0.06331221 -0.75355155
38 0.10889593 0.06331221
39 1.00943216 0.10889593
40 2.35152326 1.00943216
41 1.94817177 2.35152326
42 0.72873866 1.94817177
43 0.63433250 0.72873866
44 0.45166297 0.63433250
45 -0.05349762 0.45166297
46 0.78033466 -0.05349762
47 -0.01705724 0.78033466
48 -0.29839425 -0.01705724
49 0.03265786 -0.29839425
50 0.35070930 0.03265786
51 -0.27704328 0.35070930
52 0.40166340 -0.27704328
53 -0.18066476 0.40166340
54 -1.22394348 -0.18066476
55 -1.02160831 -1.22394348
56 -0.58579267 -1.02160831
57 -0.56394721 -0.58579267
> 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/715nn1290545353.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/815nn1290545353.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/915nn1290545353.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/10temq1290545353.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/11ff3v1290545353.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/12ix111290545353.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/13e7za1290545353.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/1408xy1290545353.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/15l8wm1290545353.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/16z0uu1290545353.tab")
+ }
>
> try(system("convert tmp/1mvpe1290545353.ps tmp/1mvpe1290545353.png",intern=TRUE))
character(0)
> try(system("convert tmp/2mvpe1290545353.ps tmp/2mvpe1290545353.png",intern=TRUE))
character(0)
> try(system("convert tmp/3fmoz1290545353.ps tmp/3fmoz1290545353.png",intern=TRUE))
character(0)
> try(system("convert tmp/4fmoz1290545353.ps tmp/4fmoz1290545353.png",intern=TRUE))
character(0)
> try(system("convert tmp/5fmoz1290545353.ps tmp/5fmoz1290545353.png",intern=TRUE))
character(0)
> try(system("convert tmp/68wo21290545353.ps tmp/68wo21290545353.png",intern=TRUE))
character(0)
> try(system("convert tmp/715nn1290545353.ps tmp/715nn1290545353.png",intern=TRUE))
character(0)
> try(system("convert tmp/815nn1290545353.ps tmp/815nn1290545353.png",intern=TRUE))
character(0)
> try(system("convert tmp/915nn1290545353.ps tmp/915nn1290545353.png",intern=TRUE))
character(0)
> try(system("convert tmp/10temq1290545353.ps tmp/10temq1290545353.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.408 1.565 6.025