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(8
+ ,78
+ ,284
+ ,9.100000381
+ ,109
+ ,9.300000191
+ ,68
+ ,433
+ ,8.699999809
+ ,144
+ ,7.5
+ ,70
+ ,739
+ ,7.199999809
+ ,113
+ ,8.899999619
+ ,96
+ ,1792
+ ,8.899999619
+ ,97
+ ,10.19999981
+ ,74
+ ,477
+ ,8.300000191
+ ,206
+ ,8.300000191
+ ,111
+ ,362
+ ,10.89999962
+ ,124
+ ,8.800000191
+ ,77
+ ,671
+ ,10
+ ,152
+ ,8.800000191
+ ,168
+ ,636
+ ,9.100000381
+ ,162
+ ,10.69999981
+ ,82
+ ,329
+ ,8.699999809
+ ,150
+ ,11.69999981
+ ,89
+ ,634
+ ,7.599999905
+ ,134
+ ,8.5
+ ,149
+ ,631
+ ,10.80000019
+ ,292
+ ,8.300000191
+ ,60
+ ,257
+ ,9.5
+ ,108
+ ,8.199999809
+ ,96
+ ,284
+ ,8.800000191
+ ,111
+ ,7.900000095
+ ,83
+ ,603
+ ,9.5
+ ,182
+ ,10.30000019
+ ,130
+ ,686
+ ,8.699999809
+ ,129
+ ,7.400000095
+ ,145
+ ,345
+ ,11.19999981
+ ,158
+ ,9.600000381
+ ,112
+ ,1357
+ ,9.699999809
+ ,186
+ ,9.300000191
+ ,131
+ ,544
+ ,9.600000381
+ ,177
+ ,10.60000038
+ ,80
+ ,205
+ ,9.100000381
+ ,127
+ ,9.699999809
+ ,130
+ ,1264
+ ,9.199999809
+ ,179
+ ,11.60000038
+ ,140
+ ,688
+ ,8.300000191
+ ,80
+ ,8.100000381
+ ,154
+ ,354
+ ,8.399999619
+ ,103
+ ,9.800000191
+ ,118
+ ,1632
+ ,9.399999619
+ ,101
+ ,7.400000095
+ ,94
+ ,348
+ ,9.800000191
+ ,117
+ ,9.399999619
+ ,119
+ ,370
+ ,10.39999962
+ ,88
+ ,11.19999981
+ ,153
+ ,648
+ ,9.899999619
+ ,78
+ ,9.100000381
+ ,116
+ ,366
+ ,9.199999809
+ ,102
+ ,10.5
+ ,97
+ ,540
+ ,10.30000019
+ ,95
+ ,11.89999962
+ ,176
+ ,680
+ ,8.899999619
+ ,80
+ ,8.399999619
+ ,75
+ ,345
+ ,9.600000381
+ ,92
+ ,5
+ ,134
+ ,525
+ ,10.30000019
+ ,126
+ ,9.800000191
+ ,161
+ ,870
+ ,10.39999962
+ ,108
+ ,9.800000191
+ ,111
+ ,669
+ ,9.699999809
+ ,77
+ ,10.80000019
+ ,114
+ ,452
+ ,9.600000381
+ ,60
+ ,10.10000038
+ ,142
+ ,430
+ ,10.69999981
+ ,71
+ ,10.89999962
+ ,238
+ ,822
+ ,10.30000019
+ ,86
+ ,9.199999809
+ ,78
+ ,190
+ ,10.69999981
+ ,93
+ ,8.300000191
+ ,196
+ ,867
+ ,9.600000381
+ ,106
+ ,7.300000191
+ ,125
+ ,969
+ ,10.5
+ ,162
+ ,9.399999619
+ ,82
+ ,499
+ ,7.699999809
+ ,95
+ ,9.399999619
+ ,125
+ ,925
+ ,10.19999981
+ ,91
+ ,9.800000191
+ ,129
+ ,353
+ ,9.899999619
+ ,52
+ ,3.599999905
+ ,84
+ ,288
+ ,8.399999619
+ ,110
+ ,8.399999619
+ ,183
+ ,718
+ ,10.39999962
+ ,69
+ ,10.80000019
+ ,119
+ ,540
+ ,9.199999809
+ ,57
+ ,10.10000038
+ ,180
+ ,668
+ ,13
+ ,106
+ ,9
+ ,82
+ ,347
+ ,8.800000191
+ ,40
+ ,10
+ ,71
+ ,345
+ ,9.199999809
+ ,50
+ ,11.30000019
+ ,118
+ ,463
+ ,7.800000191
+ ,35
+ ,11.30000019
+ ,121
+ ,728
+ ,8.199999809
+ ,86
+ ,12.80000019
+ ,68
+ ,383
+ ,7.400000095
+ ,57
+ ,10
+ ,112
+ ,316
+ ,10.39999962
+ ,57
+ ,6.699999809
+ ,109
+ ,388
+ ,8.899999619
+ ,94)
+ ,dim=c(5
+ ,53)
+ ,dimnames=list(c('death'
+ ,'doctor'
+ ,'hospital'
+ ,'income'
+ ,'population')
+ ,1:53))
> y <- array(NA,dim=c(5,53),dimnames=list(c('death','doctor','hospital','income','population'),1:53))
> 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
doctor death hospital income population
1 78 8.0 284 9.1 109
2 68 9.3 433 8.7 144
3 70 7.5 739 7.2 113
4 96 8.9 1792 8.9 97
5 74 10.2 477 8.3 206
6 111 8.3 362 10.9 124
7 77 8.8 671 10.0 152
8 168 8.8 636 9.1 162
9 82 10.7 329 8.7 150
10 89 11.7 634 7.6 134
11 149 8.5 631 10.8 292
12 60 8.3 257 9.5 108
13 96 8.2 284 8.8 111
14 83 7.9 603 9.5 182
15 130 10.3 686 8.7 129
16 145 7.4 345 11.2 158
17 112 9.6 1357 9.7 186
18 131 9.3 544 9.6 177
19 80 10.6 205 9.1 127
20 130 9.7 1264 9.2 179
21 140 11.6 688 8.3 80
22 154 8.1 354 8.4 103
23 118 9.8 1632 9.4 101
24 94 7.4 348 9.8 117
25 119 9.4 370 10.4 88
26 153 11.2 648 9.9 78
27 116 9.1 366 9.2 102
28 97 10.5 540 10.3 95
29 176 11.9 680 8.9 80
30 75 8.4 345 9.6 92
31 134 5.0 525 10.3 126
32 161 9.8 870 10.4 108
33 111 9.8 669 9.7 77
34 114 10.8 452 9.6 60
35 142 10.1 430 10.7 71
36 238 10.9 822 10.3 86
37 78 9.2 190 10.7 93
38 196 8.3 867 9.6 106
39 125 7.3 969 10.5 162
40 82 9.4 499 7.7 95
41 125 9.4 925 10.2 91
42 129 9.8 353 9.9 52
43 84 3.6 288 8.4 110
44 183 8.4 718 10.4 69
45 119 10.8 540 9.2 57
46 180 10.1 668 13.0 106
47 82 9.0 347 8.8 40
48 71 10.0 345 9.2 50
49 118 11.3 463 7.8 35
50 121 11.3 728 8.2 86
51 68 12.8 383 7.4 57
52 112 10.0 316 10.4 57
53 109 6.7 388 8.9 94
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) death hospital income population
-77.11802 3.12905 0.03251 16.24810 -0.07585
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-50.233 -20.272 -4.423 18.384 93.459
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -77.11802 54.13533 -1.425 0.160759
death 3.12905 2.93518 1.066 0.291734
hospital 0.03251 0.01420 2.289 0.026498 *
income 16.24810 4.33021 3.752 0.000473 ***
population -0.07585 0.10382 -0.731 0.468587
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 32.95 on 48 degrees of freedom
Multiple R-squared: 0.302, Adjusted R-squared: 0.2438
F-statistic: 5.192 on 4 and 48 DF, p-value: 0.001476
> 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.8209078 0.35818433 0.17909216
[2,] 0.7745060 0.45098805 0.22549402
[3,] 0.7062329 0.58753429 0.29376714
[4,] 0.5986636 0.80267277 0.40133638
[5,] 0.5499019 0.90019611 0.45009806
[6,] 0.4661375 0.93227510 0.53386245
[7,] 0.4178425 0.83568492 0.58215754
[8,] 0.4406627 0.88132539 0.55933730
[9,] 0.4147102 0.82942040 0.58528980
[10,] 0.3640316 0.72806327 0.63596836
[11,] 0.3126256 0.62525111 0.68737444
[12,] 0.2437675 0.48753502 0.75623249
[13,] 0.1915897 0.38317946 0.80841027
[14,] 0.2451794 0.49035874 0.75482063
[15,] 0.5178764 0.96424727 0.48212363
[16,] 0.7540984 0.49180324 0.24590162
[17,] 0.6831318 0.63373645 0.31686823
[18,] 0.6070414 0.78591721 0.39295860
[19,] 0.5696370 0.86072595 0.43036297
[20,] 0.5388788 0.92224233 0.46112117
[21,] 0.5218981 0.95620388 0.47810194
[22,] 0.6835345 0.63293103 0.31646551
[23,] 0.6450741 0.70985178 0.35492589
[24,] 0.6200829 0.75983413 0.37991707
[25,] 0.5589687 0.88206266 0.44103133
[26,] 0.5390693 0.92186141 0.46093071
[27,] 0.4496039 0.89920783 0.55039609
[28,] 0.3719050 0.74381009 0.62809495
[29,] 0.8632971 0.27340589 0.13670295
[30,] 0.8332011 0.33359775 0.16679888
[31,] 0.9489201 0.10215981 0.05107991
[32,] 0.9291842 0.14163156 0.07081578
[33,] 0.8781646 0.24367083 0.12183542
[34,] 0.9713847 0.05723056 0.02861528
[35,] 0.9671754 0.06564919 0.03282460
[36,] 0.9263043 0.14739149 0.07369575
[37,] 0.8697088 0.26058231 0.13029115
[38,] 0.7373993 0.52520132 0.26260066
> postscript(file="/var/wessaorg/rcomp/tmp/1plgi1322153156.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/2lwf71322153156.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/336i91322153156.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/4a9af1322153156.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/5g1sh1322153156.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 = 53
Frequency = 1
1 2 3 4 5 6
-18.7362045 -28.4934096 -8.7873703 -50.2328645 -15.5378388 -17.3192409
7 8 9 10 11 12
-46.1811840 61.3383442 -15.0382816 -4.4225433 25.6784580 -43.3723273
13 14 15 16 17 18
3.6641230 -24.7550099 21.0156103 17.7540155 -28.5305729 18.7781417
19 20 21 22 23 24
-20.9383867 -0.2272686 29.6653689 65.5940066 -33.6686009 -13.7060595
25 26 27 28 29 30
-7.6278298 19.0685946 11.0005441 -36.4401440 55.2378573 -34.3842439
31 32 33 34 35 36
20.6085937 18.3842756 -16.0596134 -8.7993967 5.0675376 93.4587202
37 38 39 40 41 42
-46.6460134 71.0221334 -10.5401195 -4.4203465 -16.1918171 7.0665751
43 44 45 46 47 48
12.3511098 46.7477527 -0.3882793 0.6151154 -20.2724559 -40.0772137
49 50 51 52 53
20.6288200 12.3837427 -23.2962560 -17.1012861 15.0627381
> postscript(file="/var/wessaorg/rcomp/tmp/6ho8e1322153156.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 = 53
Frequency = 1
lag(myerror, k = 1) myerror
0 -18.7362045 NA
1 -28.4934096 -18.7362045
2 -8.7873703 -28.4934096
3 -50.2328645 -8.7873703
4 -15.5378388 -50.2328645
5 -17.3192409 -15.5378388
6 -46.1811840 -17.3192409
7 61.3383442 -46.1811840
8 -15.0382816 61.3383442
9 -4.4225433 -15.0382816
10 25.6784580 -4.4225433
11 -43.3723273 25.6784580
12 3.6641230 -43.3723273
13 -24.7550099 3.6641230
14 21.0156103 -24.7550099
15 17.7540155 21.0156103
16 -28.5305729 17.7540155
17 18.7781417 -28.5305729
18 -20.9383867 18.7781417
19 -0.2272686 -20.9383867
20 29.6653689 -0.2272686
21 65.5940066 29.6653689
22 -33.6686009 65.5940066
23 -13.7060595 -33.6686009
24 -7.6278298 -13.7060595
25 19.0685946 -7.6278298
26 11.0005441 19.0685946
27 -36.4401440 11.0005441
28 55.2378573 -36.4401440
29 -34.3842439 55.2378573
30 20.6085937 -34.3842439
31 18.3842756 20.6085937
32 -16.0596134 18.3842756
33 -8.7993967 -16.0596134
34 5.0675376 -8.7993967
35 93.4587202 5.0675376
36 -46.6460134 93.4587202
37 71.0221334 -46.6460134
38 -10.5401195 71.0221334
39 -4.4203465 -10.5401195
40 -16.1918171 -4.4203465
41 7.0665751 -16.1918171
42 12.3511098 7.0665751
43 46.7477527 12.3511098
44 -0.3882793 46.7477527
45 0.6151154 -0.3882793
46 -20.2724559 0.6151154
47 -40.0772137 -20.2724559
48 20.6288200 -40.0772137
49 12.3837427 20.6288200
50 -23.2962560 12.3837427
51 -17.1012861 -23.2962560
52 15.0627381 -17.1012861
53 NA 15.0627381
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -28.4934096 -18.7362045
[2,] -8.7873703 -28.4934096
[3,] -50.2328645 -8.7873703
[4,] -15.5378388 -50.2328645
[5,] -17.3192409 -15.5378388
[6,] -46.1811840 -17.3192409
[7,] 61.3383442 -46.1811840
[8,] -15.0382816 61.3383442
[9,] -4.4225433 -15.0382816
[10,] 25.6784580 -4.4225433
[11,] -43.3723273 25.6784580
[12,] 3.6641230 -43.3723273
[13,] -24.7550099 3.6641230
[14,] 21.0156103 -24.7550099
[15,] 17.7540155 21.0156103
[16,] -28.5305729 17.7540155
[17,] 18.7781417 -28.5305729
[18,] -20.9383867 18.7781417
[19,] -0.2272686 -20.9383867
[20,] 29.6653689 -0.2272686
[21,] 65.5940066 29.6653689
[22,] -33.6686009 65.5940066
[23,] -13.7060595 -33.6686009
[24,] -7.6278298 -13.7060595
[25,] 19.0685946 -7.6278298
[26,] 11.0005441 19.0685946
[27,] -36.4401440 11.0005441
[28,] 55.2378573 -36.4401440
[29,] -34.3842439 55.2378573
[30,] 20.6085937 -34.3842439
[31,] 18.3842756 20.6085937
[32,] -16.0596134 18.3842756
[33,] -8.7993967 -16.0596134
[34,] 5.0675376 -8.7993967
[35,] 93.4587202 5.0675376
[36,] -46.6460134 93.4587202
[37,] 71.0221334 -46.6460134
[38,] -10.5401195 71.0221334
[39,] -4.4203465 -10.5401195
[40,] -16.1918171 -4.4203465
[41,] 7.0665751 -16.1918171
[42,] 12.3511098 7.0665751
[43,] 46.7477527 12.3511098
[44,] -0.3882793 46.7477527
[45,] 0.6151154 -0.3882793
[46,] -20.2724559 0.6151154
[47,] -40.0772137 -20.2724559
[48,] 20.6288200 -40.0772137
[49,] 12.3837427 20.6288200
[50,] -23.2962560 12.3837427
[51,] -17.1012861 -23.2962560
[52,] 15.0627381 -17.1012861
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -28.4934096 -18.7362045
2 -8.7873703 -28.4934096
3 -50.2328645 -8.7873703
4 -15.5378388 -50.2328645
5 -17.3192409 -15.5378388
6 -46.1811840 -17.3192409
7 61.3383442 -46.1811840
8 -15.0382816 61.3383442
9 -4.4225433 -15.0382816
10 25.6784580 -4.4225433
11 -43.3723273 25.6784580
12 3.6641230 -43.3723273
13 -24.7550099 3.6641230
14 21.0156103 -24.7550099
15 17.7540155 21.0156103
16 -28.5305729 17.7540155
17 18.7781417 -28.5305729
18 -20.9383867 18.7781417
19 -0.2272686 -20.9383867
20 29.6653689 -0.2272686
21 65.5940066 29.6653689
22 -33.6686009 65.5940066
23 -13.7060595 -33.6686009
24 -7.6278298 -13.7060595
25 19.0685946 -7.6278298
26 11.0005441 19.0685946
27 -36.4401440 11.0005441
28 55.2378573 -36.4401440
29 -34.3842439 55.2378573
30 20.6085937 -34.3842439
31 18.3842756 20.6085937
32 -16.0596134 18.3842756
33 -8.7993967 -16.0596134
34 5.0675376 -8.7993967
35 93.4587202 5.0675376
36 -46.6460134 93.4587202
37 71.0221334 -46.6460134
38 -10.5401195 71.0221334
39 -4.4203465 -10.5401195
40 -16.1918171 -4.4203465
41 7.0665751 -16.1918171
42 12.3511098 7.0665751
43 46.7477527 12.3511098
44 -0.3882793 46.7477527
45 0.6151154 -0.3882793
46 -20.2724559 0.6151154
47 -40.0772137 -20.2724559
48 20.6288200 -40.0772137
49 12.3837427 20.6288200
50 -23.2962560 12.3837427
51 -17.1012861 -23.2962560
52 15.0627381 -17.1012861
> 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/7qduq1322153156.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/8pxhf1322153156.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/9iyx41322153156.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/10ucbx1322153156.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/11jt7u1322153156.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/123iaj1322153156.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/13j2101322153157.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/14q9fg1322153157.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/15n97d1322153157.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/167osp1322153157.tab")
+ }
>
> try(system("convert tmp/1plgi1322153156.ps tmp/1plgi1322153156.png",intern=TRUE))
character(0)
> try(system("convert tmp/2lwf71322153156.ps tmp/2lwf71322153156.png",intern=TRUE))
character(0)
> try(system("convert tmp/336i91322153156.ps tmp/336i91322153156.png",intern=TRUE))
character(0)
> try(system("convert tmp/4a9af1322153156.ps tmp/4a9af1322153156.png",intern=TRUE))
character(0)
> try(system("convert tmp/5g1sh1322153156.ps tmp/5g1sh1322153156.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ho8e1322153156.ps tmp/6ho8e1322153156.png",intern=TRUE))
character(0)
> try(system("convert tmp/7qduq1322153156.ps tmp/7qduq1322153156.png",intern=TRUE))
character(0)
> try(system("convert tmp/8pxhf1322153156.ps tmp/8pxhf1322153156.png",intern=TRUE))
character(0)
> try(system("convert tmp/9iyx41322153156.ps tmp/9iyx41322153156.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ucbx1322153156.ps tmp/10ucbx1322153156.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.068 0.497 3.648