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('deatch'
+ ,'doctor'
+ ,'hospital'
+ ,'income'
+ ,'population')
+ ,1:53))
> y <- array(NA,dim=c(5,53),dimnames=list(c('deatch','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 = '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
deatch doctor hospital income population t
1 8.0 78 284 9.1 109 1
2 9.3 68 433 8.7 144 2
3 7.5 70 739 7.2 113 3
4 8.9 96 1792 8.9 97 4
5 10.2 74 477 8.3 206 5
6 8.3 111 362 10.9 124 6
7 8.8 77 671 10.0 152 7
8 8.8 168 636 9.1 162 8
9 10.7 82 329 8.7 150 9
10 11.7 89 634 7.6 134 10
11 8.5 149 631 10.8 292 11
12 8.3 60 257 9.5 108 12
13 8.2 96 284 8.8 111 13
14 7.9 83 603 9.5 182 14
15 10.3 130 686 8.7 129 15
16 7.4 145 345 11.2 158 16
17 9.6 112 1357 9.7 186 17
18 9.3 131 544 9.6 177 18
19 10.6 80 205 9.1 127 19
20 9.7 130 1264 9.2 179 20
21 11.6 140 688 8.3 80 21
22 8.1 154 354 8.4 103 22
23 9.8 118 1632 9.4 101 23
24 7.4 94 348 9.8 117 24
25 9.4 119 370 10.4 88 25
26 11.2 153 648 9.9 78 26
27 9.1 116 366 9.2 102 27
28 10.5 97 540 10.3 95 28
29 11.9 176 680 8.9 80 29
30 8.4 75 345 9.6 92 30
31 5.0 134 525 10.3 126 31
32 9.8 161 870 10.4 108 32
33 9.8 111 669 9.7 77 33
34 10.8 114 452 9.6 60 34
35 10.1 142 430 10.7 71 35
36 10.9 238 822 10.3 86 36
37 9.2 78 190 10.7 93 37
38 8.3 196 867 9.6 106 38
39 7.3 125 969 10.5 162 39
40 9.4 82 499 7.7 95 40
41 9.4 125 925 10.2 91 41
42 9.8 129 353 9.9 52 42
43 3.6 84 288 8.4 110 43
44 8.4 183 718 10.4 69 44
45 10.8 119 540 9.2 57 45
46 10.1 180 668 13.0 106 46
47 9.0 82 347 8.8 40 47
48 10.0 71 345 9.2 50 48
49 11.3 118 463 7.8 35 49
50 11.3 121 728 8.2 86 50
51 12.8 68 383 7.4 57 51
52 10.0 112 316 10.4 57 52
53 6.7 109 388 8.9 94 53
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) doctor hospital income population t
12.5335133 0.0081609 0.0005474 -0.3122527 -0.0115497 -0.0101436
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.4471 -0.7061 0.2474 0.9873 2.9882
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.5335133 2.1016372 5.964 3.05e-07 ***
doctor 0.0081609 0.0071471 1.142 0.2593
hospital 0.0005474 0.0007310 0.749 0.4577
income -0.3122527 0.2389660 -1.307 0.1977
population -0.0115497 0.0063913 -1.807 0.0772 .
t -0.0101436 0.0198016 -0.512 0.6109
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.614 on 47 degrees of freedom
Multiple R-squared: 0.1485, Adjusted R-squared: 0.05787
F-statistic: 1.639 on 5 and 47 DF, p-value: 0.1683
> 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.2495768 0.4991537 0.7504232
[2,] 0.2251862 0.4503723 0.7748138
[3,] 0.2296805 0.4593609 0.7703195
[4,] 0.3595046 0.7190092 0.6404954
[5,] 0.3427543 0.6855086 0.6572457
[6,] 0.3287386 0.6574771 0.6712614
[7,] 0.2613618 0.5227235 0.7386382
[8,] 0.1945983 0.3891965 0.8054017
[9,] 0.1353861 0.2707722 0.8646139
[10,] 0.1070959 0.2141919 0.8929041
[11,] 0.1280550 0.2561100 0.8719450
[12,] 0.1216007 0.2432015 0.8783993
[13,] 0.1142644 0.2285288 0.8857356
[14,] 0.1670266 0.3340532 0.8329734
[15,] 0.1758467 0.3516933 0.8241533
[16,] 0.1878720 0.3757440 0.8121280
[17,] 0.1465928 0.2931857 0.8534072
[18,] 0.1588542 0.3177083 0.8411458
[19,] 0.1301796 0.2603592 0.8698204
[20,] 0.1187409 0.2374818 0.8812591
[21,] 0.1442829 0.2885658 0.8557171
[22,] 0.1194985 0.2389969 0.8805015
[23,] 0.3683245 0.7366491 0.6316755
[24,] 0.2908963 0.5817925 0.7091037
[25,] 0.2261759 0.4523519 0.7738241
[26,] 0.1781446 0.3562892 0.8218554
[27,] 0.1334298 0.2668596 0.8665702
[28,] 0.1325520 0.2651041 0.8674480
[29,] 0.1801910 0.3603820 0.8198090
[30,] 0.1567007 0.3134015 0.8432993
[31,] 0.1196360 0.2392719 0.8803640
[32,] 0.1858446 0.3716892 0.8141554
[33,] 0.1689314 0.3378628 0.8310686
[34,] 0.2577953 0.5155907 0.7422047
[35,] 0.3307469 0.6614937 0.6692531
[36,] 0.2728042 0.5456083 0.7271958
> postscript(file="/var/wessaorg/rcomp/tmp/1ynkm1322154978.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/2bu6r1322154978.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/3jc4d1322154978.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/4rvzm1322154978.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/5r5yg1322154978.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
-1.214956249 0.374578012 -2.425513912 -1.457895306 1.823139478 -0.440943258
7 8 9 10 11 12
0.219900943 -0.658973006 1.857553575 2.115351244 1.261543385 -0.528350019
13 14 15 16 17 18
-1.110707025 -0.430474265 0.688737607 -1.021308463 0.759226275 0.624146583
19 20 21 22 23 24
1.802442780 0.556690138 1.276058510 -1.848363017 -0.254804229 -1.436286099
25 26 27 28 29 30
0.210201598 1.319082913 -0.255845584 1.476744124 1.555142953 -0.569918528
31 32 33 34 35 36
-3.928528409 0.295759843 0.247352581 1.124221282 0.688425469 0.548896446
37 38 39 40 41 42
0.716472737 -1.700270756 -1.238721178 -0.168534810 -0.008055611 0.138422722
43 44 45 46 47 48
-5.447109430 -1.529297757 0.987277305 1.482037440 -0.706085885 0.635320702
49 50 51 52 53
0.886912298 1.441457503 2.988227665 0.812721347 -2.533102664
> postscript(file="/var/wessaorg/rcomp/tmp/6nze21322154978.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 -1.214956249 NA
1 0.374578012 -1.214956249
2 -2.425513912 0.374578012
3 -1.457895306 -2.425513912
4 1.823139478 -1.457895306
5 -0.440943258 1.823139478
6 0.219900943 -0.440943258
7 -0.658973006 0.219900943
8 1.857553575 -0.658973006
9 2.115351244 1.857553575
10 1.261543385 2.115351244
11 -0.528350019 1.261543385
12 -1.110707025 -0.528350019
13 -0.430474265 -1.110707025
14 0.688737607 -0.430474265
15 -1.021308463 0.688737607
16 0.759226275 -1.021308463
17 0.624146583 0.759226275
18 1.802442780 0.624146583
19 0.556690138 1.802442780
20 1.276058510 0.556690138
21 -1.848363017 1.276058510
22 -0.254804229 -1.848363017
23 -1.436286099 -0.254804229
24 0.210201598 -1.436286099
25 1.319082913 0.210201598
26 -0.255845584 1.319082913
27 1.476744124 -0.255845584
28 1.555142953 1.476744124
29 -0.569918528 1.555142953
30 -3.928528409 -0.569918528
31 0.295759843 -3.928528409
32 0.247352581 0.295759843
33 1.124221282 0.247352581
34 0.688425469 1.124221282
35 0.548896446 0.688425469
36 0.716472737 0.548896446
37 -1.700270756 0.716472737
38 -1.238721178 -1.700270756
39 -0.168534810 -1.238721178
40 -0.008055611 -0.168534810
41 0.138422722 -0.008055611
42 -5.447109430 0.138422722
43 -1.529297757 -5.447109430
44 0.987277305 -1.529297757
45 1.482037440 0.987277305
46 -0.706085885 1.482037440
47 0.635320702 -0.706085885
48 0.886912298 0.635320702
49 1.441457503 0.886912298
50 2.988227665 1.441457503
51 0.812721347 2.988227665
52 -2.533102664 0.812721347
53 NA -2.533102664
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.374578012 -1.214956249
[2,] -2.425513912 0.374578012
[3,] -1.457895306 -2.425513912
[4,] 1.823139478 -1.457895306
[5,] -0.440943258 1.823139478
[6,] 0.219900943 -0.440943258
[7,] -0.658973006 0.219900943
[8,] 1.857553575 -0.658973006
[9,] 2.115351244 1.857553575
[10,] 1.261543385 2.115351244
[11,] -0.528350019 1.261543385
[12,] -1.110707025 -0.528350019
[13,] -0.430474265 -1.110707025
[14,] 0.688737607 -0.430474265
[15,] -1.021308463 0.688737607
[16,] 0.759226275 -1.021308463
[17,] 0.624146583 0.759226275
[18,] 1.802442780 0.624146583
[19,] 0.556690138 1.802442780
[20,] 1.276058510 0.556690138
[21,] -1.848363017 1.276058510
[22,] -0.254804229 -1.848363017
[23,] -1.436286099 -0.254804229
[24,] 0.210201598 -1.436286099
[25,] 1.319082913 0.210201598
[26,] -0.255845584 1.319082913
[27,] 1.476744124 -0.255845584
[28,] 1.555142953 1.476744124
[29,] -0.569918528 1.555142953
[30,] -3.928528409 -0.569918528
[31,] 0.295759843 -3.928528409
[32,] 0.247352581 0.295759843
[33,] 1.124221282 0.247352581
[34,] 0.688425469 1.124221282
[35,] 0.548896446 0.688425469
[36,] 0.716472737 0.548896446
[37,] -1.700270756 0.716472737
[38,] -1.238721178 -1.700270756
[39,] -0.168534810 -1.238721178
[40,] -0.008055611 -0.168534810
[41,] 0.138422722 -0.008055611
[42,] -5.447109430 0.138422722
[43,] -1.529297757 -5.447109430
[44,] 0.987277305 -1.529297757
[45,] 1.482037440 0.987277305
[46,] -0.706085885 1.482037440
[47,] 0.635320702 -0.706085885
[48,] 0.886912298 0.635320702
[49,] 1.441457503 0.886912298
[50,] 2.988227665 1.441457503
[51,] 0.812721347 2.988227665
[52,] -2.533102664 0.812721347
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.374578012 -1.214956249
2 -2.425513912 0.374578012
3 -1.457895306 -2.425513912
4 1.823139478 -1.457895306
5 -0.440943258 1.823139478
6 0.219900943 -0.440943258
7 -0.658973006 0.219900943
8 1.857553575 -0.658973006
9 2.115351244 1.857553575
10 1.261543385 2.115351244
11 -0.528350019 1.261543385
12 -1.110707025 -0.528350019
13 -0.430474265 -1.110707025
14 0.688737607 -0.430474265
15 -1.021308463 0.688737607
16 0.759226275 -1.021308463
17 0.624146583 0.759226275
18 1.802442780 0.624146583
19 0.556690138 1.802442780
20 1.276058510 0.556690138
21 -1.848363017 1.276058510
22 -0.254804229 -1.848363017
23 -1.436286099 -0.254804229
24 0.210201598 -1.436286099
25 1.319082913 0.210201598
26 -0.255845584 1.319082913
27 1.476744124 -0.255845584
28 1.555142953 1.476744124
29 -0.569918528 1.555142953
30 -3.928528409 -0.569918528
31 0.295759843 -3.928528409
32 0.247352581 0.295759843
33 1.124221282 0.247352581
34 0.688425469 1.124221282
35 0.548896446 0.688425469
36 0.716472737 0.548896446
37 -1.700270756 0.716472737
38 -1.238721178 -1.700270756
39 -0.168534810 -1.238721178
40 -0.008055611 -0.168534810
41 0.138422722 -0.008055611
42 -5.447109430 0.138422722
43 -1.529297757 -5.447109430
44 0.987277305 -1.529297757
45 1.482037440 0.987277305
46 -0.706085885 1.482037440
47 0.635320702 -0.706085885
48 0.886912298 0.635320702
49 1.441457503 0.886912298
50 2.988227665 1.441457503
51 0.812721347 2.988227665
52 -2.533102664 0.812721347
> 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/7bbbd1322154978.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/8kspi1322154978.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/9g95s1322154978.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/104bjo1322154978.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/11no1f1322154978.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/12j01b1322154978.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/1360ru1322154978.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/14z6wb1322154978.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/15i7o21322154978.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/16zdkf1322154978.tab")
+ }
>
> try(system("convert tmp/1ynkm1322154978.ps tmp/1ynkm1322154978.png",intern=TRUE))
character(0)
> try(system("convert tmp/2bu6r1322154978.ps tmp/2bu6r1322154978.png",intern=TRUE))
character(0)
> try(system("convert tmp/3jc4d1322154978.ps tmp/3jc4d1322154978.png",intern=TRUE))
character(0)
> try(system("convert tmp/4rvzm1322154978.ps tmp/4rvzm1322154978.png",intern=TRUE))
character(0)
> try(system("convert tmp/5r5yg1322154978.ps tmp/5r5yg1322154978.png",intern=TRUE))
character(0)
> try(system("convert tmp/6nze21322154978.ps tmp/6nze21322154978.png",intern=TRUE))
character(0)
> try(system("convert tmp/7bbbd1322154978.ps tmp/7bbbd1322154978.png",intern=TRUE))
character(0)
> try(system("convert tmp/8kspi1322154978.ps tmp/8kspi1322154978.png",intern=TRUE))
character(0)
> try(system("convert tmp/9g95s1322154978.ps tmp/9g95s1322154978.png",intern=TRUE))
character(0)
> try(system("convert tmp/104bjo1322154978.ps tmp/104bjo1322154978.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.061 0.483 3.570