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(888
+ ,51
+ ,256
+ ,5
+ ,10345
+ ,545
+ ,24
+ ,160
+ ,7
+ ,17607
+ ,186
+ ,17
+ ,70
+ ,0
+ ,1423
+ ,1405
+ ,66
+ ,360
+ ,12
+ ,20050
+ ,2047
+ ,85
+ ,721
+ ,15
+ ,21212
+ ,3626
+ ,130
+ ,938
+ ,16
+ ,93979
+ ,845
+ ,36
+ ,287
+ ,12
+ ,15524
+ ,643
+ ,32
+ ,149
+ ,13
+ ,16182
+ ,1181
+ ,33
+ ,311
+ ,15
+ ,19238
+ ,1835
+ ,64
+ ,617
+ ,13
+ ,28909
+ ,855
+ ,35
+ ,262
+ ,6
+ ,22357
+ ,1245
+ ,46
+ ,385
+ ,16
+ ,25560
+ ,993
+ ,69
+ ,369
+ ,7
+ ,9954
+ ,1685
+ ,61
+ ,558
+ ,12
+ ,18490
+ ,741
+ ,24
+ ,220
+ ,9
+ ,17777
+ ,854
+ ,38
+ ,313
+ ,10
+ ,25268
+ ,949
+ ,34
+ ,229
+ ,16
+ ,37525
+ ,331
+ ,20
+ ,88
+ ,5
+ ,6023
+ ,1602
+ ,54
+ ,494
+ ,20
+ ,25042
+ ,510
+ ,16
+ ,153
+ ,7
+ ,35713
+ ,628
+ ,37
+ ,234
+ ,13
+ ,7039
+ ,1279
+ ,51
+ ,361
+ ,13
+ ,40841
+ ,767
+ ,28
+ ,280
+ ,11
+ ,9214
+ ,1156
+ ,32
+ ,331
+ ,9
+ ,17446
+ ,1120
+ ,51
+ ,378
+ ,10
+ ,10295
+ ,622
+ ,12
+ ,227
+ ,7
+ ,13206
+ ,1203
+ ,98
+ ,396
+ ,13
+ ,26093
+ ,745
+ ,53
+ ,179
+ ,15
+ ,20744
+ ,1535
+ ,61
+ ,509
+ ,13
+ ,68013
+ ,1234
+ ,25
+ ,504
+ ,7
+ ,12840
+ ,757
+ ,28
+ ,225
+ ,14
+ ,12672
+ ,1087
+ ,23
+ ,366
+ ,11
+ ,10872
+ ,1105
+ ,60
+ ,341
+ ,3
+ ,21325
+ ,592
+ ,40
+ ,171
+ ,8
+ ,24542
+ ,1305
+ ,29
+ ,437
+ ,12
+ ,16401
+ ,0
+ ,0
+ ,0
+ ,0
+ ,0
+ ,705
+ ,33
+ ,313
+ ,12
+ ,12821
+ ,1188
+ ,34
+ ,366
+ ,8
+ ,14662
+ ,1110
+ ,21
+ ,232
+ ,20
+ ,22190
+ ,1094
+ ,35
+ ,389
+ ,18
+ ,37929
+ ,1062
+ ,34
+ ,340
+ ,9
+ ,18009
+ ,748
+ ,26
+ ,316
+ ,14
+ ,11076
+ ,404
+ ,12
+ ,140
+ ,7
+ ,24981
+ ,1076
+ ,45
+ ,419
+ ,13
+ ,30691
+ ,673
+ ,29
+ ,226
+ ,11
+ ,29164
+ ,517
+ ,36
+ ,161
+ ,11
+ ,13985
+ ,354
+ ,13
+ ,103
+ ,14
+ ,7588
+ ,1011
+ ,54
+ ,356
+ ,9
+ ,20023
+ ,890
+ ,39
+ ,293
+ ,12
+ ,25524
+ ,1067
+ ,28
+ ,414
+ ,11
+ ,14717
+ ,518
+ ,21
+ ,156
+ ,17
+ ,6832
+ ,697
+ ,36
+ ,189
+ ,10
+ ,9624
+ ,1095
+ ,44
+ ,442
+ ,11
+ ,24300
+ ,928
+ ,44
+ ,321
+ ,12
+ ,21790
+ ,1008
+ ,33
+ ,367
+ ,17
+ ,16493
+ ,951
+ ,30
+ ,309
+ ,6
+ ,9269
+ ,779
+ ,27
+ ,235
+ ,8
+ ,20105
+ ,439
+ ,12
+ ,137
+ ,12
+ ,11216
+ ,568
+ ,37
+ ,194
+ ,13
+ ,15569
+ ,614
+ ,24
+ ,220
+ ,14
+ ,21799
+ ,499
+ ,21
+ ,149
+ ,17
+ ,3772)
+ ,dim=c(5
+ ,61)
+ ,dimnames=list(c('Y_t'
+ ,'X_1t'
+ ,'X_2t'
+ ,'X_3t'
+ ,'X_4t')
+ ,1:61))
> y <- array(NA,dim=c(5,61),dimnames=list(c('Y_t','X_1t','X_2t','X_3t','X_4t'),1:61))
> 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_t X_1t X_2t X_3t X_4t
1 888 51 256 5 10345
2 545 24 160 7 17607
3 186 17 70 0 1423
4 1405 66 360 12 20050
5 2047 85 721 15 21212
6 3626 130 938 16 93979
7 845 36 287 12 15524
8 643 32 149 13 16182
9 1181 33 311 15 19238
10 1835 64 617 13 28909
11 855 35 262 6 22357
12 1245 46 385 16 25560
13 993 69 369 7 9954
14 1685 61 558 12 18490
15 741 24 220 9 17777
16 854 38 313 10 25268
17 949 34 229 16 37525
18 331 20 88 5 6023
19 1602 54 494 20 25042
20 510 16 153 7 35713
21 628 37 234 13 7039
22 1279 51 361 13 40841
23 767 28 280 11 9214
24 1156 32 331 9 17446
25 1120 51 378 10 10295
26 622 12 227 7 13206
27 1203 98 396 13 26093
28 745 53 179 15 20744
29 1535 61 509 13 68013
30 1234 25 504 7 12840
31 757 28 225 14 12672
32 1087 23 366 11 10872
33 1105 60 341 3 21325
34 592 40 171 8 24542
35 1305 29 437 12 16401
36 0 0 0 0 0
37 705 33 313 12 12821
38 1188 34 366 8 14662
39 1110 21 232 20 22190
40 1094 35 389 18 37929
41 1062 34 340 9 18009
42 748 26 316 14 11076
43 404 12 140 7 24981
44 1076 45 419 13 30691
45 673 29 226 11 29164
46 517 36 161 11 13985
47 354 13 103 14 7588
48 1011 54 356 9 20023
49 890 39 293 12 25524
50 1067 28 414 11 14717
51 518 21 156 17 6832
52 697 36 189 10 9624
53 1095 44 442 11 24300
54 928 44 321 12 21790
55 1008 33 367 17 16493
56 951 30 309 6 9269
57 779 27 235 8 20105
58 439 12 137 12 11216
59 568 37 194 13 15569
60 614 24 220 14 21799
61 499 21 149 17 3772
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X_1t X_2t X_3t X_4t
-89.689118 3.596519 2.370210 5.441360 0.005912
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-273.98 -82.36 17.42 68.85 382.20
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -89.689118 52.953339 -1.694 0.095871 .
X_1t 3.596519 1.323199 2.718 0.008724 **
X_2t 2.370210 0.181401 13.066 < 2e-16 ***
X_3t 5.441360 4.379563 1.242 0.219251
X_4t 0.005912 0.001580 3.742 0.000431 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 134.1 on 56 degrees of freedom
Multiple R-squared: 0.9384, Adjusted R-squared: 0.934
F-statistic: 213.1 on 4 and 56 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.03561566 7.123132e-02 9.643843e-01
[2,] 0.59039404 8.192119e-01 4.096060e-01
[3,] 0.44872757 8.974551e-01 5.512724e-01
[4,] 0.33320086 6.664017e-01 6.667991e-01
[5,] 0.27524937 5.504987e-01 7.247506e-01
[6,] 0.42645587 8.529117e-01 5.735441e-01
[7,] 0.42242632 8.448526e-01 5.775737e-01
[8,] 0.32692141 6.538428e-01 6.730786e-01
[9,] 0.54786882 9.042624e-01 4.521312e-01
[10,] 0.55950740 8.809852e-01 4.404926e-01
[11,] 0.48783788 9.756758e-01 5.121621e-01
[12,] 0.44794669 8.958934e-01 5.520533e-01
[13,] 0.48176478 9.635296e-01 5.182352e-01
[14,] 0.46652071 9.330414e-01 5.334793e-01
[15,] 0.48178106 9.635621e-01 5.182189e-01
[16,] 0.40317306 8.063461e-01 5.968269e-01
[17,] 0.60097914 7.980417e-01 3.990209e-01
[18,] 0.52590143 9.481971e-01 4.740986e-01
[19,] 0.44898170 8.979634e-01 5.510183e-01
[20,] 0.59773373 8.045325e-01 4.022663e-01
[21,] 0.53436792 9.312642e-01 4.656321e-01
[22,] 0.82886347 3.422731e-01 1.711365e-01
[23,] 0.80303473 3.939305e-01 1.969653e-01
[24,] 0.75269887 4.946023e-01 2.473011e-01
[25,] 0.71106950 5.778610e-01 2.889305e-01
[26,] 0.69118885 6.176223e-01 3.088112e-01
[27,] 0.62946624 7.410675e-01 3.705338e-01
[28,] 0.60005553 7.998889e-01 3.999445e-01
[29,] 0.55163453 8.967309e-01 4.483655e-01
[30,] 0.67028674 6.594265e-01 3.297133e-01
[31,] 0.74525378 5.094924e-01 2.547462e-01
[32,] 0.99923829 1.523416e-03 7.617082e-04
[33,] 0.99975545 4.891041e-04 2.445520e-04
[34,] 0.99986437 2.712615e-04 1.356307e-04
[35,] 0.99995499 9.001298e-05 4.500649e-05
[36,] 0.99993218 1.356431e-04 6.782157e-05
[37,] 0.99987253 2.549441e-04 1.274721e-04
[38,] 0.99964296 7.140831e-04 3.570416e-04
[39,] 0.99939194 1.216120e-03 6.080600e-04
[40,] 0.99867967 2.640651e-03 1.320326e-03
[41,] 0.99675851 6.482987e-03 3.241493e-03
[42,] 0.99523151 9.536977e-03 4.768489e-03
[43,] 0.98855572 2.288856e-02 1.144428e-02
[44,] 0.96855314 6.289371e-02 3.144686e-02
[45,] 0.93518924 1.296215e-01 6.481076e-02
[46,] 0.96650128 6.699745e-02 3.349872e-02
> postscript(file="/var/wessaorg/rcomp/tmp/1yq161322145706.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/2pxmc1322145706.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/3vnni1322145706.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/4q4sq1322145706.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/5llq41322145706.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 = 61
Frequency = 1
1 2 3 4 5 6
99.124360 26.953498 40.220519 220.207343 -84.966511 382.200555
7 8 9 10 11 12
-32.113146 98.030379 219.509397 -9.561168 32.988773 18.540837
13 14 15 16 17 18
-137.017805 58.111304 68.853083 -138.657346 64.712490 77.364288
19 20 21 22 23 24
69.712832 -69.729050 -82.364914 17.423144 -22.002223 193.944516
25 26 27 28 29 30
15.047505 14.327192 -223.377579 15.543026 -273.979298 -74.811966
31 32 33 34 35 36
61.590905 102.339862 28.254644 -56.105538 92.345936 89.689118
37 38 39 40 41 42
-206.968401 157.695021 334.254650 -186.388935 68.091020 -146.469341
43 44 45 46 47 48
-67.080581 -241.461182 -109.555638 -46.926413 31.761945 -104.669881
49 50 51 52 53 54
-71.245898 -72.145214 29.514184 97.932131 -224.711955 -95.518253
55 56 57 58 59 60
-80.876133 112.950266 52.188145 29.204593 -97.987508 -109.132556
61
45.196968
> postscript(file="/var/wessaorg/rcomp/tmp/6a4e41322145706.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 99.124360 NA
1 26.953498 99.124360
2 40.220519 26.953498
3 220.207343 40.220519
4 -84.966511 220.207343
5 382.200555 -84.966511
6 -32.113146 382.200555
7 98.030379 -32.113146
8 219.509397 98.030379
9 -9.561168 219.509397
10 32.988773 -9.561168
11 18.540837 32.988773
12 -137.017805 18.540837
13 58.111304 -137.017805
14 68.853083 58.111304
15 -138.657346 68.853083
16 64.712490 -138.657346
17 77.364288 64.712490
18 69.712832 77.364288
19 -69.729050 69.712832
20 -82.364914 -69.729050
21 17.423144 -82.364914
22 -22.002223 17.423144
23 193.944516 -22.002223
24 15.047505 193.944516
25 14.327192 15.047505
26 -223.377579 14.327192
27 15.543026 -223.377579
28 -273.979298 15.543026
29 -74.811966 -273.979298
30 61.590905 -74.811966
31 102.339862 61.590905
32 28.254644 102.339862
33 -56.105538 28.254644
34 92.345936 -56.105538
35 89.689118 92.345936
36 -206.968401 89.689118
37 157.695021 -206.968401
38 334.254650 157.695021
39 -186.388935 334.254650
40 68.091020 -186.388935
41 -146.469341 68.091020
42 -67.080581 -146.469341
43 -241.461182 -67.080581
44 -109.555638 -241.461182
45 -46.926413 -109.555638
46 31.761945 -46.926413
47 -104.669881 31.761945
48 -71.245898 -104.669881
49 -72.145214 -71.245898
50 29.514184 -72.145214
51 97.932131 29.514184
52 -224.711955 97.932131
53 -95.518253 -224.711955
54 -80.876133 -95.518253
55 112.950266 -80.876133
56 52.188145 112.950266
57 29.204593 52.188145
58 -97.987508 29.204593
59 -109.132556 -97.987508
60 45.196968 -109.132556
61 NA 45.196968
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 26.953498 99.124360
[2,] 40.220519 26.953498
[3,] 220.207343 40.220519
[4,] -84.966511 220.207343
[5,] 382.200555 -84.966511
[6,] -32.113146 382.200555
[7,] 98.030379 -32.113146
[8,] 219.509397 98.030379
[9,] -9.561168 219.509397
[10,] 32.988773 -9.561168
[11,] 18.540837 32.988773
[12,] -137.017805 18.540837
[13,] 58.111304 -137.017805
[14,] 68.853083 58.111304
[15,] -138.657346 68.853083
[16,] 64.712490 -138.657346
[17,] 77.364288 64.712490
[18,] 69.712832 77.364288
[19,] -69.729050 69.712832
[20,] -82.364914 -69.729050
[21,] 17.423144 -82.364914
[22,] -22.002223 17.423144
[23,] 193.944516 -22.002223
[24,] 15.047505 193.944516
[25,] 14.327192 15.047505
[26,] -223.377579 14.327192
[27,] 15.543026 -223.377579
[28,] -273.979298 15.543026
[29,] -74.811966 -273.979298
[30,] 61.590905 -74.811966
[31,] 102.339862 61.590905
[32,] 28.254644 102.339862
[33,] -56.105538 28.254644
[34,] 92.345936 -56.105538
[35,] 89.689118 92.345936
[36,] -206.968401 89.689118
[37,] 157.695021 -206.968401
[38,] 334.254650 157.695021
[39,] -186.388935 334.254650
[40,] 68.091020 -186.388935
[41,] -146.469341 68.091020
[42,] -67.080581 -146.469341
[43,] -241.461182 -67.080581
[44,] -109.555638 -241.461182
[45,] -46.926413 -109.555638
[46,] 31.761945 -46.926413
[47,] -104.669881 31.761945
[48,] -71.245898 -104.669881
[49,] -72.145214 -71.245898
[50,] 29.514184 -72.145214
[51,] 97.932131 29.514184
[52,] -224.711955 97.932131
[53,] -95.518253 -224.711955
[54,] -80.876133 -95.518253
[55,] 112.950266 -80.876133
[56,] 52.188145 112.950266
[57,] 29.204593 52.188145
[58,] -97.987508 29.204593
[59,] -109.132556 -97.987508
[60,] 45.196968 -109.132556
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 26.953498 99.124360
2 40.220519 26.953498
3 220.207343 40.220519
4 -84.966511 220.207343
5 382.200555 -84.966511
6 -32.113146 382.200555
7 98.030379 -32.113146
8 219.509397 98.030379
9 -9.561168 219.509397
10 32.988773 -9.561168
11 18.540837 32.988773
12 -137.017805 18.540837
13 58.111304 -137.017805
14 68.853083 58.111304
15 -138.657346 68.853083
16 64.712490 -138.657346
17 77.364288 64.712490
18 69.712832 77.364288
19 -69.729050 69.712832
20 -82.364914 -69.729050
21 17.423144 -82.364914
22 -22.002223 17.423144
23 193.944516 -22.002223
24 15.047505 193.944516
25 14.327192 15.047505
26 -223.377579 14.327192
27 15.543026 -223.377579
28 -273.979298 15.543026
29 -74.811966 -273.979298
30 61.590905 -74.811966
31 102.339862 61.590905
32 28.254644 102.339862
33 -56.105538 28.254644
34 92.345936 -56.105538
35 89.689118 92.345936
36 -206.968401 89.689118
37 157.695021 -206.968401
38 334.254650 157.695021
39 -186.388935 334.254650
40 68.091020 -186.388935
41 -146.469341 68.091020
42 -67.080581 -146.469341
43 -241.461182 -67.080581
44 -109.555638 -241.461182
45 -46.926413 -109.555638
46 31.761945 -46.926413
47 -104.669881 31.761945
48 -71.245898 -104.669881
49 -72.145214 -71.245898
50 29.514184 -72.145214
51 97.932131 29.514184
52 -224.711955 97.932131
53 -95.518253 -224.711955
54 -80.876133 -95.518253
55 112.950266 -80.876133
56 52.188145 112.950266
57 29.204593 52.188145
58 -97.987508 29.204593
59 -109.132556 -97.987508
60 45.196968 -109.132556
> 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/746hq1322145706.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/83u121322145706.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/9k38a1322145706.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/100unv1322145706.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/11iae91322145706.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/12f6001322145706.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/13rq411322145706.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/14vj6z1322145706.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/15qnqt1322145706.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/16irv21322145706.tab")
+ }
>
> try(system("convert tmp/1yq161322145706.ps tmp/1yq161322145706.png",intern=TRUE))
character(0)
> try(system("convert tmp/2pxmc1322145706.ps tmp/2pxmc1322145706.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vnni1322145706.ps tmp/3vnni1322145706.png",intern=TRUE))
character(0)
> try(system("convert tmp/4q4sq1322145706.ps tmp/4q4sq1322145706.png",intern=TRUE))
character(0)
> try(system("convert tmp/5llq41322145706.ps tmp/5llq41322145706.png",intern=TRUE))
character(0)
> try(system("convert tmp/6a4e41322145706.ps tmp/6a4e41322145706.png",intern=TRUE))
character(0)
> try(system("convert tmp/746hq1322145706.ps tmp/746hq1322145706.png",intern=TRUE))
character(0)
> try(system("convert tmp/83u121322145706.ps tmp/83u121322145706.png",intern=TRUE))
character(0)
> try(system("convert tmp/9k38a1322145706.ps tmp/9k38a1322145706.png",intern=TRUE))
character(0)
> try(system("convert tmp/100unv1322145706.ps tmp/100unv1322145706.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.223 0.471 3.731