R version 2.13.0 (2011-04-13)
Copyright (C) 2011 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(210907
+ ,79
+ ,30
+ ,94
+ ,112285
+ ,144
+ ,145
+ ,120982
+ ,58
+ ,28
+ ,103
+ ,84786
+ ,103
+ ,101
+ ,176508
+ ,60
+ ,38
+ ,93
+ ,83123
+ ,98
+ ,98
+ ,179321
+ ,108
+ ,30
+ ,103
+ ,101193
+ ,135
+ ,132
+ ,123185
+ ,49
+ ,22
+ ,51
+ ,38361
+ ,61
+ ,60
+ ,52746
+ ,0
+ ,26
+ ,70
+ ,68504
+ ,39
+ ,38
+ ,385534
+ ,121
+ ,25
+ ,91
+ ,119182
+ ,150
+ ,144
+ ,33170
+ ,1
+ ,18
+ ,22
+ ,22807
+ ,5
+ ,5
+ ,101645
+ ,20
+ ,11
+ ,38
+ ,17140
+ ,28
+ ,28
+ ,149061
+ ,43
+ ,26
+ ,93
+ ,116174
+ ,84
+ ,84
+ ,165446
+ ,69
+ ,25
+ ,60
+ ,57635
+ ,80
+ ,79
+ ,237213
+ ,78
+ ,38
+ ,123
+ ,66198
+ ,130
+ ,127
+ ,173326
+ ,86
+ ,44
+ ,148
+ ,71701
+ ,82
+ ,78
+ ,133131
+ ,44
+ ,30
+ ,90
+ ,57793
+ ,60
+ ,60
+ ,258873
+ ,104
+ ,40
+ ,124
+ ,80444
+ ,131
+ ,131
+ ,180083
+ ,63
+ ,34
+ ,70
+ ,53855
+ ,84
+ ,84
+ ,324799
+ ,158
+ ,47
+ ,168
+ ,97668
+ ,140
+ ,133
+ ,230964
+ ,102
+ ,30
+ ,115
+ ,133824
+ ,151
+ ,150
+ ,236785
+ ,77
+ ,31
+ ,71
+ ,101481
+ ,91
+ ,91
+ ,135473
+ ,82
+ ,23
+ ,66
+ ,99645
+ ,138
+ ,132
+ ,202925
+ ,115
+ ,36
+ ,134
+ ,114789
+ ,150
+ ,136
+ ,215147
+ ,101
+ ,36
+ ,117
+ ,99052
+ ,124
+ ,124
+ ,344297
+ ,80
+ ,30
+ ,108
+ ,67654
+ ,119
+ ,118
+ ,153935
+ ,50
+ ,25
+ ,84
+ ,65553
+ ,73
+ ,70
+ ,132943
+ ,83
+ ,39
+ ,156
+ ,97500
+ ,110
+ ,107
+ ,174724
+ ,123
+ ,34
+ ,120
+ ,69112
+ ,123
+ ,119
+ ,174415
+ ,73
+ ,31
+ ,114
+ ,82753
+ ,90
+ ,89
+ ,225548
+ ,81
+ ,31
+ ,94
+ ,85323
+ ,116
+ ,112
+ ,223632
+ ,105
+ ,33
+ ,120
+ ,72654
+ ,113
+ ,108
+ ,124817
+ ,47
+ ,25
+ ,81
+ ,30727
+ ,56
+ ,52
+ ,221698
+ ,105
+ ,33
+ ,110
+ ,77873
+ ,115
+ ,112
+ ,210767
+ ,94
+ ,35
+ ,133
+ ,117478
+ ,119
+ ,116
+ ,170266
+ ,44
+ ,42
+ ,122
+ ,74007
+ ,129
+ ,123
+ ,260561
+ ,114
+ ,43
+ ,158
+ ,90183
+ ,127
+ ,125
+ ,84853
+ ,38
+ ,30
+ ,109
+ ,61542
+ ,27
+ ,27
+ ,294424
+ ,107
+ ,33
+ ,124
+ ,101494
+ ,175
+ ,162
+ ,101011
+ ,30
+ ,13
+ ,39
+ ,27570
+ ,35
+ ,32
+ ,215641
+ ,71
+ ,32
+ ,92
+ ,55813
+ ,64
+ ,64
+ ,325107
+ ,84
+ ,36
+ ,126
+ ,79215
+ ,96
+ ,92
+ ,7176
+ ,0
+ ,0
+ ,0
+ ,1423
+ ,0
+ ,0
+ ,167542
+ ,59
+ ,28
+ ,70
+ ,55461
+ ,84
+ ,83
+ ,106408
+ ,33
+ ,14
+ ,37
+ ,31081
+ ,41
+ ,41
+ ,96560
+ ,42
+ ,17
+ ,38
+ ,22996
+ ,47
+ ,47
+ ,265769
+ ,96
+ ,32
+ ,120
+ ,83122
+ ,126
+ ,120
+ ,269651
+ ,106
+ ,30
+ ,93
+ ,70106
+ ,105
+ ,105
+ ,149112
+ ,56
+ ,35
+ ,95
+ ,60578
+ ,80
+ ,79
+ ,175824
+ ,57
+ ,20
+ ,77
+ ,39992
+ ,70
+ ,65
+ ,152871
+ ,59
+ ,28
+ ,90
+ ,79892
+ ,73
+ ,70)
+ ,dim=c(7
+ ,48)
+ ,dimnames=list(c('Y'
+ ,'X1'
+ ,'X2'
+ ,'X3'
+ ,'X4'
+ ,'X5'
+ ,'X6')
+ ,1:48))
> y <- array(NA,dim=c(7,48),dimnames=list(c('Y','X1','X2','X3','X4','X5','X6'),1:48))
> 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 = '6'
> 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
X5 Y X1 X2 X3 X4 X6
1 144 210907 79 30 94 112285 145
2 103 120982 58 28 103 84786 101
3 98 176508 60 38 93 83123 98
4 135 179321 108 30 103 101193 132
5 61 123185 49 22 51 38361 60
6 39 52746 0 26 70 68504 38
7 150 385534 121 25 91 119182 144
8 5 33170 1 18 22 22807 5
9 28 101645 20 11 38 17140 28
10 84 149061 43 26 93 116174 84
11 80 165446 69 25 60 57635 79
12 130 237213 78 38 123 66198 127
13 82 173326 86 44 148 71701 78
14 60 133131 44 30 90 57793 60
15 131 258873 104 40 124 80444 131
16 84 180083 63 34 70 53855 84
17 140 324799 158 47 168 97668 133
18 151 230964 102 30 115 133824 150
19 91 236785 77 31 71 101481 91
20 138 135473 82 23 66 99645 132
21 150 202925 115 36 134 114789 136
22 124 215147 101 36 117 99052 124
23 119 344297 80 30 108 67654 118
24 73 153935 50 25 84 65553 70
25 110 132943 83 39 156 97500 107
26 123 174724 123 34 120 69112 119
27 90 174415 73 31 114 82753 89
28 116 225548 81 31 94 85323 112
29 113 223632 105 33 120 72654 108
30 56 124817 47 25 81 30727 52
31 115 221698 105 33 110 77873 112
32 119 210767 94 35 133 117478 116
33 129 170266 44 42 122 74007 123
34 127 260561 114 43 158 90183 125
35 27 84853 38 30 109 61542 27
36 175 294424 107 33 124 101494 162
37 35 101011 30 13 39 27570 32
38 64 215641 71 32 92 55813 64
39 96 325107 84 36 126 79215 92
40 0 7176 0 0 0 1423 0
41 84 167542 59 28 70 55461 83
42 41 106408 33 14 37 31081 41
43 47 96560 42 17 38 22996 47
44 126 265769 96 32 120 83122 120
45 105 269651 106 30 93 70106 105
46 80 149112 56 35 95 60578 79
47 70 175824 57 20 77 39992 65
48 73 152871 59 28 90 79892 70
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Y X1 X2 X3 X4
8.224e-01 -4.450e-06 1.059e-02 -2.219e-01 6.679e-02 -2.921e-05
X6
1.045e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4.6160 -1.4121 -0.3749 1.5141 9.1065
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.224e-01 1.479e+00 0.556 0.5812
Y -4.450e-06 8.984e-06 -0.495 0.6230
X1 1.059e-02 2.759e-02 0.384 0.7030
X2 -2.219e-01 1.050e-01 -2.114 0.0407 *
X3 6.679e-02 3.012e-02 2.218 0.0322 *
X4 -2.921e-05 2.529e-05 -1.155 0.2547
X6 1.045e+00 2.636e-02 39.652 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.726 on 41 degrees of freedom
Multiple R-squared: 0.9961, Adjusted R-squared: 0.9956
F-statistic: 1767 on 6 and 41 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.22523940 0.45047880 0.77476060
[2,] 0.11693417 0.23386835 0.88306583
[3,] 0.05302719 0.10605438 0.94697281
[4,] 0.04309243 0.08618486 0.95690757
[5,] 0.03228572 0.06457145 0.96771428
[6,] 0.05642187 0.11284374 0.94357813
[7,] 0.02824429 0.05648857 0.97175571
[8,] 0.01603648 0.03207296 0.98396352
[9,] 0.02084085 0.04168169 0.97915915
[10,] 0.01391958 0.02783916 0.98608042
[11,] 0.09442442 0.18884885 0.90557558
[12,] 0.87839303 0.24321393 0.12160697
[13,] 0.91860737 0.16278526 0.08139263
[14,] 0.98715821 0.02568359 0.01284179
[15,] 0.97697842 0.04604316 0.02302158
[16,] 0.97173922 0.05652157 0.02826078
[17,] 0.96522594 0.06954811 0.03477406
[18,] 0.97113208 0.05773584 0.02886792
[19,] 0.95411938 0.09176124 0.04588062
[20,] 0.93867263 0.12265474 0.06132737
[21,] 0.94958429 0.10083141 0.05041571
[22,] 0.92776567 0.14446867 0.07223433
[23,] 0.94965982 0.10068036 0.05034018
[24,] 0.96716718 0.06566563 0.03283282
[25,] 0.96552216 0.06895567 0.03447784
[26,] 0.93381969 0.13236063 0.06618031
[27,] 0.97814611 0.04370777 0.02185389
[28,] 0.98483386 0.03033228 0.01516614
[29,] 0.94594598 0.10810804 0.05405402
> postscript(file="/var/wessaorg/rcomp/tmp/1wgs21324629049.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/wessaorg/rcomp/tmp/2znhh1324629049.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/wessaorg/rcomp/tmp/3igzo1324629049.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/wessaorg/rcomp/tmp/4yb2y1324629049.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/wessaorg/rcomp/tmp/5m3341324629049.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 = 48
Frequency = 1
1 2 3 4 5 6
-4.61598468 -1.65351600 -0.45397714 -1.40137970 0.09047536 1.78902567
7 8 9 10 11 12
2.05405378 2.27900313 -1.44419954 -1.46035812 -0.16445105 -1.18323107
13 14 15 16 17 18
1.48414310 -1.07475383 -3.74990229 -0.04328521 1.99839475 -3.76980398
19 20 21 22 23 24
0.40323626 2.55198102 9.10646855 -3.47274169 -3.05160147 1.02048698
25 26 27 28 29 30
-1.86442976 -1.17837766 -2.16082439 1.35358464 0.60862661 1.91920838
31 32 33 34 35 36
-0.76034586 -0.80891749 3.24206080 -2.89810243 -0.89435336 7.03743927
37 38 39 40 41 42
1.94790535 -0.92197318 1.46218815 -0.74887944 -0.29582704 -1.00855424
43 44 45 46 47 48
-1.05621362 1.43309476 -2.99815339 -0.13264788 1.88129943 1.60410950
> postscript(file="/var/wessaorg/rcomp/tmp/60f2o1324629049.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 = 48
Frequency = 1
lag(myerror, k = 1) myerror
0 -4.61598468 NA
1 -1.65351600 -4.61598468
2 -0.45397714 -1.65351600
3 -1.40137970 -0.45397714
4 0.09047536 -1.40137970
5 1.78902567 0.09047536
6 2.05405378 1.78902567
7 2.27900313 2.05405378
8 -1.44419954 2.27900313
9 -1.46035812 -1.44419954
10 -0.16445105 -1.46035812
11 -1.18323107 -0.16445105
12 1.48414310 -1.18323107
13 -1.07475383 1.48414310
14 -3.74990229 -1.07475383
15 -0.04328521 -3.74990229
16 1.99839475 -0.04328521
17 -3.76980398 1.99839475
18 0.40323626 -3.76980398
19 2.55198102 0.40323626
20 9.10646855 2.55198102
21 -3.47274169 9.10646855
22 -3.05160147 -3.47274169
23 1.02048698 -3.05160147
24 -1.86442976 1.02048698
25 -1.17837766 -1.86442976
26 -2.16082439 -1.17837766
27 1.35358464 -2.16082439
28 0.60862661 1.35358464
29 1.91920838 0.60862661
30 -0.76034586 1.91920838
31 -0.80891749 -0.76034586
32 3.24206080 -0.80891749
33 -2.89810243 3.24206080
34 -0.89435336 -2.89810243
35 7.03743927 -0.89435336
36 1.94790535 7.03743927
37 -0.92197318 1.94790535
38 1.46218815 -0.92197318
39 -0.74887944 1.46218815
40 -0.29582704 -0.74887944
41 -1.00855424 -0.29582704
42 -1.05621362 -1.00855424
43 1.43309476 -1.05621362
44 -2.99815339 1.43309476
45 -0.13264788 -2.99815339
46 1.88129943 -0.13264788
47 1.60410950 1.88129943
48 NA 1.60410950
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.65351600 -4.61598468
[2,] -0.45397714 -1.65351600
[3,] -1.40137970 -0.45397714
[4,] 0.09047536 -1.40137970
[5,] 1.78902567 0.09047536
[6,] 2.05405378 1.78902567
[7,] 2.27900313 2.05405378
[8,] -1.44419954 2.27900313
[9,] -1.46035812 -1.44419954
[10,] -0.16445105 -1.46035812
[11,] -1.18323107 -0.16445105
[12,] 1.48414310 -1.18323107
[13,] -1.07475383 1.48414310
[14,] -3.74990229 -1.07475383
[15,] -0.04328521 -3.74990229
[16,] 1.99839475 -0.04328521
[17,] -3.76980398 1.99839475
[18,] 0.40323626 -3.76980398
[19,] 2.55198102 0.40323626
[20,] 9.10646855 2.55198102
[21,] -3.47274169 9.10646855
[22,] -3.05160147 -3.47274169
[23,] 1.02048698 -3.05160147
[24,] -1.86442976 1.02048698
[25,] -1.17837766 -1.86442976
[26,] -2.16082439 -1.17837766
[27,] 1.35358464 -2.16082439
[28,] 0.60862661 1.35358464
[29,] 1.91920838 0.60862661
[30,] -0.76034586 1.91920838
[31,] -0.80891749 -0.76034586
[32,] 3.24206080 -0.80891749
[33,] -2.89810243 3.24206080
[34,] -0.89435336 -2.89810243
[35,] 7.03743927 -0.89435336
[36,] 1.94790535 7.03743927
[37,] -0.92197318 1.94790535
[38,] 1.46218815 -0.92197318
[39,] -0.74887944 1.46218815
[40,] -0.29582704 -0.74887944
[41,] -1.00855424 -0.29582704
[42,] -1.05621362 -1.00855424
[43,] 1.43309476 -1.05621362
[44,] -2.99815339 1.43309476
[45,] -0.13264788 -2.99815339
[46,] 1.88129943 -0.13264788
[47,] 1.60410950 1.88129943
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.65351600 -4.61598468
2 -0.45397714 -1.65351600
3 -1.40137970 -0.45397714
4 0.09047536 -1.40137970
5 1.78902567 0.09047536
6 2.05405378 1.78902567
7 2.27900313 2.05405378
8 -1.44419954 2.27900313
9 -1.46035812 -1.44419954
10 -0.16445105 -1.46035812
11 -1.18323107 -0.16445105
12 1.48414310 -1.18323107
13 -1.07475383 1.48414310
14 -3.74990229 -1.07475383
15 -0.04328521 -3.74990229
16 1.99839475 -0.04328521
17 -3.76980398 1.99839475
18 0.40323626 -3.76980398
19 2.55198102 0.40323626
20 9.10646855 2.55198102
21 -3.47274169 9.10646855
22 -3.05160147 -3.47274169
23 1.02048698 -3.05160147
24 -1.86442976 1.02048698
25 -1.17837766 -1.86442976
26 -2.16082439 -1.17837766
27 1.35358464 -2.16082439
28 0.60862661 1.35358464
29 1.91920838 0.60862661
30 -0.76034586 1.91920838
31 -0.80891749 -0.76034586
32 3.24206080 -0.80891749
33 -2.89810243 3.24206080
34 -0.89435336 -2.89810243
35 7.03743927 -0.89435336
36 1.94790535 7.03743927
37 -0.92197318 1.94790535
38 1.46218815 -0.92197318
39 -0.74887944 1.46218815
40 -0.29582704 -0.74887944
41 -1.00855424 -0.29582704
42 -1.05621362 -1.00855424
43 1.43309476 -1.05621362
44 -2.99815339 1.43309476
45 -0.13264788 -2.99815339
46 1.88129943 -0.13264788
47 1.60410950 1.88129943
> 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/wessaorg/rcomp/tmp/7xsx01324629049.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/wessaorg/rcomp/tmp/8lm4i1324629049.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/wessaorg/rcomp/tmp/92xk11324629049.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/wessaorg/rcomp/tmp/10l1j31324629049.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11iirg1324629049.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/wessaorg/rcomp/tmp/12uufh1324629049.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/wessaorg/rcomp/tmp/13yily1324629049.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/wessaorg/rcomp/tmp/14xomv1324629049.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/wessaorg/rcomp/tmp/15hjw21324629049.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/wessaorg/rcomp/tmp/1661k81324629050.tab")
+ }
>
> try(system("convert tmp/1wgs21324629049.ps tmp/1wgs21324629049.png",intern=TRUE))
character(0)
> try(system("convert tmp/2znhh1324629049.ps tmp/2znhh1324629049.png",intern=TRUE))
character(0)
> try(system("convert tmp/3igzo1324629049.ps tmp/3igzo1324629049.png",intern=TRUE))
character(0)
> try(system("convert tmp/4yb2y1324629049.ps tmp/4yb2y1324629049.png",intern=TRUE))
character(0)
> try(system("convert tmp/5m3341324629049.ps tmp/5m3341324629049.png",intern=TRUE))
character(0)
> try(system("convert tmp/60f2o1324629049.ps tmp/60f2o1324629049.png",intern=TRUE))
character(0)
> try(system("convert tmp/7xsx01324629049.ps tmp/7xsx01324629049.png",intern=TRUE))
character(0)
> try(system("convert tmp/8lm4i1324629049.ps tmp/8lm4i1324629049.png",intern=TRUE))
character(0)
> try(system("convert tmp/92xk11324629049.ps tmp/92xk11324629049.png",intern=TRUE))
character(0)
> try(system("convert tmp/10l1j31324629049.ps tmp/10l1j31324629049.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.038 0.589 3.879