R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(107.11
+ ,236.67
+ ,8.92
+ ,1
+ ,122.23
+ ,258.1
+ ,9.32
+ ,2
+ ,134.69
+ ,241.52
+ ,8.9
+ ,3
+ ,128.79
+ ,190.71
+ ,8.53
+ ,4
+ ,126.16
+ ,200.32
+ ,8.51
+ ,5
+ ,119.98
+ ,223.41
+ ,9.03
+ ,6
+ ,108.45
+ ,201.38
+ ,9.6
+ ,7
+ ,108.43
+ ,211.83
+ ,9.88
+ ,8
+ ,98.17
+ ,224.41
+ ,10.81
+ ,9
+ ,106.09
+ ,211.57
+ ,11.61
+ ,10
+ ,108.81
+ ,194.77
+ ,11.81
+ ,11
+ ,103.03
+ ,201.86
+ ,13.93
+ ,12
+ ,124.36
+ ,225
+ ,16.19
+ ,13
+ ,118.52
+ ,278.9
+ ,18.05
+ ,14
+ ,112.2
+ ,259.74
+ ,17.08
+ ,15
+ ,114.71
+ ,230.45
+ ,17.46
+ ,16
+ ,107.96
+ ,238.26
+ ,16.9
+ ,17
+ ,101.21
+ ,250.14
+ ,15.69
+ ,18
+ ,102.77
+ ,263.81
+ ,15.86
+ ,19
+ ,112.13
+ ,247.22
+ ,12.98
+ ,20
+ ,109.36
+ ,229.81
+ ,12.31
+ ,21
+ ,110.91
+ ,224.27
+ ,11.51
+ ,22
+ ,123.57
+ ,213.23
+ ,11.73
+ ,23
+ ,129.95
+ ,239.57
+ ,11.7
+ ,24
+ ,124.46
+ ,249.7
+ ,10.9
+ ,25
+ ,122.34
+ ,212.5
+ ,10.57
+ ,26
+ ,116.61
+ ,203.27
+ ,10.37
+ ,27
+ ,114.59
+ ,192.05
+ ,9.59
+ ,28
+ ,112.52
+ ,190.04
+ ,9.09
+ ,29
+ ,118.67
+ ,202.05
+ ,9.26
+ ,30
+ ,116.8
+ ,211.91
+ ,9.9
+ ,31
+ ,123.63
+ ,210.39
+ ,9.61
+ ,32
+ ,128.04
+ ,231.25
+ ,9.85
+ ,33
+ ,134.57
+ ,224.3
+ ,9.99
+ ,34
+ ,130.33
+ ,209.64
+ ,9.9
+ ,35
+ ,136.47
+ ,206.05
+ ,10.45
+ ,36
+ ,139.05
+ ,229.7
+ ,11.66
+ ,37
+ ,158.21
+ ,264.67
+ ,13.61
+ ,38
+ ,148.07
+ ,246.29
+ ,12.88
+ ,39
+ ,137.74
+ ,260.91
+ ,12.52
+ ,40
+ ,139.74
+ ,265.14
+ ,10.93
+ ,41
+ ,144.08
+ ,284.52
+ ,12.07
+ ,42
+ ,145.35
+ ,287.48
+ ,13.21
+ ,43
+ ,145.77
+ ,321.9
+ ,13.68
+ ,44
+ ,140.56
+ ,321.59
+ ,14.02
+ ,45
+ ,121.41
+ ,282.39
+ ,11.7
+ ,46
+ ,120.44
+ ,241
+ ,11.83
+ ,47
+ ,116.97
+ ,228.48
+ ,11.32
+ ,48
+ ,128.03
+ ,261.59
+ ,12.24
+ ,49
+ ,128.51
+ ,270
+ ,13.31
+ ,50
+ ,127.76
+ ,262.86
+ ,12.93
+ ,51
+ ,134.58
+ ,277.41
+ ,13.47
+ ,52
+ ,147.64
+ ,288
+ ,15.47
+ ,53
+ ,144.46
+ ,287.14
+ ,16.58
+ ,54
+ ,137.6
+ ,337.65
+ ,17.8
+ ,55
+ ,146.87
+ ,328.38
+ ,21.72
+ ,56
+ ,145.67
+ ,374.41
+ ,23.45
+ ,57
+ ,151.95
+ ,344.77
+ ,23.16
+ ,58
+ ,150.23
+ ,361.05
+ ,22.77
+ ,59
+ ,155.86
+ ,374.22
+ ,24.9
+ ,60)
+ ,dim=c(4
+ ,60)
+ ,dimnames=list(c('Coffee'
+ ,'Tea'
+ ,'Sugar'
+ ,'Month')
+ ,1:60))
> y <- array(NA,dim=c(4,60),dimnames=list(c('Coffee','Tea','Sugar','Month'),1:60))
> 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
> 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
Coffee Tea Sugar Month
1 107.11 236.67 8.92 1
2 122.23 258.10 9.32 2
3 134.69 241.52 8.90 3
4 128.79 190.71 8.53 4
5 126.16 200.32 8.51 5
6 119.98 223.41 9.03 6
7 108.45 201.38 9.60 7
8 108.43 211.83 9.88 8
9 98.17 224.41 10.81 9
10 106.09 211.57 11.61 10
11 108.81 194.77 11.81 11
12 103.03 201.86 13.93 12
13 124.36 225.00 16.19 13
14 118.52 278.90 18.05 14
15 112.20 259.74 17.08 15
16 114.71 230.45 17.46 16
17 107.96 238.26 16.90 17
18 101.21 250.14 15.69 18
19 102.77 263.81 15.86 19
20 112.13 247.22 12.98 20
21 109.36 229.81 12.31 21
22 110.91 224.27 11.51 22
23 123.57 213.23 11.73 23
24 129.95 239.57 11.70 24
25 124.46 249.70 10.90 25
26 122.34 212.50 10.57 26
27 116.61 203.27 10.37 27
28 114.59 192.05 9.59 28
29 112.52 190.04 9.09 29
30 118.67 202.05 9.26 30
31 116.80 211.91 9.90 31
32 123.63 210.39 9.61 32
33 128.04 231.25 9.85 33
34 134.57 224.30 9.99 34
35 130.33 209.64 9.90 35
36 136.47 206.05 10.45 36
37 139.05 229.70 11.66 37
38 158.21 264.67 13.61 38
39 148.07 246.29 12.88 39
40 137.74 260.91 12.52 40
41 139.74 265.14 10.93 41
42 144.08 284.52 12.07 42
43 145.35 287.48 13.21 43
44 145.77 321.90 13.68 44
45 140.56 321.59 14.02 45
46 121.41 282.39 11.70 46
47 120.44 241.00 11.83 47
48 116.97 228.48 11.32 48
49 128.03 261.59 12.24 49
50 128.51 270.00 13.31 50
51 127.76 262.86 12.93 51
52 134.58 277.41 13.47 52
53 147.64 288.00 15.47 53
54 144.46 287.14 16.58 54
55 137.60 337.65 17.80 55
56 146.87 328.38 21.72 56
57 145.67 374.41 23.45 57
58 151.95 344.77 23.16 58
59 150.23 361.05 22.77 59
60 155.86 374.22 24.90 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Tea Sugar Month
81.5752 0.1855 -1.1106 0.4169
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-18.7197 -6.1637 -0.4629 5.9862 26.8231
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 81.57518 8.38126 9.733 1.21e-13 ***
Tea 0.18546 0.05627 3.296 0.001706 **
Sugar -1.11064 0.56067 -1.981 0.052518 .
Month 0.41689 0.10377 4.018 0.000177 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.07 on 56 degrees of freedom
Multiple R-squared: 0.5912, Adjusted R-squared: 0.5693
F-statistic: 26.99 on 3 and 56 DF, p-value: 6.25e-11
> 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.6537039 0.692592174 0.346296087
[2,] 0.4883469 0.976693792 0.511653104
[3,] 0.3757211 0.751442141 0.624278929
[4,] 0.6038689 0.792262247 0.396131123
[5,] 0.5923909 0.815218108 0.407609054
[6,] 0.5452392 0.909521532 0.454760766
[7,] 0.7833480 0.433304054 0.216652027
[8,] 0.7252902 0.549419632 0.274709816
[9,] 0.6396362 0.720727653 0.360363826
[10,] 0.5593758 0.881248372 0.440624186
[11,] 0.4696365 0.939272975 0.530363513
[12,] 0.4338394 0.867678716 0.566160642
[13,] 0.4611160 0.922232043 0.538883978
[14,] 0.5675127 0.864974663 0.432487332
[15,] 0.5987035 0.802592902 0.401296451
[16,] 0.6383617 0.723276571 0.361638286
[17,] 0.7076921 0.584615810 0.292307905
[18,] 0.7549358 0.490128307 0.245064154
[19,] 0.7335760 0.532848018 0.266424009
[20,] 0.6886539 0.622692263 0.311346132
[21,] 0.6700081 0.659983746 0.329991873
[22,] 0.6675176 0.664964817 0.332482409
[23,] 0.7149729 0.570054262 0.285027131
[24,] 0.7284620 0.543076031 0.271538016
[25,] 0.8448504 0.310299201 0.155149601
[26,] 0.8826042 0.234791559 0.117395779
[27,] 0.9132921 0.173415757 0.086707879
[28,] 0.9172192 0.165561633 0.082780816
[29,] 0.9206362 0.158727574 0.079363787
[30,] 0.9205880 0.158823915 0.079411957
[31,] 0.9243443 0.151311320 0.075655660
[32,] 0.9730101 0.053979838 0.026989919
[33,] 0.9735848 0.052830385 0.026415193
[34,] 0.9573761 0.085247799 0.042623900
[35,] 0.9391124 0.121775208 0.060887604
[36,] 0.9276442 0.144711523 0.072355761
[37,] 0.9380798 0.123840465 0.061920232
[38,] 0.9560490 0.087901923 0.043950961
[39,] 0.9933388 0.013322452 0.006661226
[40,] 0.9946777 0.010644596 0.005322298
[41,] 0.9899573 0.020085386 0.010042693
[42,] 0.9960173 0.007965445 0.003982723
[43,] 0.9899173 0.020165360 0.010082680
[44,] 0.9756963 0.048607399 0.024303699
[45,] 0.9860458 0.027908448 0.013954224
[46,] 0.9756119 0.048776218 0.024388109
[47,] 0.9950498 0.009900348 0.004950174
> postscript(file="/var/www/rcomp/tmp/15zq71292276289.ps",horizontal=F,onefile=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/25zq71292276289.ps",horizontal=F,onefile=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/3g87r1292276289.ps",horizontal=F,onefile=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/4g87r1292276289.ps",horizontal=F,onefile=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/5g87r1292276289.ps",horizontal=F,onefile=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 = 60
Frequency = 1
1 2 3 4 5 6
-8.8679769 2.3049798 16.9565472 19.6519453 14.8005703 4.4989404
7 8 9 10 11 12
-2.7291987 -4.7931674 -16.7702485 -5.9973179 -0.3563508 -5.5135915
13 14 15 16 17 18
13.6180234 -0.5693706 -4.8301695 3.1171096 -6.1201842 -16.8342181
19 20 21 22 23 24
-18.0375390 -9.2162974 -9.9184587 -8.6464146 5.8885152 6.9332867
25 26 27 28 29 30
-1.7408287 2.2548834 -2.4023393 -3.6246692 -6.2941065 -2.5995636
31 32 33 34 35 36
-6.0042796 0.3686426 0.7596088 8.3171555 6.2791516 13.2789159
37 38 39 40 41 42
12.3997716 26.8230947 18.8641912 5.0060427 4.0387343 5.6337596
43 44 45 46 47 48
7.2040394 1.3456147 -3.8461651 -18.7197113 -12.2860260 -14.4173842
49 50 51 52 53 54
-8.8930671 -9.2012897 -9.4660397 -5.1616277 7.7387442 5.5341622
55 56 57 58 59 60
-9.7553333 5.1707093 -3.0614974 7.9765622 2.3872308 7.5234992
> postscript(file="/var/www/rcomp/tmp/68hoc1292276289.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -8.8679769 NA
1 2.3049798 -8.8679769
2 16.9565472 2.3049798
3 19.6519453 16.9565472
4 14.8005703 19.6519453
5 4.4989404 14.8005703
6 -2.7291987 4.4989404
7 -4.7931674 -2.7291987
8 -16.7702485 -4.7931674
9 -5.9973179 -16.7702485
10 -0.3563508 -5.9973179
11 -5.5135915 -0.3563508
12 13.6180234 -5.5135915
13 -0.5693706 13.6180234
14 -4.8301695 -0.5693706
15 3.1171096 -4.8301695
16 -6.1201842 3.1171096
17 -16.8342181 -6.1201842
18 -18.0375390 -16.8342181
19 -9.2162974 -18.0375390
20 -9.9184587 -9.2162974
21 -8.6464146 -9.9184587
22 5.8885152 -8.6464146
23 6.9332867 5.8885152
24 -1.7408287 6.9332867
25 2.2548834 -1.7408287
26 -2.4023393 2.2548834
27 -3.6246692 -2.4023393
28 -6.2941065 -3.6246692
29 -2.5995636 -6.2941065
30 -6.0042796 -2.5995636
31 0.3686426 -6.0042796
32 0.7596088 0.3686426
33 8.3171555 0.7596088
34 6.2791516 8.3171555
35 13.2789159 6.2791516
36 12.3997716 13.2789159
37 26.8230947 12.3997716
38 18.8641912 26.8230947
39 5.0060427 18.8641912
40 4.0387343 5.0060427
41 5.6337596 4.0387343
42 7.2040394 5.6337596
43 1.3456147 7.2040394
44 -3.8461651 1.3456147
45 -18.7197113 -3.8461651
46 -12.2860260 -18.7197113
47 -14.4173842 -12.2860260
48 -8.8930671 -14.4173842
49 -9.2012897 -8.8930671
50 -9.4660397 -9.2012897
51 -5.1616277 -9.4660397
52 7.7387442 -5.1616277
53 5.5341622 7.7387442
54 -9.7553333 5.5341622
55 5.1707093 -9.7553333
56 -3.0614974 5.1707093
57 7.9765622 -3.0614974
58 2.3872308 7.9765622
59 7.5234992 2.3872308
60 NA 7.5234992
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2.3049798 -8.8679769
[2,] 16.9565472 2.3049798
[3,] 19.6519453 16.9565472
[4,] 14.8005703 19.6519453
[5,] 4.4989404 14.8005703
[6,] -2.7291987 4.4989404
[7,] -4.7931674 -2.7291987
[8,] -16.7702485 -4.7931674
[9,] -5.9973179 -16.7702485
[10,] -0.3563508 -5.9973179
[11,] -5.5135915 -0.3563508
[12,] 13.6180234 -5.5135915
[13,] -0.5693706 13.6180234
[14,] -4.8301695 -0.5693706
[15,] 3.1171096 -4.8301695
[16,] -6.1201842 3.1171096
[17,] -16.8342181 -6.1201842
[18,] -18.0375390 -16.8342181
[19,] -9.2162974 -18.0375390
[20,] -9.9184587 -9.2162974
[21,] -8.6464146 -9.9184587
[22,] 5.8885152 -8.6464146
[23,] 6.9332867 5.8885152
[24,] -1.7408287 6.9332867
[25,] 2.2548834 -1.7408287
[26,] -2.4023393 2.2548834
[27,] -3.6246692 -2.4023393
[28,] -6.2941065 -3.6246692
[29,] -2.5995636 -6.2941065
[30,] -6.0042796 -2.5995636
[31,] 0.3686426 -6.0042796
[32,] 0.7596088 0.3686426
[33,] 8.3171555 0.7596088
[34,] 6.2791516 8.3171555
[35,] 13.2789159 6.2791516
[36,] 12.3997716 13.2789159
[37,] 26.8230947 12.3997716
[38,] 18.8641912 26.8230947
[39,] 5.0060427 18.8641912
[40,] 4.0387343 5.0060427
[41,] 5.6337596 4.0387343
[42,] 7.2040394 5.6337596
[43,] 1.3456147 7.2040394
[44,] -3.8461651 1.3456147
[45,] -18.7197113 -3.8461651
[46,] -12.2860260 -18.7197113
[47,] -14.4173842 -12.2860260
[48,] -8.8930671 -14.4173842
[49,] -9.2012897 -8.8930671
[50,] -9.4660397 -9.2012897
[51,] -5.1616277 -9.4660397
[52,] 7.7387442 -5.1616277
[53,] 5.5341622 7.7387442
[54,] -9.7553333 5.5341622
[55,] 5.1707093 -9.7553333
[56,] -3.0614974 5.1707093
[57,] 7.9765622 -3.0614974
[58,] 2.3872308 7.9765622
[59,] 7.5234992 2.3872308
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2.3049798 -8.8679769
2 16.9565472 2.3049798
3 19.6519453 16.9565472
4 14.8005703 19.6519453
5 4.4989404 14.8005703
6 -2.7291987 4.4989404
7 -4.7931674 -2.7291987
8 -16.7702485 -4.7931674
9 -5.9973179 -16.7702485
10 -0.3563508 -5.9973179
11 -5.5135915 -0.3563508
12 13.6180234 -5.5135915
13 -0.5693706 13.6180234
14 -4.8301695 -0.5693706
15 3.1171096 -4.8301695
16 -6.1201842 3.1171096
17 -16.8342181 -6.1201842
18 -18.0375390 -16.8342181
19 -9.2162974 -18.0375390
20 -9.9184587 -9.2162974
21 -8.6464146 -9.9184587
22 5.8885152 -8.6464146
23 6.9332867 5.8885152
24 -1.7408287 6.9332867
25 2.2548834 -1.7408287
26 -2.4023393 2.2548834
27 -3.6246692 -2.4023393
28 -6.2941065 -3.6246692
29 -2.5995636 -6.2941065
30 -6.0042796 -2.5995636
31 0.3686426 -6.0042796
32 0.7596088 0.3686426
33 8.3171555 0.7596088
34 6.2791516 8.3171555
35 13.2789159 6.2791516
36 12.3997716 13.2789159
37 26.8230947 12.3997716
38 18.8641912 26.8230947
39 5.0060427 18.8641912
40 4.0387343 5.0060427
41 5.6337596 4.0387343
42 7.2040394 5.6337596
43 1.3456147 7.2040394
44 -3.8461651 1.3456147
45 -18.7197113 -3.8461651
46 -12.2860260 -18.7197113
47 -14.4173842 -12.2860260
48 -8.8930671 -14.4173842
49 -9.2012897 -8.8930671
50 -9.4660397 -9.2012897
51 -5.1616277 -9.4660397
52 7.7387442 -5.1616277
53 5.5341622 7.7387442
54 -9.7553333 5.5341622
55 5.1707093 -9.7553333
56 -3.0614974 5.1707093
57 7.9765622 -3.0614974
58 2.3872308 7.9765622
59 7.5234992 2.3872308
> 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/7196f1292276289.ps",horizontal=F,onefile=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/8196f1292276289.ps",horizontal=F,onefile=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/9196f1292276289.ps",horizontal=F,onefile=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/10u0501292276289.ps",horizontal=F,onefile=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/11filo1292276289.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/12i12u1292276289.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/13etz31292276289.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/14pkhn1292276289.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/15luxw1292276289.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/167cvk1292276289.tab")
+ }
>
> try(system("convert tmp/15zq71292276289.ps tmp/15zq71292276289.png",intern=TRUE))
character(0)
> try(system("convert tmp/25zq71292276289.ps tmp/25zq71292276289.png",intern=TRUE))
character(0)
> try(system("convert tmp/3g87r1292276289.ps tmp/3g87r1292276289.png",intern=TRUE))
character(0)
> try(system("convert tmp/4g87r1292276289.ps tmp/4g87r1292276289.png",intern=TRUE))
character(0)
> try(system("convert tmp/5g87r1292276289.ps tmp/5g87r1292276289.png",intern=TRUE))
character(0)
> try(system("convert tmp/68hoc1292276289.ps tmp/68hoc1292276289.png",intern=TRUE))
character(0)
> try(system("convert tmp/7196f1292276289.ps tmp/7196f1292276289.png",intern=TRUE))
character(0)
> try(system("convert tmp/8196f1292276289.ps tmp/8196f1292276289.png",intern=TRUE))
character(0)
> try(system("convert tmp/9196f1292276289.ps tmp/9196f1292276289.png",intern=TRUE))
character(0)
> try(system("convert tmp/10u0501292276289.ps tmp/10u0501292276289.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.040 1.780 4.792