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(9
+ ,13
+ ,2528
+ ,80
+ ,15.3
+ ,9
+ ,12
+ ,3333
+ ,83
+ ,19.4
+ ,9
+ ,54
+ ,19611
+ ,96
+ ,19.5
+ ,9
+ ,19
+ ,3570
+ ,56
+ ,17
+ ,9
+ ,37
+ ,1722
+ ,43
+ ,19.3
+ ,9
+ ,2
+ ,583
+ ,51
+ ,12.9
+ ,9
+ ,72
+ ,4790
+ ,91
+ ,16.7
+ ,9
+ ,164
+ ,35971
+ ,81
+ ,13.8
+ ,9
+ ,18
+ ,25440
+ ,120
+ ,13.7
+ ,9
+ ,1
+ ,2217
+ ,46
+ ,14.3
+ ,9
+ ,53
+ ,1971
+ ,56
+ ,22.2
+ ,9
+ ,16
+ ,12620
+ ,37
+ ,16.8
+ ,9
+ ,32
+ ,19046
+ ,120
+ ,18
+ ,10
+ ,21
+ ,8612
+ ,103
+ ,15
+ ,10
+ ,23
+ ,3896
+ ,105
+ ,18.4
+ ,10
+ ,18
+ ,6298
+ ,42
+ ,18.2
+ ,10
+ ,112
+ ,27350
+ ,65
+ ,24.1
+ ,10
+ ,25
+ ,4145
+ ,51
+ ,23
+ ,10
+ ,5
+ ,1175
+ ,57
+ ,21.8
+ ,10
+ ,26
+ ,8297
+ ,60
+ ,19.1
+ ,10
+ ,8
+ ,7814
+ ,160
+ ,22.6
+ ,10
+ ,15
+ ,1745
+ ,48
+ ,14.3
+ ,10
+ ,11
+ ,5046
+ ,109
+ ,19
+ ,10
+ ,11
+ ,18943
+ ,50
+ ,18.5
+ ,10
+ ,87
+ ,8624
+ ,78
+ ,21.3
+ ,10
+ ,33
+ ,2225
+ ,41
+ ,20.5
+ ,11
+ ,22
+ ,12659
+ ,65
+ ,18
+ ,11
+ ,98
+ ,1967
+ ,50
+ ,16.8
+ ,11
+ ,1
+ ,1172
+ ,73
+ ,20.5
+ ,11
+ ,5
+ ,639
+ ,26
+ ,20.1
+ ,11
+ ,1
+ ,7056
+ ,60
+ ,24.5
+ ,11
+ ,38
+ ,1934
+ ,85
+ ,12
+ ,11
+ ,30
+ ,6260
+ ,133
+ ,21
+ ,11
+ ,12
+ ,424
+ ,62
+ ,20.2
+ ,11
+ ,24
+ ,3488
+ ,44
+ ,24
+ ,11
+ ,6
+ ,3330
+ ,67
+ ,14.9
+ ,11
+ ,15
+ ,2227
+ ,54
+ ,24
+ ,11
+ ,38
+ ,8115
+ ,110
+ ,20.5
+ ,11
+ ,84
+ ,1600
+ ,56
+ ,19.5
+ ,12
+ ,3
+ ,15305
+ ,85
+ ,17.5
+ ,12
+ ,18
+ ,7121
+ ,58
+ ,16
+ ,12
+ ,63
+ ,5794
+ ,34
+ ,17.5
+ ,12
+ ,239
+ ,8636
+ ,150
+ ,18.1
+ ,12
+ ,234
+ ,4803
+ ,93
+ ,24.3
+ ,12
+ ,6
+ ,1097
+ ,53
+ ,13.1
+ ,12
+ ,76
+ ,9765
+ ,130
+ ,16.9
+ ,12
+ ,25
+ ,4266
+ ,68
+ ,17
+ ,12
+ ,8
+ ,1507
+ ,51
+ ,21
+ ,12
+ ,23
+ ,3836
+ ,121
+ ,18.5
+ ,12
+ ,16
+ ,17419
+ ,48
+ ,18
+ ,12
+ ,6
+ ,8735
+ ,63
+ ,20.8
+ ,12
+ ,100
+ ,22550
+ ,107
+ ,23
+ ,13
+ ,80
+ ,9961
+ ,79
+ ,21.8
+ ,13
+ ,28
+ ,4706
+ ,61
+ ,19.7
+ ,13
+ ,48
+ ,4011
+ ,52
+ ,18.9
+ ,13
+ ,18
+ ,6949
+ ,100
+ ,18.6
+ ,13
+ ,36
+ ,11405
+ ,70
+ ,23.6
+ ,13
+ ,19
+ ,904
+ ,39
+ ,19.2
+ ,13
+ ,32
+ ,3332
+ ,73
+ ,17.7
+ ,13
+ ,3
+ ,575
+ ,33
+ ,18
+ ,13
+ ,106
+ ,29708
+ ,73
+ ,21.4
+ ,13
+ ,62
+ ,2511
+ ,60
+ ,17.7
+ ,13
+ ,23
+ ,18422
+ ,45
+ ,20.1
+ ,13
+ ,2
+ ,6311
+ ,46
+ ,18.5
+ ,13
+ ,26
+ ,1450
+ ,60
+ ,18.6
+ ,14
+ ,20
+ ,4106
+ ,96
+ ,15.4
+ ,14
+ ,38
+ ,10274
+ ,90
+ ,15
+ ,14
+ ,19
+ ,510
+ ,82
+ ,26.5)
+ ,dim=c(5
+ ,68)
+ ,dimnames=list(c('Month'
+ ,'Fish'
+ ,'Acre'
+ ,'Depth'
+ ,'Temp')
+ ,1:68))
> y <- array(NA,dim=c(5,68),dimnames=list(c('Month','Fish','Acre','Depth','Temp'),1:68))
> 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 = '2'
> 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
Fish Month Acre Depth Temp
1 13 9 2528 80 15.3
2 12 9 3333 83 19.4
3 54 9 19611 96 19.5
4 19 9 3570 56 17.0
5 37 9 1722 43 19.3
6 2 9 583 51 12.9
7 72 9 4790 91 16.7
8 164 9 35971 81 13.8
9 18 9 25440 120 13.7
10 1 9 2217 46 14.3
11 53 9 1971 56 22.2
12 16 9 12620 37 16.8
13 32 9 19046 120 18.0
14 21 10 8612 103 15.0
15 23 10 3896 105 18.4
16 18 10 6298 42 18.2
17 112 10 27350 65 24.1
18 25 10 4145 51 23.0
19 5 10 1175 57 21.8
20 26 10 8297 60 19.1
21 8 10 7814 160 22.6
22 15 10 1745 48 14.3
23 11 10 5046 109 19.0
24 11 10 18943 50 18.5
25 87 10 8624 78 21.3
26 33 10 2225 41 20.5
27 22 11 12659 65 18.0
28 98 11 1967 50 16.8
29 1 11 1172 73 20.5
30 5 11 639 26 20.1
31 1 11 7056 60 24.5
32 38 11 1934 85 12.0
33 30 11 6260 133 21.0
34 12 11 424 62 20.2
35 24 11 3488 44 24.0
36 6 11 3330 67 14.9
37 15 11 2227 54 24.0
38 38 11 8115 110 20.5
39 84 11 1600 56 19.5
40 3 12 15305 85 17.5
41 18 12 7121 58 16.0
42 63 12 5794 34 17.5
43 239 12 8636 150 18.1
44 234 12 4803 93 24.3
45 6 12 1097 53 13.1
46 76 12 9765 130 16.9
47 25 12 4266 68 17.0
48 8 12 1507 51 21.0
49 23 12 3836 121 18.5
50 16 12 17419 48 18.0
51 6 12 8735 63 20.8
52 100 12 22550 107 23.0
53 80 13 9961 79 21.8
54 28 13 4706 61 19.7
55 48 13 4011 52 18.9
56 18 13 6949 100 18.6
57 36 13 11405 70 23.6
58 19 13 904 39 19.2
59 32 13 3332 73 17.7
60 3 13 575 33 18.0
61 106 13 29708 73 21.4
62 62 13 2511 60 17.7
63 23 13 18422 45 20.1
64 2 13 6311 46 18.5
65 26 13 1450 60 18.6
66 20 14 4106 96 15.4
67 38 14 10274 90 15.0
68 19 14 510 82 26.5
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Month Acre Depth Temp
-63.642389 2.312471 0.001793 0.374300 1.907180
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-68.482 -21.567 -7.543 13.893 180.127
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.364e+01 4.944e+01 -1.287 0.2027
Month 2.312e+00 3.582e+00 0.646 0.5209
Acre 1.793e-03 7.178e-04 2.498 0.0151 *
Depth 3.743e-01 1.889e-01 1.982 0.0519 .
Temp 1.907e+00 1.728e+00 1.104 0.2738
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 43.76 on 63 degrees of freedom
Multiple R-squared: 0.1947, Adjusted R-squared: 0.1436
F-statistic: 3.808 on 4 and 63 DF, p-value: 0.00778
> 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.402228188 0.8044563758 0.5977718121
[2,] 0.551718555 0.8965628894 0.4482814447
[3,] 0.448029000 0.8960580005 0.5519709998
[4,] 0.325185345 0.6503706894 0.6748146553
[5,] 0.371001074 0.7420021477 0.6289989262
[6,] 0.297462160 0.5949243193 0.7025378404
[7,] 0.208701399 0.4174027981 0.7912986010
[8,] 0.141063522 0.2821270443 0.8589364778
[9,] 0.104708029 0.2094160581 0.8952919710
[10,] 0.070487273 0.1409745456 0.9295127272
[11,] 0.044845459 0.0896909176 0.9551545412
[12,] 0.028679470 0.0573589390 0.9713205305
[13,] 0.016424029 0.0328480575 0.9835759712
[14,] 0.014095470 0.0281909397 0.9859045301
[15,] 0.008104295 0.0162085895 0.9918957053
[16,] 0.005339409 0.0106788171 0.9946605914
[17,] 0.009150744 0.0183014883 0.9908492558
[18,] 0.015151854 0.0303037079 0.9848481461
[19,] 0.009189191 0.0183783825 0.9908108087
[20,] 0.005392695 0.0107853906 0.9946073047
[21,] 0.037779032 0.0755580644 0.9622209678
[22,] 0.029347187 0.0586943732 0.9706528134
[23,] 0.020649964 0.0412999279 0.9793500361
[24,] 0.020457644 0.0409152883 0.9795423558
[25,] 0.015295363 0.0305907253 0.9847046374
[26,] 0.016348743 0.0326974869 0.9836512566
[27,] 0.010852022 0.0217040434 0.9891479783
[28,] 0.006771367 0.0135427332 0.9932286334
[29,] 0.004460156 0.0089203119 0.9955398441
[30,] 0.003444699 0.0068893978 0.9965553011
[31,] 0.004488106 0.0089762116 0.9955118942
[32,] 0.007186771 0.0143735414 0.9928132293
[33,] 0.011700800 0.0234015997 0.9882992002
[34,] 0.007345073 0.0146901451 0.9926549274
[35,] 0.007867066 0.0157341330 0.9921329335
[36,] 0.511693106 0.9766137875 0.4883068938
[37,] 0.999827549 0.0003449016 0.0001724508
[38,] 0.999621801 0.0007563980 0.0003781990
[39,] 0.999391862 0.0012162770 0.0006081385
[40,] 0.998746497 0.0025070058 0.0012535029
[41,] 0.997490287 0.0050194262 0.0025097131
[42,] 0.995777740 0.0084445204 0.0042222602
[43,] 0.995486062 0.0090278763 0.0045139381
[44,] 0.997358006 0.0052839872 0.0026419936
[45,] 0.993989136 0.0120217287 0.0060108644
[46,] 0.993797565 0.0124048693 0.0062024346
[47,] 0.986118393 0.0277632133 0.0138816067
[48,] 0.980091637 0.0398167256 0.0199083628
[49,] 0.986422539 0.0271549212 0.0135774606
[50,] 0.980154128 0.0396917436 0.0198458718
[51,] 0.959501101 0.0809977975 0.0404988988
[52,] 0.942990398 0.1140192046 0.0570096023
[53,] 0.907757085 0.1844858310 0.0922429155
> postscript(file="/var/www/rcomp/tmp/1wiyt1321982219.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/2rzlv1321982219.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/3u7pj1321982219.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/4szon1321982219.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/51a3y1321982219.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 = 68
Frequency = 1
1 2 3 4 5 6
-7.82594011 -19.21148932 -11.45143268 2.04695181 23.83944811 0.09300767
7 8 9 10 11 12
40.33137502 85.70366290 -55.82326757 -4.63499147 28.99631457 -9.68481864
13 14 15 16 17 18
-38.56092989 -21.08260301 -17.86073480 -3.20471296 33.19183349 -4.86796410
19 20 21 22 23 24
-19.50051008 -7.24240043 -68.48162044 7.15014198 -34.56397233 -36.44132129
25 26 27 28 29 30
42.23815234 15.08517854 -21.14870160 81.92311656 -29.31707144 -6.00652555
31 32 33 34 35 36
-42.62877031 18.03623840 -32.85047633 -12.28659554 -6.28963996 -15.25994070
37 38 39 40 41 42
-16.77191104 -18.61364090 61.18588967 -53.73735522 -11.09814279 42.40334654
43 44 45 46 47 48
168.74506091 180.12748479 -4.89594619 13.49560034 -4.62985858 -17.94912010
49 50 51 52 53 54
-28.55762977 -31.63183493 -37.01769981 11.54965899 24.57586388 -7.26045069
55 56 57 58 59 60
18.87999616 -33.78152625 -22.07717324 -0.25600241 -1.47437455 -11.13175158
61 62 63 64 65 66
18.18196374 34.86342258 -31.62467030 -28.23478850 -0.95087060 -21.39686843
67 68
-11.44623317 -31.87942936
> postscript(file="/var/www/rcomp/tmp/6i6pn1321982219.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 = 68
Frequency = 1
lag(myerror, k = 1) myerror
0 -7.82594011 NA
1 -19.21148932 -7.82594011
2 -11.45143268 -19.21148932
3 2.04695181 -11.45143268
4 23.83944811 2.04695181
5 0.09300767 23.83944811
6 40.33137502 0.09300767
7 85.70366290 40.33137502
8 -55.82326757 85.70366290
9 -4.63499147 -55.82326757
10 28.99631457 -4.63499147
11 -9.68481864 28.99631457
12 -38.56092989 -9.68481864
13 -21.08260301 -38.56092989
14 -17.86073480 -21.08260301
15 -3.20471296 -17.86073480
16 33.19183349 -3.20471296
17 -4.86796410 33.19183349
18 -19.50051008 -4.86796410
19 -7.24240043 -19.50051008
20 -68.48162044 -7.24240043
21 7.15014198 -68.48162044
22 -34.56397233 7.15014198
23 -36.44132129 -34.56397233
24 42.23815234 -36.44132129
25 15.08517854 42.23815234
26 -21.14870160 15.08517854
27 81.92311656 -21.14870160
28 -29.31707144 81.92311656
29 -6.00652555 -29.31707144
30 -42.62877031 -6.00652555
31 18.03623840 -42.62877031
32 -32.85047633 18.03623840
33 -12.28659554 -32.85047633
34 -6.28963996 -12.28659554
35 -15.25994070 -6.28963996
36 -16.77191104 -15.25994070
37 -18.61364090 -16.77191104
38 61.18588967 -18.61364090
39 -53.73735522 61.18588967
40 -11.09814279 -53.73735522
41 42.40334654 -11.09814279
42 168.74506091 42.40334654
43 180.12748479 168.74506091
44 -4.89594619 180.12748479
45 13.49560034 -4.89594619
46 -4.62985858 13.49560034
47 -17.94912010 -4.62985858
48 -28.55762977 -17.94912010
49 -31.63183493 -28.55762977
50 -37.01769981 -31.63183493
51 11.54965899 -37.01769981
52 24.57586388 11.54965899
53 -7.26045069 24.57586388
54 18.87999616 -7.26045069
55 -33.78152625 18.87999616
56 -22.07717324 -33.78152625
57 -0.25600241 -22.07717324
58 -1.47437455 -0.25600241
59 -11.13175158 -1.47437455
60 18.18196374 -11.13175158
61 34.86342258 18.18196374
62 -31.62467030 34.86342258
63 -28.23478850 -31.62467030
64 -0.95087060 -28.23478850
65 -21.39686843 -0.95087060
66 -11.44623317 -21.39686843
67 -31.87942936 -11.44623317
68 NA -31.87942936
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -19.21148932 -7.82594011
[2,] -11.45143268 -19.21148932
[3,] 2.04695181 -11.45143268
[4,] 23.83944811 2.04695181
[5,] 0.09300767 23.83944811
[6,] 40.33137502 0.09300767
[7,] 85.70366290 40.33137502
[8,] -55.82326757 85.70366290
[9,] -4.63499147 -55.82326757
[10,] 28.99631457 -4.63499147
[11,] -9.68481864 28.99631457
[12,] -38.56092989 -9.68481864
[13,] -21.08260301 -38.56092989
[14,] -17.86073480 -21.08260301
[15,] -3.20471296 -17.86073480
[16,] 33.19183349 -3.20471296
[17,] -4.86796410 33.19183349
[18,] -19.50051008 -4.86796410
[19,] -7.24240043 -19.50051008
[20,] -68.48162044 -7.24240043
[21,] 7.15014198 -68.48162044
[22,] -34.56397233 7.15014198
[23,] -36.44132129 -34.56397233
[24,] 42.23815234 -36.44132129
[25,] 15.08517854 42.23815234
[26,] -21.14870160 15.08517854
[27,] 81.92311656 -21.14870160
[28,] -29.31707144 81.92311656
[29,] -6.00652555 -29.31707144
[30,] -42.62877031 -6.00652555
[31,] 18.03623840 -42.62877031
[32,] -32.85047633 18.03623840
[33,] -12.28659554 -32.85047633
[34,] -6.28963996 -12.28659554
[35,] -15.25994070 -6.28963996
[36,] -16.77191104 -15.25994070
[37,] -18.61364090 -16.77191104
[38,] 61.18588967 -18.61364090
[39,] -53.73735522 61.18588967
[40,] -11.09814279 -53.73735522
[41,] 42.40334654 -11.09814279
[42,] 168.74506091 42.40334654
[43,] 180.12748479 168.74506091
[44,] -4.89594619 180.12748479
[45,] 13.49560034 -4.89594619
[46,] -4.62985858 13.49560034
[47,] -17.94912010 -4.62985858
[48,] -28.55762977 -17.94912010
[49,] -31.63183493 -28.55762977
[50,] -37.01769981 -31.63183493
[51,] 11.54965899 -37.01769981
[52,] 24.57586388 11.54965899
[53,] -7.26045069 24.57586388
[54,] 18.87999616 -7.26045069
[55,] -33.78152625 18.87999616
[56,] -22.07717324 -33.78152625
[57,] -0.25600241 -22.07717324
[58,] -1.47437455 -0.25600241
[59,] -11.13175158 -1.47437455
[60,] 18.18196374 -11.13175158
[61,] 34.86342258 18.18196374
[62,] -31.62467030 34.86342258
[63,] -28.23478850 -31.62467030
[64,] -0.95087060 -28.23478850
[65,] -21.39686843 -0.95087060
[66,] -11.44623317 -21.39686843
[67,] -31.87942936 -11.44623317
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -19.21148932 -7.82594011
2 -11.45143268 -19.21148932
3 2.04695181 -11.45143268
4 23.83944811 2.04695181
5 0.09300767 23.83944811
6 40.33137502 0.09300767
7 85.70366290 40.33137502
8 -55.82326757 85.70366290
9 -4.63499147 -55.82326757
10 28.99631457 -4.63499147
11 -9.68481864 28.99631457
12 -38.56092989 -9.68481864
13 -21.08260301 -38.56092989
14 -17.86073480 -21.08260301
15 -3.20471296 -17.86073480
16 33.19183349 -3.20471296
17 -4.86796410 33.19183349
18 -19.50051008 -4.86796410
19 -7.24240043 -19.50051008
20 -68.48162044 -7.24240043
21 7.15014198 -68.48162044
22 -34.56397233 7.15014198
23 -36.44132129 -34.56397233
24 42.23815234 -36.44132129
25 15.08517854 42.23815234
26 -21.14870160 15.08517854
27 81.92311656 -21.14870160
28 -29.31707144 81.92311656
29 -6.00652555 -29.31707144
30 -42.62877031 -6.00652555
31 18.03623840 -42.62877031
32 -32.85047633 18.03623840
33 -12.28659554 -32.85047633
34 -6.28963996 -12.28659554
35 -15.25994070 -6.28963996
36 -16.77191104 -15.25994070
37 -18.61364090 -16.77191104
38 61.18588967 -18.61364090
39 -53.73735522 61.18588967
40 -11.09814279 -53.73735522
41 42.40334654 -11.09814279
42 168.74506091 42.40334654
43 180.12748479 168.74506091
44 -4.89594619 180.12748479
45 13.49560034 -4.89594619
46 -4.62985858 13.49560034
47 -17.94912010 -4.62985858
48 -28.55762977 -17.94912010
49 -31.63183493 -28.55762977
50 -37.01769981 -31.63183493
51 11.54965899 -37.01769981
52 24.57586388 11.54965899
53 -7.26045069 24.57586388
54 18.87999616 -7.26045069
55 -33.78152625 18.87999616
56 -22.07717324 -33.78152625
57 -0.25600241 -22.07717324
58 -1.47437455 -0.25600241
59 -11.13175158 -1.47437455
60 18.18196374 -11.13175158
61 34.86342258 18.18196374
62 -31.62467030 34.86342258
63 -28.23478850 -31.62467030
64 -0.95087060 -28.23478850
65 -21.39686843 -0.95087060
66 -11.44623317 -21.39686843
67 -31.87942936 -11.44623317
> 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/79jvb1321982219.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/8m5gj1321982219.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/9nfde1321982219.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/10fya91321982219.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/115zty1321982219.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/12zil01321982219.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/13sozb1321982219.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/14voz41321982219.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/158ugu1321982219.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/160xen1321982219.tab")
+ }
>
> try(system("convert tmp/1wiyt1321982219.ps tmp/1wiyt1321982219.png",intern=TRUE))
character(0)
> try(system("convert tmp/2rzlv1321982219.ps tmp/2rzlv1321982219.png",intern=TRUE))
character(0)
> try(system("convert tmp/3u7pj1321982219.ps tmp/3u7pj1321982219.png",intern=TRUE))
character(0)
> try(system("convert tmp/4szon1321982219.ps tmp/4szon1321982219.png",intern=TRUE))
character(0)
> try(system("convert tmp/51a3y1321982219.ps tmp/51a3y1321982219.png",intern=TRUE))
character(0)
> try(system("convert tmp/6i6pn1321982219.ps tmp/6i6pn1321982219.png",intern=TRUE))
character(0)
> try(system("convert tmp/79jvb1321982219.ps tmp/79jvb1321982219.png",intern=TRUE))
character(0)
> try(system("convert tmp/8m5gj1321982219.ps tmp/8m5gj1321982219.png",intern=TRUE))
character(0)
> try(system("convert tmp/9nfde1321982219.ps tmp/9nfde1321982219.png",intern=TRUE))
character(0)
> try(system("convert tmp/10fya91321982219.ps tmp/10fya91321982219.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.410 0.400 4.833