R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(521
+ ,18308
+ ,185
+ ,4.041
+ ,7.2
+ ,367
+ ,1148
+ ,600
+ ,0.55
+ ,8.5
+ ,443
+ ,18068
+ ,372
+ ,3.665
+ ,5.7
+ ,365
+ ,7729
+ ,142
+ ,2.351
+ ,7.3
+ ,614
+ ,100484
+ ,432
+ ,29.76
+ ,7.5
+ ,385
+ ,16728
+ ,290
+ ,3.294
+ ,5
+ ,286
+ ,14630
+ ,346
+ ,3.287
+ ,6.7
+ ,397
+ ,4008
+ ,328
+ ,0.666
+ ,6.2
+ ,764
+ ,38927
+ ,354
+ ,12.938
+ ,7.3
+ ,427
+ ,22322
+ ,266
+ ,6.478
+ ,5
+ ,153
+ ,3711
+ ,320
+ ,1.108
+ ,2.8
+ ,231
+ ,3136
+ ,197
+ ,1.007
+ ,6.1
+ ,524
+ ,50508
+ ,266
+ ,11.431
+ ,7.1
+ ,328
+ ,28886
+ ,173
+ ,5.544
+ ,5.9
+ ,240
+ ,16996
+ ,190
+ ,2.777
+ ,4.6
+ ,286
+ ,13035
+ ,239
+ ,2.478
+ ,4.4
+ ,285
+ ,12973
+ ,190
+ ,3.685
+ ,7.4
+ ,569
+ ,16309
+ ,241
+ ,4.22
+ ,7.1
+ ,96
+ ,5227
+ ,189
+ ,1.228
+ ,7.5
+ ,498
+ ,19235
+ ,358
+ ,4.781
+ ,5.9
+ ,481
+ ,44487
+ ,315
+ ,6.016
+ ,9
+ ,468
+ ,44213
+ ,303
+ ,9.295
+ ,9.2
+ ,177
+ ,23619
+ ,228
+ ,4.375
+ ,5.1
+ ,198
+ ,9106
+ ,134
+ ,2.573
+ ,8.6
+ ,458
+ ,24917
+ ,189
+ ,5.117
+ ,6.6
+ ,108
+ ,3872
+ ,196
+ ,0.799
+ ,6.9
+ ,246
+ ,8945
+ ,183
+ ,1.578
+ ,2.7
+ ,291
+ ,2373
+ ,417
+ ,1.202
+ ,5.5
+ ,68
+ ,7128
+ ,233
+ ,1.109
+ ,7.2
+ ,311
+ ,23624
+ ,349
+ ,7.73
+ ,6.6
+ ,606
+ ,5242
+ ,284
+ ,1.515
+ ,6.9
+ ,512
+ ,92629
+ ,499
+ ,17.99
+ ,7.2
+ ,426
+ ,28795
+ ,231
+ ,6.629
+ ,5.8
+ ,47
+ ,4487
+ ,143
+ ,0.639
+ ,4.1
+ ,265
+ ,48799
+ ,249
+ ,10.847
+ ,6.4
+ ,370
+ ,14067
+ ,195
+ ,3.146
+ ,6.7
+ ,312
+ ,12693
+ ,288
+ ,2.842
+ ,6
+ ,222
+ ,62184
+ ,229
+ ,11.882
+ ,6.9
+ ,280
+ ,9153
+ ,287
+ ,1.003
+ ,8.5
+ ,759
+ ,14250
+ ,224
+ ,3.487
+ ,6.2
+ ,114
+ ,3680
+ ,161
+ ,0.696
+ ,3.4
+ ,419
+ ,18063
+ ,221
+ ,4.877
+ ,6.6
+ ,435
+ ,65112
+ ,237
+ ,16.987
+ ,6.6
+ ,186
+ ,11340
+ ,220
+ ,1.723
+ ,4.9
+ ,87
+ ,4553
+ ,185
+ ,0.563
+ ,6.4
+ ,188
+ ,28960
+ ,260
+ ,6.187
+ ,5.8
+ ,303
+ ,19201
+ ,261
+ ,4.867
+ ,6.3
+ ,102
+ ,7533
+ ,118
+ ,1.793
+ ,10.5
+ ,127
+ ,26343
+ ,268
+ ,4.892
+ ,5.4
+ ,251
+ ,1641
+ ,300
+ ,0.454
+ ,5.1
+ ,205
+ ,145360
+ ,237
+ ,10.379
+ ,6.8
+ ,453
+ ,9066420
+ ,240
+ ,82.422
+ ,5.6
+ ,320
+ ,1038933
+ ,185
+ ,16.491
+ ,3.8
+ ,405
+ ,2739420
+ ,201
+ ,60.876
+ ,8.2
+ ,89
+ ,61620
+ ,193
+ ,0.474
+ ,4.1
+ ,74
+ ,827530
+ ,254
+ ,7.523
+ ,2.8
+ ,101
+ ,534100
+ ,230
+ ,5.45
+ ,6.3
+ ,321
+ ,328755
+ ,197
+ ,10.605
+ ,11.4
+ ,315
+ ,1413895
+ ,248
+ ,40.397
+ ,19.4
+ ,229
+ ,2909136
+ ,258
+ ,60.607
+ ,5.8
+ ,302
+ ,3604246
+ ,206
+ ,58.133
+ ,6.9
+ ,216
+ ,917504
+ ,199
+ ,8.192
+ ,3.5)
+ ,dim=c(5
+ ,62)
+ ,dimnames=list(c('Assaults'
+ ,'BachDegrees'
+ ,'PoliceExp'
+ ,'Popul'
+ ,'Unempl')
+ ,1:62))
> y <- array(NA,dim=c(5,62),dimnames=list(c('Assaults','BachDegrees','PoliceExp','Popul','Unempl'),1:62))
> 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
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, as.Date.numeric
> 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
Assaults BachDegrees PoliceExp Popul Unempl
1 521 18308 185 4.041 7.2
2 367 1148 600 0.550 8.5
3 443 18068 372 3.665 5.7
4 365 7729 142 2.351 7.3
5 614 100484 432 29.760 7.5
6 385 16728 290 3.294 5.0
7 286 14630 346 3.287 6.7
8 397 4008 328 0.666 6.2
9 764 38927 354 12.938 7.3
10 427 22322 266 6.478 5.0
11 153 3711 320 1.108 2.8
12 231 3136 197 1.007 6.1
13 524 50508 266 11.431 7.1
14 328 28886 173 5.544 5.9
15 240 16996 190 2.777 4.6
16 286 13035 239 2.478 4.4
17 285 12973 190 3.685 7.4
18 569 16309 241 4.220 7.1
19 96 5227 189 1.228 7.5
20 498 19235 358 4.781 5.9
21 481 44487 315 6.016 9.0
22 468 44213 303 9.295 9.2
23 177 23619 228 4.375 5.1
24 198 9106 134 2.573 8.6
25 458 24917 189 5.117 6.6
26 108 3872 196 0.799 6.9
27 246 8945 183 1.578 2.7
28 291 2373 417 1.202 5.5
29 68 7128 233 1.109 7.2
30 311 23624 349 7.730 6.6
31 606 5242 284 1.515 6.9
32 512 92629 499 17.990 7.2
33 426 28795 231 6.629 5.8
34 47 4487 143 0.639 4.1
35 265 48799 249 10.847 6.4
36 370 14067 195 3.146 6.7
37 312 12693 288 2.842 6.0
38 222 62184 229 11.882 6.9
39 280 9153 287 1.003 8.5
40 759 14250 224 3.487 6.2
41 114 3680 161 0.696 3.4
42 419 18063 221 4.877 6.6
43 435 65112 237 16.987 6.6
44 186 11340 220 1.723 4.9
45 87 4553 185 0.563 6.4
46 188 28960 260 6.187 5.8
47 303 19201 261 4.867 6.3
48 102 7533 118 1.793 10.5
49 127 26343 268 4.892 5.4
50 251 1641 300 0.454 5.1
51 205 145360 237 10.379 6.8
52 453 9066420 240 82.422 5.6
53 320 1038933 185 16.491 3.8
54 405 2739420 201 60.876 8.2
55 89 61620 193 0.474 4.1
56 74 827530 254 7.523 2.8
57 101 534100 230 5.450 6.3
58 321 328755 197 10.605 11.4
59 315 1413895 248 40.397 19.4
60 229 2909136 258 60.607 5.8
61 302 3604246 206 58.133 6.9
62 216 917504 199 8.192 3.5
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) BachDegrees PoliceExp Popul Unempl
6.984e+01 -3.012e-05 6.863e-01 3.750e+00 7.334e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-218.50 -111.39 -13.82 91.57 477.31
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.984e+01 8.072e+01 0.865 0.39052
BachDegrees -3.012e-05 3.354e-05 -0.898 0.37290
PoliceExp 6.863e-01 2.306e-01 2.976 0.00428 **
Popul 3.750e+00 2.753e+00 1.362 0.17850
Unempl 7.334e+00 8.986e+00 0.816 0.41780
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 152.6 on 57 degrees of freedom
Multiple R-squared: 0.2146, Adjusted R-squared: 0.1595
F-statistic: 3.894 on 4 and 57 DF, p-value: 0.007283
> 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.1882753 0.37655066 0.81172467
[2,] 0.4532425 0.90648500 0.54675750
[3,] 0.3682861 0.73657226 0.63171387
[4,] 0.4093291 0.81865820 0.59067090
[5,] 0.4235461 0.84709221 0.57645390
[6,] 0.3377318 0.67546369 0.66226816
[7,] 0.2661484 0.53229683 0.73385159
[8,] 0.1954344 0.39086873 0.80456563
[9,] 0.1307566 0.26151325 0.86924337
[10,] 0.1327394 0.26547878 0.86726061
[11,] 0.1836808 0.36736160 0.81631920
[12,] 0.3715427 0.74308533 0.62845733
[13,] 0.3448088 0.68961759 0.65519121
[14,] 0.2827644 0.56552872 0.71723564
[15,] 0.2340548 0.46810967 0.76594516
[16,] 0.2277855 0.45557096 0.77221452
[17,] 0.2210043 0.44200857 0.77899572
[18,] 0.2332798 0.46655961 0.76672020
[19,] 0.2562802 0.51256037 0.74371981
[20,] 0.1988200 0.39763998 0.80118001
[21,] 0.1592175 0.31843497 0.84078251
[22,] 0.2334241 0.46684828 0.76657586
[23,] 0.1965286 0.39305719 0.80347141
[24,] 0.4131976 0.82639522 0.58680239
[25,] 0.3611438 0.72228760 0.63885620
[26,] 0.3498861 0.69977215 0.65011392
[27,] 0.3548660 0.70973206 0.64513397
[28,] 0.3178138 0.63562756 0.68218622
[29,] 0.2891057 0.57821146 0.71089427
[30,] 0.2276107 0.45522135 0.77238932
[31,] 0.2056003 0.41120053 0.79439974
[32,] 0.1563950 0.31278997 0.84360502
[33,] 0.9359577 0.12808459 0.06404229
[34,] 0.9097122 0.18057569 0.09028785
[35,] 0.9556885 0.08862294 0.04431147
[36,] 0.9853735 0.02925306 0.01462653
[37,] 0.9745024 0.05099517 0.02549758
[38,] 0.9687089 0.06258218 0.03129109
[39,] 0.9490761 0.10184780 0.05092390
[40,] 0.9503666 0.09926687 0.04963344
[41,] 0.9616681 0.07666375 0.03833188
[42,] 0.9357443 0.12851135 0.06425567
[43,] 0.9871753 0.02564942 0.01282471
[44,] 0.9815091 0.03698189 0.01849095
[45,] 0.9641971 0.07160581 0.03580291
[46,] 0.9495100 0.10097992 0.05048996
[47,] 0.8964077 0.20718464 0.10359232
> postscript(file="/var/wessaorg/rcomp/tmp/1xj3n1355044503.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/28auo1355044503.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/3rexq1355044503.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/4i7e71355044503.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/5u9gn1355044503.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 = 62
Frequency = 1
1 2 3 4 5 6
256.782529 -178.990344 62.848540 135.579151 84.090522 67.609737
7 8 9 10 11 12
-82.327648 54.202032 350.320277 114.308766 -161.038631 -22.464575
13 14 15 16 17 18
178.181680 76.234327 -3.879924 10.960256 17.058625 266.351726
19 20 21 22 23 24
-163.007553 121.839855 107.744483 89.208128 -102.419570 -36.255686
25 26 27 28 29 30
191.601566 -149.843240 25.112687 -109.804090 -218.500957 -75.044365
31 32 33 34 35 36
285.118634 -17.788073 131.090720 -153.315773 -61.878834 105.815215
37 38 39 40 41 42
-9.778076 -98.298150 -52.636707 477.306242 -93.773434 131.333627
43 44 45 46 47 48
92.355100 -76.886715 -158.721419 -125.149375 -9.846337 -132.331674
49 50 51 52 53 54
-183.928496 -63.790506 -111.913358 141.353213 64.769044 -8.715347
55 56 57 58 59 60
-143.291075 -193.987230 -177.248831 2.478599 -176.236005 -200.112229
61 62
-69.277567 -19.173487
> postscript(file="/var/wessaorg/rcomp/tmp/67cwu1355044503.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 = 62
Frequency = 1
lag(myerror, k = 1) myerror
0 256.782529 NA
1 -178.990344 256.782529
2 62.848540 -178.990344
3 135.579151 62.848540
4 84.090522 135.579151
5 67.609737 84.090522
6 -82.327648 67.609737
7 54.202032 -82.327648
8 350.320277 54.202032
9 114.308766 350.320277
10 -161.038631 114.308766
11 -22.464575 -161.038631
12 178.181680 -22.464575
13 76.234327 178.181680
14 -3.879924 76.234327
15 10.960256 -3.879924
16 17.058625 10.960256
17 266.351726 17.058625
18 -163.007553 266.351726
19 121.839855 -163.007553
20 107.744483 121.839855
21 89.208128 107.744483
22 -102.419570 89.208128
23 -36.255686 -102.419570
24 191.601566 -36.255686
25 -149.843240 191.601566
26 25.112687 -149.843240
27 -109.804090 25.112687
28 -218.500957 -109.804090
29 -75.044365 -218.500957
30 285.118634 -75.044365
31 -17.788073 285.118634
32 131.090720 -17.788073
33 -153.315773 131.090720
34 -61.878834 -153.315773
35 105.815215 -61.878834
36 -9.778076 105.815215
37 -98.298150 -9.778076
38 -52.636707 -98.298150
39 477.306242 -52.636707
40 -93.773434 477.306242
41 131.333627 -93.773434
42 92.355100 131.333627
43 -76.886715 92.355100
44 -158.721419 -76.886715
45 -125.149375 -158.721419
46 -9.846337 -125.149375
47 -132.331674 -9.846337
48 -183.928496 -132.331674
49 -63.790506 -183.928496
50 -111.913358 -63.790506
51 141.353213 -111.913358
52 64.769044 141.353213
53 -8.715347 64.769044
54 -143.291075 -8.715347
55 -193.987230 -143.291075
56 -177.248831 -193.987230
57 2.478599 -177.248831
58 -176.236005 2.478599
59 -200.112229 -176.236005
60 -69.277567 -200.112229
61 -19.173487 -69.277567
62 NA -19.173487
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -178.990344 256.782529
[2,] 62.848540 -178.990344
[3,] 135.579151 62.848540
[4,] 84.090522 135.579151
[5,] 67.609737 84.090522
[6,] -82.327648 67.609737
[7,] 54.202032 -82.327648
[8,] 350.320277 54.202032
[9,] 114.308766 350.320277
[10,] -161.038631 114.308766
[11,] -22.464575 -161.038631
[12,] 178.181680 -22.464575
[13,] 76.234327 178.181680
[14,] -3.879924 76.234327
[15,] 10.960256 -3.879924
[16,] 17.058625 10.960256
[17,] 266.351726 17.058625
[18,] -163.007553 266.351726
[19,] 121.839855 -163.007553
[20,] 107.744483 121.839855
[21,] 89.208128 107.744483
[22,] -102.419570 89.208128
[23,] -36.255686 -102.419570
[24,] 191.601566 -36.255686
[25,] -149.843240 191.601566
[26,] 25.112687 -149.843240
[27,] -109.804090 25.112687
[28,] -218.500957 -109.804090
[29,] -75.044365 -218.500957
[30,] 285.118634 -75.044365
[31,] -17.788073 285.118634
[32,] 131.090720 -17.788073
[33,] -153.315773 131.090720
[34,] -61.878834 -153.315773
[35,] 105.815215 -61.878834
[36,] -9.778076 105.815215
[37,] -98.298150 -9.778076
[38,] -52.636707 -98.298150
[39,] 477.306242 -52.636707
[40,] -93.773434 477.306242
[41,] 131.333627 -93.773434
[42,] 92.355100 131.333627
[43,] -76.886715 92.355100
[44,] -158.721419 -76.886715
[45,] -125.149375 -158.721419
[46,] -9.846337 -125.149375
[47,] -132.331674 -9.846337
[48,] -183.928496 -132.331674
[49,] -63.790506 -183.928496
[50,] -111.913358 -63.790506
[51,] 141.353213 -111.913358
[52,] 64.769044 141.353213
[53,] -8.715347 64.769044
[54,] -143.291075 -8.715347
[55,] -193.987230 -143.291075
[56,] -177.248831 -193.987230
[57,] 2.478599 -177.248831
[58,] -176.236005 2.478599
[59,] -200.112229 -176.236005
[60,] -69.277567 -200.112229
[61,] -19.173487 -69.277567
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -178.990344 256.782529
2 62.848540 -178.990344
3 135.579151 62.848540
4 84.090522 135.579151
5 67.609737 84.090522
6 -82.327648 67.609737
7 54.202032 -82.327648
8 350.320277 54.202032
9 114.308766 350.320277
10 -161.038631 114.308766
11 -22.464575 -161.038631
12 178.181680 -22.464575
13 76.234327 178.181680
14 -3.879924 76.234327
15 10.960256 -3.879924
16 17.058625 10.960256
17 266.351726 17.058625
18 -163.007553 266.351726
19 121.839855 -163.007553
20 107.744483 121.839855
21 89.208128 107.744483
22 -102.419570 89.208128
23 -36.255686 -102.419570
24 191.601566 -36.255686
25 -149.843240 191.601566
26 25.112687 -149.843240
27 -109.804090 25.112687
28 -218.500957 -109.804090
29 -75.044365 -218.500957
30 285.118634 -75.044365
31 -17.788073 285.118634
32 131.090720 -17.788073
33 -153.315773 131.090720
34 -61.878834 -153.315773
35 105.815215 -61.878834
36 -9.778076 105.815215
37 -98.298150 -9.778076
38 -52.636707 -98.298150
39 477.306242 -52.636707
40 -93.773434 477.306242
41 131.333627 -93.773434
42 92.355100 131.333627
43 -76.886715 92.355100
44 -158.721419 -76.886715
45 -125.149375 -158.721419
46 -9.846337 -125.149375
47 -132.331674 -9.846337
48 -183.928496 -132.331674
49 -63.790506 -183.928496
50 -111.913358 -63.790506
51 141.353213 -111.913358
52 64.769044 141.353213
53 -8.715347 64.769044
54 -143.291075 -8.715347
55 -193.987230 -143.291075
56 -177.248831 -193.987230
57 2.478599 -177.248831
58 -176.236005 2.478599
59 -200.112229 -176.236005
60 -69.277567 -200.112229
61 -19.173487 -69.277567
> 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/70dz31355044503.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/8p85t1355044503.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/9qsp41355044503.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/103wte1355044503.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/118omr1355044503.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/12i7l41355044503.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/13pq761355044503.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/14w7kj1355044503.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/15ny451355044503.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/168q7j1355044503.tab")
+ }
>
> try(system("convert tmp/1xj3n1355044503.ps tmp/1xj3n1355044503.png",intern=TRUE))
character(0)
> try(system("convert tmp/28auo1355044503.ps tmp/28auo1355044503.png",intern=TRUE))
character(0)
> try(system("convert tmp/3rexq1355044503.ps tmp/3rexq1355044503.png",intern=TRUE))
character(0)
> try(system("convert tmp/4i7e71355044503.ps tmp/4i7e71355044503.png",intern=TRUE))
character(0)
> try(system("convert tmp/5u9gn1355044503.ps tmp/5u9gn1355044503.png",intern=TRUE))
character(0)
> try(system("convert tmp/67cwu1355044503.ps tmp/67cwu1355044503.png",intern=TRUE))
character(0)
> try(system("convert tmp/70dz31355044503.ps tmp/70dz31355044503.png",intern=TRUE))
character(0)
> try(system("convert tmp/8p85t1355044503.ps tmp/8p85t1355044503.png",intern=TRUE))
character(0)
> try(system("convert tmp/9qsp41355044503.ps tmp/9qsp41355044503.png",intern=TRUE))
character(0)
> try(system("convert tmp/103wte1355044503.ps tmp/103wte1355044503.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
8.460 1.334 9.845