R version 2.12.1 (2010-12-16)
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(13
+ ,2528
+ ,80
+ ,15.3
+ ,12
+ ,3333
+ ,83
+ ,19.4
+ ,54
+ ,19611
+ ,96
+ ,19.5
+ ,19
+ ,3570
+ ,56
+ ,17
+ ,37
+ ,1722
+ ,43
+ ,19.3
+ ,2
+ ,583
+ ,51
+ ,12.9
+ ,72
+ ,4790
+ ,91
+ ,16.7
+ ,164
+ ,35971
+ ,81
+ ,13.8
+ ,18
+ ,25440
+ ,120
+ ,13.7
+ ,1
+ ,2217
+ ,46
+ ,14.3
+ ,53
+ ,1971
+ ,56
+ ,22.2
+ ,16
+ ,12620
+ ,37
+ ,16.8
+ ,32
+ ,19046
+ ,120
+ ,18
+ ,21
+ ,8612
+ ,103
+ ,15
+ ,23
+ ,3896
+ ,105
+ ,18.4
+ ,18
+ ,6298
+ ,42
+ ,18.2
+ ,112
+ ,27350
+ ,65
+ ,24.1
+ ,25
+ ,4145
+ ,51
+ ,23
+ ,5
+ ,1175
+ ,57
+ ,21.8
+ ,26
+ ,8297
+ ,60
+ ,19.1
+ ,8
+ ,7814
+ ,160
+ ,22.6
+ ,15
+ ,1745
+ ,48
+ ,14.3
+ ,11
+ ,5046
+ ,109
+ ,19
+ ,11
+ ,18943
+ ,50
+ ,18.5
+ ,87
+ ,8624
+ ,78
+ ,21.3
+ ,33
+ ,2225
+ ,41
+ ,20.5
+ ,22
+ ,12659
+ ,65
+ ,18
+ ,98
+ ,1967
+ ,50
+ ,16.8
+ ,1
+ ,1172
+ ,73
+ ,20.5
+ ,5
+ ,639
+ ,26
+ ,20.1
+ ,1
+ ,7056
+ ,60
+ ,24.5
+ ,38
+ ,1934
+ ,85
+ ,12
+ ,30
+ ,6260
+ ,133
+ ,21
+ ,12
+ ,424
+ ,62
+ ,20.2
+ ,24
+ ,3488
+ ,44
+ ,24
+ ,6
+ ,3330
+ ,67
+ ,14.9
+ ,15
+ ,2227
+ ,54
+ ,24
+ ,38
+ ,8115
+ ,110
+ ,20.5
+ ,84
+ ,1600
+ ,56
+ ,19.5
+ ,3
+ ,15305
+ ,85
+ ,17.5
+ ,18
+ ,7121
+ ,58
+ ,16
+ ,63
+ ,5794
+ ,34
+ ,17.5
+ ,239
+ ,8636
+ ,150
+ ,18.1
+ ,234
+ ,4803
+ ,93
+ ,24.3
+ ,6
+ ,1097
+ ,53
+ ,13.1
+ ,76
+ ,9765
+ ,130
+ ,16.9
+ ,25
+ ,4266
+ ,68
+ ,17
+ ,8
+ ,1507
+ ,51
+ ,21
+ ,23
+ ,3836
+ ,121
+ ,18.5
+ ,16
+ ,17419
+ ,48
+ ,18
+ ,6
+ ,8735
+ ,63
+ ,20.8
+ ,100
+ ,22550
+ ,107
+ ,23
+ ,80
+ ,9961
+ ,79
+ ,21.8
+ ,28
+ ,4706
+ ,61
+ ,19.7
+ ,48
+ ,4011
+ ,52
+ ,18.9
+ ,18
+ ,6949
+ ,100
+ ,18.6
+ ,36
+ ,11405
+ ,70
+ ,23.6
+ ,19
+ ,904
+ ,39
+ ,19.2
+ ,32
+ ,3332
+ ,73
+ ,17.7
+ ,3
+ ,575
+ ,33
+ ,18
+ ,106
+ ,29708
+ ,73
+ ,21.4
+ ,62
+ ,2511
+ ,60
+ ,17.7
+ ,23
+ ,18422
+ ,45
+ ,20.1
+ ,2
+ ,6311
+ ,46
+ ,18.5
+ ,26
+ ,1450
+ ,60
+ ,18.6
+ ,20
+ ,4106
+ ,96
+ ,15.4
+ ,38
+ ,10274
+ ,90
+ ,15
+ ,19
+ ,510
+ ,82
+ ,26.5)
+ ,dim=c(4
+ ,68)
+ ,dimnames=list(c('fish'
+ ,'acre'
+ ,'depth'
+ ,'temp')
+ ,1:68))
> y <- array(NA,dim=c(4,68),dimnames=list(c('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 = '1'
> 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 acre depth temp
1 13 2528 80 15.3
2 12 3333 83 19.4
3 54 19611 96 19.5
4 19 3570 56 17.0
5 37 1722 43 19.3
6 2 583 51 12.9
7 72 4790 91 16.7
8 164 35971 81 13.8
9 18 25440 120 13.7
10 1 2217 46 14.3
11 53 1971 56 22.2
12 16 12620 37 16.8
13 32 19046 120 18.0
14 21 8612 103 15.0
15 23 3896 105 18.4
16 18 6298 42 18.2
17 112 27350 65 24.1
18 25 4145 51 23.0
19 5 1175 57 21.8
20 26 8297 60 19.1
21 8 7814 160 22.6
22 15 1745 48 14.3
23 11 5046 109 19.0
24 11 18943 50 18.5
25 87 8624 78 21.3
26 33 2225 41 20.5
27 22 12659 65 18.0
28 98 1967 50 16.8
29 1 1172 73 20.5
30 5 639 26 20.1
31 1 7056 60 24.5
32 38 1934 85 12.0
33 30 6260 133 21.0
34 12 424 62 20.2
35 24 3488 44 24.0
36 6 3330 67 14.9
37 15 2227 54 24.0
38 38 8115 110 20.5
39 84 1600 56 19.5
40 3 15305 85 17.5
41 18 7121 58 16.0
42 63 5794 34 17.5
43 239 8636 150 18.1
44 234 4803 93 24.3
45 6 1097 53 13.1
46 76 9765 130 16.9
47 25 4266 68 17.0
48 8 1507 51 21.0
49 23 3836 121 18.5
50 16 17419 48 18.0
51 6 8735 63 20.8
52 100 22550 107 23.0
53 80 9961 79 21.8
54 28 4706 61 19.7
55 48 4011 52 18.9
56 18 6949 100 18.6
57 36 11405 70 23.6
58 19 904 39 19.2
59 32 3332 73 17.7
60 3 575 33 18.0
61 106 29708 73 21.4
62 62 2511 60 17.7
63 23 18422 45 20.1
64 2 6311 46 18.5
65 26 1450 60 18.6
66 20 4106 96 15.4
67 38 10274 90 15.0
68 19 510 82 26.5
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) acre depth temp
-41.239164 0.001745 0.373837 2.105029
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-71.786 -23.236 -8.514 14.106 180.938
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.124e+01 3.505e+01 -1.177 0.2437
acre 1.745e-03 7.107e-04 2.456 0.0168 *
depth 3.738e-01 1.880e-01 1.988 0.0510 .
temp 2.105e+00 1.692e+00 1.244 0.2181
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 43.56 on 64 degrees of freedom
Multiple R-squared: 0.1894, Adjusted R-squared: 0.1514
F-statistic: 4.984 on 3 and 64 DF, p-value: 0.003605
> 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.242424779 0.4848495571 7.575752e-01
[2,] 0.244951265 0.4899025295 7.550487e-01
[3,] 0.414471792 0.8289435846 5.855282e-01
[4,] 0.329183124 0.6583662471 6.708169e-01
[5,] 0.226714986 0.4534299721 7.732850e-01
[6,] 0.276470107 0.5529402137 7.235299e-01
[7,] 0.221431701 0.4428634012 7.785683e-01
[8,] 0.152915429 0.3058308572 8.470846e-01
[9,] 0.104104792 0.2082095840 8.958952e-01
[10,] 0.074561482 0.1491229636 9.254385e-01
[11,] 0.048433502 0.0968670045 9.515665e-01
[12,] 0.030978038 0.0619560753 9.690220e-01
[13,] 0.020583941 0.0411678811 9.794161e-01
[14,] 0.012122794 0.0242455873 9.878772e-01
[15,] 0.011093591 0.0221871823 9.889064e-01
[16,] 0.006083825 0.0121676501 9.939162e-01
[17,] 0.004015668 0.0080313362 9.959843e-01
[18,] 0.009168437 0.0183368748 9.908316e-01
[19,] 0.013610825 0.0272216497 9.863892e-01
[20,] 0.008181767 0.0163635334 9.918182e-01
[21,] 0.005665452 0.0113309043 9.943345e-01
[22,] 0.032931449 0.0658628988 9.670686e-01
[23,] 0.025233355 0.0504667097 9.747666e-01
[24,] 0.017737726 0.0354754518 9.822623e-01
[25,] 0.018154528 0.0363090554 9.818455e-01
[26,] 0.013433372 0.0268667433 9.865666e-01
[27,] 0.014108242 0.0282164831 9.858918e-01
[28,] 0.009054016 0.0181080318 9.909460e-01
[29,] 0.005458831 0.0109176620 9.945412e-01
[30,] 0.003397499 0.0067949988 9.966025e-01
[31,] 0.002199644 0.0043992871 9.978004e-01
[32,] 0.001822762 0.0036455250 9.981772e-01
[33,] 0.004408641 0.0088172818 9.955914e-01
[34,] 0.006126631 0.0122532623 9.938734e-01
[35,] 0.003704930 0.0074098606 9.962951e-01
[36,] 0.005018448 0.0100368962 9.949816e-01
[37,] 0.541983431 0.9160331374 4.580166e-01
[38,] 0.999947050 0.0001058996 5.294978e-05
[39,] 0.999869322 0.0002613568 1.306784e-04
[40,] 0.999765355 0.0004692891 2.346446e-04
[41,] 0.999454276 0.0010914475 5.457237e-04
[42,] 0.998858809 0.0022823814 1.141191e-03
[43,] 0.998146304 0.0037073920 1.853696e-03
[44,] 0.997840087 0.0043198265 2.159913e-03
[45,] 0.997856686 0.0042866279 2.143314e-03
[46,] 0.996031928 0.0079361444 3.968072e-03
[47,] 0.996945321 0.0061093577 3.054679e-03
[48,] 0.992838958 0.0143220831 7.161042e-03
[49,] 0.990275819 0.0194483623 9.724181e-03
[50,] 0.984829892 0.0303402161 1.517011e-02
[51,] 0.967758145 0.0644837091 3.224185e-02
[52,] 0.935642306 0.1287153878 6.435769e-02
[53,] 0.874157269 0.2516854624 1.258427e-01
[54,] 0.766999336 0.4660013289 2.330007e-01
[55,] 0.820859578 0.3582808441 1.791404e-01
> postscript(file="/var/www/rcomp/tmp/13myr1321913885.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/25kkl1321913885.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/3pqkt1321913885.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/4hpce1321913885.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/536rc1321913885.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
-12.286810 -24.443882 -15.923734 -2.711831 18.531741 -3.998910 35.706143
8 9 10 11 12 13 14
83.129919 -58.859824 -8.928532 23.132699 -13.982585 -42.752211 -22.871767
15 16 17 18 19 20 21
-20.545860 -5.765233 30.475466 -8.476349 -23.009894 -9.877621 -71.785997
22 23 24 25 26 27 28
5.147560 -37.311285 -38.456347 39.191539 11.875512 -21.044132 81.749863
29 30 31 32 33 34 35
-30.249518 -6.906922 -44.078902 18.827289 -33.612195 -13.200338 -7.817876
36 37 38 39 40 41 42
-14.984617 -18.355467 -19.198892 60.463772 -51.086343 -8.551921 44.578602
43 44 45 46 47 48 49
170.990399 180.937558 -2.064658 16.022775 -2.412586 -16.662275 -26.633045
50 51 52 53 54 55 56
-28.996370 -35.342110 13.467112 28.431765 -3.247221 23.014300 -29.425981
57 58 59 60 61 62 63
-18.512918 3.665223 2.874787 -6.991525 23.053005 39.167537 -27.045922
64 65 66 67 68
-23.914779 3.124739 -14.232742 -3.912516 -27.088863
> postscript(file="/var/www/rcomp/tmp/6jehp1321913885.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 -12.286810 NA
1 -24.443882 -12.286810
2 -15.923734 -24.443882
3 -2.711831 -15.923734
4 18.531741 -2.711831
5 -3.998910 18.531741
6 35.706143 -3.998910
7 83.129919 35.706143
8 -58.859824 83.129919
9 -8.928532 -58.859824
10 23.132699 -8.928532
11 -13.982585 23.132699
12 -42.752211 -13.982585
13 -22.871767 -42.752211
14 -20.545860 -22.871767
15 -5.765233 -20.545860
16 30.475466 -5.765233
17 -8.476349 30.475466
18 -23.009894 -8.476349
19 -9.877621 -23.009894
20 -71.785997 -9.877621
21 5.147560 -71.785997
22 -37.311285 5.147560
23 -38.456347 -37.311285
24 39.191539 -38.456347
25 11.875512 39.191539
26 -21.044132 11.875512
27 81.749863 -21.044132
28 -30.249518 81.749863
29 -6.906922 -30.249518
30 -44.078902 -6.906922
31 18.827289 -44.078902
32 -33.612195 18.827289
33 -13.200338 -33.612195
34 -7.817876 -13.200338
35 -14.984617 -7.817876
36 -18.355467 -14.984617
37 -19.198892 -18.355467
38 60.463772 -19.198892
39 -51.086343 60.463772
40 -8.551921 -51.086343
41 44.578602 -8.551921
42 170.990399 44.578602
43 180.937558 170.990399
44 -2.064658 180.937558
45 16.022775 -2.064658
46 -2.412586 16.022775
47 -16.662275 -2.412586
48 -26.633045 -16.662275
49 -28.996370 -26.633045
50 -35.342110 -28.996370
51 13.467112 -35.342110
52 28.431765 13.467112
53 -3.247221 28.431765
54 23.014300 -3.247221
55 -29.425981 23.014300
56 -18.512918 -29.425981
57 3.665223 -18.512918
58 2.874787 3.665223
59 -6.991525 2.874787
60 23.053005 -6.991525
61 39.167537 23.053005
62 -27.045922 39.167537
63 -23.914779 -27.045922
64 3.124739 -23.914779
65 -14.232742 3.124739
66 -3.912516 -14.232742
67 -27.088863 -3.912516
68 NA -27.088863
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -24.443882 -12.286810
[2,] -15.923734 -24.443882
[3,] -2.711831 -15.923734
[4,] 18.531741 -2.711831
[5,] -3.998910 18.531741
[6,] 35.706143 -3.998910
[7,] 83.129919 35.706143
[8,] -58.859824 83.129919
[9,] -8.928532 -58.859824
[10,] 23.132699 -8.928532
[11,] -13.982585 23.132699
[12,] -42.752211 -13.982585
[13,] -22.871767 -42.752211
[14,] -20.545860 -22.871767
[15,] -5.765233 -20.545860
[16,] 30.475466 -5.765233
[17,] -8.476349 30.475466
[18,] -23.009894 -8.476349
[19,] -9.877621 -23.009894
[20,] -71.785997 -9.877621
[21,] 5.147560 -71.785997
[22,] -37.311285 5.147560
[23,] -38.456347 -37.311285
[24,] 39.191539 -38.456347
[25,] 11.875512 39.191539
[26,] -21.044132 11.875512
[27,] 81.749863 -21.044132
[28,] -30.249518 81.749863
[29,] -6.906922 -30.249518
[30,] -44.078902 -6.906922
[31,] 18.827289 -44.078902
[32,] -33.612195 18.827289
[33,] -13.200338 -33.612195
[34,] -7.817876 -13.200338
[35,] -14.984617 -7.817876
[36,] -18.355467 -14.984617
[37,] -19.198892 -18.355467
[38,] 60.463772 -19.198892
[39,] -51.086343 60.463772
[40,] -8.551921 -51.086343
[41,] 44.578602 -8.551921
[42,] 170.990399 44.578602
[43,] 180.937558 170.990399
[44,] -2.064658 180.937558
[45,] 16.022775 -2.064658
[46,] -2.412586 16.022775
[47,] -16.662275 -2.412586
[48,] -26.633045 -16.662275
[49,] -28.996370 -26.633045
[50,] -35.342110 -28.996370
[51,] 13.467112 -35.342110
[52,] 28.431765 13.467112
[53,] -3.247221 28.431765
[54,] 23.014300 -3.247221
[55,] -29.425981 23.014300
[56,] -18.512918 -29.425981
[57,] 3.665223 -18.512918
[58,] 2.874787 3.665223
[59,] -6.991525 2.874787
[60,] 23.053005 -6.991525
[61,] 39.167537 23.053005
[62,] -27.045922 39.167537
[63,] -23.914779 -27.045922
[64,] 3.124739 -23.914779
[65,] -14.232742 3.124739
[66,] -3.912516 -14.232742
[67,] -27.088863 -3.912516
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -24.443882 -12.286810
2 -15.923734 -24.443882
3 -2.711831 -15.923734
4 18.531741 -2.711831
5 -3.998910 18.531741
6 35.706143 -3.998910
7 83.129919 35.706143
8 -58.859824 83.129919
9 -8.928532 -58.859824
10 23.132699 -8.928532
11 -13.982585 23.132699
12 -42.752211 -13.982585
13 -22.871767 -42.752211
14 -20.545860 -22.871767
15 -5.765233 -20.545860
16 30.475466 -5.765233
17 -8.476349 30.475466
18 -23.009894 -8.476349
19 -9.877621 -23.009894
20 -71.785997 -9.877621
21 5.147560 -71.785997
22 -37.311285 5.147560
23 -38.456347 -37.311285
24 39.191539 -38.456347
25 11.875512 39.191539
26 -21.044132 11.875512
27 81.749863 -21.044132
28 -30.249518 81.749863
29 -6.906922 -30.249518
30 -44.078902 -6.906922
31 18.827289 -44.078902
32 -33.612195 18.827289
33 -13.200338 -33.612195
34 -7.817876 -13.200338
35 -14.984617 -7.817876
36 -18.355467 -14.984617
37 -19.198892 -18.355467
38 60.463772 -19.198892
39 -51.086343 60.463772
40 -8.551921 -51.086343
41 44.578602 -8.551921
42 170.990399 44.578602
43 180.937558 170.990399
44 -2.064658 180.937558
45 16.022775 -2.064658
46 -2.412586 16.022775
47 -16.662275 -2.412586
48 -26.633045 -16.662275
49 -28.996370 -26.633045
50 -35.342110 -28.996370
51 13.467112 -35.342110
52 28.431765 13.467112
53 -3.247221 28.431765
54 23.014300 -3.247221
55 -29.425981 23.014300
56 -18.512918 -29.425981
57 3.665223 -18.512918
58 2.874787 3.665223
59 -6.991525 2.874787
60 23.053005 -6.991525
61 39.167537 23.053005
62 -27.045922 39.167537
63 -23.914779 -27.045922
64 3.124739 -23.914779
65 -14.232742 3.124739
66 -3.912516 -14.232742
67 -27.088863 -3.912516
> 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/7ipy11321913885.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/8j04n1321913885.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/9n3fs1321913885.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/10o0ee1321913885.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/11qz781321913885.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/12okms1321913885.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/13qe4a1321913885.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/14sl0b1321913885.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/15bzdu1321913885.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/160hin1321913885.tab")
+ }
>
> try(system("convert tmp/13myr1321913885.ps tmp/13myr1321913885.png",intern=TRUE))
character(0)
> try(system("convert tmp/25kkl1321913885.ps tmp/25kkl1321913885.png",intern=TRUE))
character(0)
> try(system("convert tmp/3pqkt1321913885.ps tmp/3pqkt1321913885.png",intern=TRUE))
character(0)
> try(system("convert tmp/4hpce1321913885.ps tmp/4hpce1321913885.png",intern=TRUE))
character(0)
> try(system("convert tmp/536rc1321913885.ps tmp/536rc1321913885.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jehp1321913885.ps tmp/6jehp1321913885.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ipy11321913885.ps tmp/7ipy11321913885.png",intern=TRUE))
character(0)
> try(system("convert tmp/8j04n1321913885.ps tmp/8j04n1321913885.png",intern=TRUE))
character(0)
> try(system("convert tmp/9n3fs1321913885.ps tmp/9n3fs1321913885.png",intern=TRUE))
character(0)
> try(system("convert tmp/10o0ee1321913885.ps tmp/10o0ee1321913885.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.368 0.620 4.961