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 = '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
Y X1 X2 X3 X4 X5 X6
1 210907 79 30 94 112285 144 145
2 120982 58 28 103 84786 103 101
3 176508 60 38 93 83123 98 98
4 179321 108 30 103 101193 135 132
5 123185 49 22 51 38361 61 60
6 52746 0 26 70 68504 39 38
7 385534 121 25 91 119182 150 144
8 33170 1 18 22 22807 5 5
9 101645 20 11 38 17140 28 28
10 149061 43 26 93 116174 84 84
11 165446 69 25 60 57635 80 79
12 237213 78 38 123 66198 130 127
13 173326 86 44 148 71701 82 78
14 133131 44 30 90 57793 60 60
15 258873 104 40 124 80444 131 131
16 180083 63 34 70 53855 84 84
17 324799 158 47 168 97668 140 133
18 230964 102 30 115 133824 151 150
19 236785 77 31 71 101481 91 91
20 135473 82 23 66 99645 138 132
21 202925 115 36 134 114789 150 136
22 215147 101 36 117 99052 124 124
23 344297 80 30 108 67654 119 118
24 153935 50 25 84 65553 73 70
25 132943 83 39 156 97500 110 107
26 174724 123 34 120 69112 123 119
27 174415 73 31 114 82753 90 89
28 225548 81 31 94 85323 116 112
29 223632 105 33 120 72654 113 108
30 124817 47 25 81 30727 56 52
31 221698 105 33 110 77873 115 112
32 210767 94 35 133 117478 119 116
33 170266 44 42 122 74007 129 123
34 260561 114 43 158 90183 127 125
35 84853 38 30 109 61542 27 27
36 294424 107 33 124 101494 175 162
37 101011 30 13 39 27570 35 32
38 215641 71 32 92 55813 64 64
39 325107 84 36 126 79215 96 92
40 7176 0 0 0 1423 0 0
41 167542 59 28 70 55461 84 83
42 106408 33 14 37 31081 41 41
43 96560 42 17 38 22996 47 47
44 265769 96 32 120 83122 126 120
45 269651 106 30 93 70106 105 105
46 149112 56 35 95 60578 80 79
47 175824 57 20 77 39992 70 65
48 152871 59 28 90 79892 73 70
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X1 X2 X3 X4 X5
36729.7493 1330.0799 516.0858 -96.4629 -0.2918 -1337.0089
X6
2108.2066
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-97834 -19196 -2073 15799 127739
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 36729.7493 25085.8886 1.464 0.15078
X1 1330.0799 431.6830 3.081 0.00368 **
X2 516.0858 1914.1261 0.270 0.78881
X3 -96.4629 552.2423 -0.175 0.86219
X4 -0.2918 0.4430 -0.659 0.51379
X5 -1337.0089 2698.8468 -0.495 0.62296
X6 2108.2066 2846.9160 0.741 0.46320
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 47250 on 41 degrees of freedom
Multiple R-squared: 0.6872, Adjusted R-squared: 0.6415
F-statistic: 15.02 on 6 and 41 DF, p-value: 5.345e-09
> 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.7840160 0.43196806 0.21598403
[2,] 0.6447050 0.71059010 0.35529505
[3,] 0.6561054 0.68778926 0.34389463
[4,] 0.5915520 0.81689604 0.40844802
[5,] 0.5232831 0.95343381 0.47671690
[6,] 0.4715592 0.94311834 0.52844083
[7,] 0.3652524 0.73050471 0.63474764
[8,] 0.2653010 0.53060194 0.73469903
[9,] 0.1916145 0.38322893 0.80838553
[10,] 0.1890526 0.37810515 0.81094743
[11,] 0.3817639 0.76352781 0.61823610
[12,] 0.3088399 0.61767982 0.69116009
[13,] 0.2698570 0.53971403 0.73014298
[14,] 0.8511159 0.29776818 0.14888409
[15,] 0.7893164 0.42136713 0.21068356
[16,] 0.8205176 0.35896488 0.17948244
[17,] 0.9715054 0.05698924 0.02849462
[18,] 0.9497376 0.10052479 0.05026239
[19,] 0.9203586 0.15928290 0.07964145
[20,] 0.9167981 0.16640374 0.08320187
[21,] 0.9027099 0.19458013 0.09729007
[22,] 0.9134385 0.17312308 0.08656154
[23,] 0.8569308 0.28613833 0.14306917
[24,] 0.8251572 0.34968558 0.17484279
[25,] 0.7959039 0.40819220 0.20409610
[26,] 0.8900498 0.21990046 0.10995023
[27,] 0.8294209 0.34115821 0.17057910
[28,] 0.7295439 0.54091223 0.27045612
[29,] 0.5871644 0.82567130 0.41283565
> postscript(file="/var/wessaorg/rcomp/tmp/17kj31324641930.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/21mtf1324641930.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/3vgtu1324641930.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/4tnop1324641930.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/5tuqb1324641930.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
-17709.3310 -47882.9688 -1988.3147 -74862.5903 -18893.8162 1371.8177
7 8 9 10 11 12
115488.2018 -9257.9520 19710.3239 19810.3307 -12942.9878 14376.6069
13 14 15 16 17 18
-20104.4655 1669.3223 -2419.9243 -301.6766 5156.2688 -21115.2021
19 20 21 22 23 24
47923.1068 -80525.7310 -45085.6080 -29938.2667 126173.8212 15057.9603
25 26 27 28 29 30
-69318.5263 -97834.1194 -7564.2104 18022.3124 -13614.8337 -5303.0660
31 32 33 34 35 36
-20749.3350 -7390.7835 -133.9502 -2158.1349 -10251.8954 32369.9446
37 38 39 40 41 42
8809.5915 43765.3522 127739.0600 -29138.5090 -1848.9947 -419.9344
43 44 45 46 47 48
-30676.8880 36146.2352 24902.8706 -12911.8183 28612.7359 5237.9717
> postscript(file="/var/wessaorg/rcomp/tmp/6405m1324641930.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 -17709.3310 NA
1 -47882.9688 -17709.3310
2 -1988.3147 -47882.9688
3 -74862.5903 -1988.3147
4 -18893.8162 -74862.5903
5 1371.8177 -18893.8162
6 115488.2018 1371.8177
7 -9257.9520 115488.2018
8 19710.3239 -9257.9520
9 19810.3307 19710.3239
10 -12942.9878 19810.3307
11 14376.6069 -12942.9878
12 -20104.4655 14376.6069
13 1669.3223 -20104.4655
14 -2419.9243 1669.3223
15 -301.6766 -2419.9243
16 5156.2688 -301.6766
17 -21115.2021 5156.2688
18 47923.1068 -21115.2021
19 -80525.7310 47923.1068
20 -45085.6080 -80525.7310
21 -29938.2667 -45085.6080
22 126173.8212 -29938.2667
23 15057.9603 126173.8212
24 -69318.5263 15057.9603
25 -97834.1194 -69318.5263
26 -7564.2104 -97834.1194
27 18022.3124 -7564.2104
28 -13614.8337 18022.3124
29 -5303.0660 -13614.8337
30 -20749.3350 -5303.0660
31 -7390.7835 -20749.3350
32 -133.9502 -7390.7835
33 -2158.1349 -133.9502
34 -10251.8954 -2158.1349
35 32369.9446 -10251.8954
36 8809.5915 32369.9446
37 43765.3522 8809.5915
38 127739.0600 43765.3522
39 -29138.5090 127739.0600
40 -1848.9947 -29138.5090
41 -419.9344 -1848.9947
42 -30676.8880 -419.9344
43 36146.2352 -30676.8880
44 24902.8706 36146.2352
45 -12911.8183 24902.8706
46 28612.7359 -12911.8183
47 5237.9717 28612.7359
48 NA 5237.9717
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -47882.9688 -17709.3310
[2,] -1988.3147 -47882.9688
[3,] -74862.5903 -1988.3147
[4,] -18893.8162 -74862.5903
[5,] 1371.8177 -18893.8162
[6,] 115488.2018 1371.8177
[7,] -9257.9520 115488.2018
[8,] 19710.3239 -9257.9520
[9,] 19810.3307 19710.3239
[10,] -12942.9878 19810.3307
[11,] 14376.6069 -12942.9878
[12,] -20104.4655 14376.6069
[13,] 1669.3223 -20104.4655
[14,] -2419.9243 1669.3223
[15,] -301.6766 -2419.9243
[16,] 5156.2688 -301.6766
[17,] -21115.2021 5156.2688
[18,] 47923.1068 -21115.2021
[19,] -80525.7310 47923.1068
[20,] -45085.6080 -80525.7310
[21,] -29938.2667 -45085.6080
[22,] 126173.8212 -29938.2667
[23,] 15057.9603 126173.8212
[24,] -69318.5263 15057.9603
[25,] -97834.1194 -69318.5263
[26,] -7564.2104 -97834.1194
[27,] 18022.3124 -7564.2104
[28,] -13614.8337 18022.3124
[29,] -5303.0660 -13614.8337
[30,] -20749.3350 -5303.0660
[31,] -7390.7835 -20749.3350
[32,] -133.9502 -7390.7835
[33,] -2158.1349 -133.9502
[34,] -10251.8954 -2158.1349
[35,] 32369.9446 -10251.8954
[36,] 8809.5915 32369.9446
[37,] 43765.3522 8809.5915
[38,] 127739.0600 43765.3522
[39,] -29138.5090 127739.0600
[40,] -1848.9947 -29138.5090
[41,] -419.9344 -1848.9947
[42,] -30676.8880 -419.9344
[43,] 36146.2352 -30676.8880
[44,] 24902.8706 36146.2352
[45,] -12911.8183 24902.8706
[46,] 28612.7359 -12911.8183
[47,] 5237.9717 28612.7359
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -47882.9688 -17709.3310
2 -1988.3147 -47882.9688
3 -74862.5903 -1988.3147
4 -18893.8162 -74862.5903
5 1371.8177 -18893.8162
6 115488.2018 1371.8177
7 -9257.9520 115488.2018
8 19710.3239 -9257.9520
9 19810.3307 19710.3239
10 -12942.9878 19810.3307
11 14376.6069 -12942.9878
12 -20104.4655 14376.6069
13 1669.3223 -20104.4655
14 -2419.9243 1669.3223
15 -301.6766 -2419.9243
16 5156.2688 -301.6766
17 -21115.2021 5156.2688
18 47923.1068 -21115.2021
19 -80525.7310 47923.1068
20 -45085.6080 -80525.7310
21 -29938.2667 -45085.6080
22 126173.8212 -29938.2667
23 15057.9603 126173.8212
24 -69318.5263 15057.9603
25 -97834.1194 -69318.5263
26 -7564.2104 -97834.1194
27 18022.3124 -7564.2104
28 -13614.8337 18022.3124
29 -5303.0660 -13614.8337
30 -20749.3350 -5303.0660
31 -7390.7835 -20749.3350
32 -133.9502 -7390.7835
33 -2158.1349 -133.9502
34 -10251.8954 -2158.1349
35 32369.9446 -10251.8954
36 8809.5915 32369.9446
37 43765.3522 8809.5915
38 127739.0600 43765.3522
39 -29138.5090 127739.0600
40 -1848.9947 -29138.5090
41 -419.9344 -1848.9947
42 -30676.8880 -419.9344
43 36146.2352 -30676.8880
44 24902.8706 36146.2352
45 -12911.8183 24902.8706
46 28612.7359 -12911.8183
47 5237.9717 28612.7359
> 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/73blc1324641930.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/8ts311324641930.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/9dn4e1324641930.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/109yu51324641930.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/11enz81324641930.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/12anx31324641930.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/13hnfl1324641930.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/14nl8r1324641930.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/15wt4q1324641930.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/16m7rz1324641930.tab")
+ }
>
> try(system("convert tmp/17kj31324641930.ps tmp/17kj31324641930.png",intern=TRUE))
character(0)
> try(system("convert tmp/21mtf1324641930.ps tmp/21mtf1324641930.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vgtu1324641930.ps tmp/3vgtu1324641930.png",intern=TRUE))
character(0)
> try(system("convert tmp/4tnop1324641930.ps tmp/4tnop1324641930.png",intern=TRUE))
character(0)
> try(system("convert tmp/5tuqb1324641930.ps tmp/5tuqb1324641930.png",intern=TRUE))
character(0)
> try(system("convert tmp/6405m1324641930.ps tmp/6405m1324641930.png",intern=TRUE))
character(0)
> try(system("convert tmp/73blc1324641930.ps tmp/73blc1324641930.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ts311324641930.ps tmp/8ts311324641930.png",intern=TRUE))
character(0)
> try(system("convert tmp/9dn4e1324641930.ps tmp/9dn4e1324641930.png",intern=TRUE))
character(0)
> try(system("convert tmp/109yu51324641930.ps tmp/109yu51324641930.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.022 0.582 3.629