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(1
+ ,521
+ ,18308
+ ,185
+ ,4.041
+ ,7.2
+ ,1
+ ,367
+ ,1148
+ ,600
+ ,0.55
+ ,8.5
+ ,1
+ ,443
+ ,18068
+ ,372
+ ,3.665
+ ,5.7
+ ,1
+ ,365
+ ,7729
+ ,142
+ ,2.351
+ ,7.3
+ ,1
+ ,614
+ ,100484
+ ,432
+ ,29.76
+ ,7.5
+ ,1
+ ,385
+ ,16728
+ ,290
+ ,3.294
+ ,5
+ ,1
+ ,286
+ ,14630
+ ,346
+ ,3.287
+ ,6.7
+ ,1
+ ,397
+ ,4008
+ ,328
+ ,0.666
+ ,6.2
+ ,1
+ ,764
+ ,38927
+ ,354
+ ,12.938
+ ,7.3
+ ,1
+ ,427
+ ,22322
+ ,266
+ ,6.478
+ ,5
+ ,1
+ ,153
+ ,3711
+ ,320
+ ,1.108
+ ,2.8
+ ,1
+ ,231
+ ,3136
+ ,197
+ ,1.007
+ ,6.1
+ ,1
+ ,524
+ ,50508
+ ,266
+ ,11.431
+ ,7.1
+ ,1
+ ,328
+ ,28886
+ ,173
+ ,5.544
+ ,5.9
+ ,1
+ ,240
+ ,16996
+ ,190
+ ,2.777
+ ,4.6
+ ,1
+ ,286
+ ,13035
+ ,239
+ ,2.478
+ ,4.4
+ ,1
+ ,285
+ ,12973
+ ,190
+ ,3.685
+ ,7.4
+ ,1
+ ,569
+ ,16309
+ ,241
+ ,4.22
+ ,7.1
+ ,1
+ ,96
+ ,5227
+ ,189
+ ,1.228
+ ,7.5
+ ,1
+ ,498
+ ,19235
+ ,358
+ ,4.781
+ ,5.9
+ ,1
+ ,481
+ ,44487
+ ,315
+ ,6.016
+ ,9
+ ,1
+ ,468
+ ,44213
+ ,303
+ ,9.295
+ ,9.2
+ ,1
+ ,177
+ ,23619
+ ,228
+ ,4.375
+ ,5.1
+ ,1
+ ,198
+ ,9106
+ ,134
+ ,2.573
+ ,8.6
+ ,1
+ ,458
+ ,24917
+ ,189
+ ,5.117
+ ,6.6
+ ,1
+ ,108
+ ,3872
+ ,196
+ ,0.799
+ ,6.9
+ ,1
+ ,246
+ ,8945
+ ,183
+ ,1.578
+ ,2.7
+ ,1
+ ,291
+ ,2373
+ ,417
+ ,1.202
+ ,5.5
+ ,1
+ ,68
+ ,7128
+ ,233
+ ,1.109
+ ,7.2
+ ,1
+ ,311
+ ,23624
+ ,349
+ ,7.73
+ ,6.6
+ ,1
+ ,606
+ ,5242
+ ,284
+ ,1.515
+ ,6.9
+ ,1
+ ,512
+ ,92629
+ ,499
+ ,17.99
+ ,7.2
+ ,1
+ ,426
+ ,28795
+ ,231
+ ,6.629
+ ,5.8
+ ,1
+ ,47
+ ,4487
+ ,143
+ ,0.639
+ ,4.1
+ ,1
+ ,265
+ ,48799
+ ,249
+ ,10.847
+ ,6.4
+ ,1
+ ,370
+ ,14067
+ ,195
+ ,3.146
+ ,6.7
+ ,1
+ ,312
+ ,12693
+ ,288
+ ,2.842
+ ,6
+ ,1
+ ,222
+ ,62184
+ ,229
+ ,11.882
+ ,6.9
+ ,1
+ ,280
+ ,9153
+ ,287
+ ,1.003
+ ,8.5
+ ,1
+ ,759
+ ,14250
+ ,224
+ ,3.487
+ ,6.2
+ ,1
+ ,114
+ ,3680
+ ,161
+ ,0.696
+ ,3.4
+ ,1
+ ,419
+ ,18063
+ ,221
+ ,4.877
+ ,6.6
+ ,1
+ ,435
+ ,65112
+ ,237
+ ,16.987
+ ,6.6
+ ,1
+ ,186
+ ,11340
+ ,220
+ ,1.723
+ ,4.9
+ ,1
+ ,87
+ ,4553
+ ,185
+ ,0.563
+ ,6.4
+ ,1
+ ,188
+ ,28960
+ ,260
+ ,6.187
+ ,5.8
+ ,1
+ ,303
+ ,19201
+ ,261
+ ,4.867
+ ,6.3
+ ,1
+ ,102
+ ,7533
+ ,118
+ ,1.793
+ ,10.5
+ ,1
+ ,127
+ ,26343
+ ,268
+ ,4.892
+ ,5.4
+ ,1
+ ,251
+ ,1641
+ ,300
+ ,0.454
+ ,5.1
+ ,0
+ ,205
+ ,145360
+ ,237
+ ,10.379
+ ,6.8
+ ,0
+ ,453
+ ,9066420
+ ,240
+ ,82.422
+ ,5.6
+ ,0
+ ,320
+ ,1038933
+ ,185
+ ,16.491
+ ,3.8
+ ,0
+ ,405
+ ,2739420
+ ,201
+ ,60.876
+ ,8.2
+ ,0
+ ,89
+ ,61620
+ ,193
+ ,0.474
+ ,4.1
+ ,0
+ ,74
+ ,827530
+ ,254
+ ,7.523
+ ,2.8
+ ,0
+ ,101
+ ,534100
+ ,230
+ ,5.45
+ ,6.3
+ ,0
+ ,321
+ ,328755
+ ,197
+ ,10.605
+ ,11.4
+ ,0
+ ,315
+ ,1413895
+ ,248
+ ,40.397
+ ,19.4
+ ,0
+ ,229
+ ,2909136
+ ,258
+ ,60.607
+ ,5.8
+ ,0
+ ,302
+ ,3604246
+ ,206
+ ,58.133
+ ,6.9
+ ,0
+ ,216
+ ,917504
+ ,199
+ ,8.192
+ ,3.5)
+ ,dim=c(6
+ ,62)
+ ,dimnames=list(c('Pop'
+ ,'Assaults'
+ ,'BachDgr'
+ ,'PoliceExp'
+ ,'TotlPopul'
+ ,'Unempl')
+ ,1:62))
> y <- array(NA,dim=c(6,62),dimnames=list(c('Pop','Assaults','BachDgr','PoliceExp','TotlPopul','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 = '2'
> 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 Pop BachDgr PoliceExp TotlPopul Unempl
1 521 1 18308 185 4.041 7.2
2 367 1 1148 600 0.550 8.5
3 443 1 18068 372 3.665 5.7
4 365 1 7729 142 2.351 7.3
5 614 1 100484 432 29.760 7.5
6 385 1 16728 290 3.294 5.0
7 286 1 14630 346 3.287 6.7
8 397 1 4008 328 0.666 6.2
9 764 1 38927 354 12.938 7.3
10 427 1 22322 266 6.478 5.0
11 153 1 3711 320 1.108 2.8
12 231 1 3136 197 1.007 6.1
13 524 1 50508 266 11.431 7.1
14 328 1 28886 173 5.544 5.9
15 240 1 16996 190 2.777 4.6
16 286 1 13035 239 2.478 4.4
17 285 1 12973 190 3.685 7.4
18 569 1 16309 241 4.220 7.1
19 96 1 5227 189 1.228 7.5
20 498 1 19235 358 4.781 5.9
21 481 1 44487 315 6.016 9.0
22 468 1 44213 303 9.295 9.2
23 177 1 23619 228 4.375 5.1
24 198 1 9106 134 2.573 8.6
25 458 1 24917 189 5.117 6.6
26 108 1 3872 196 0.799 6.9
27 246 1 8945 183 1.578 2.7
28 291 1 2373 417 1.202 5.5
29 68 1 7128 233 1.109 7.2
30 311 1 23624 349 7.730 6.6
31 606 1 5242 284 1.515 6.9
32 512 1 92629 499 17.990 7.2
33 426 1 28795 231 6.629 5.8
34 47 1 4487 143 0.639 4.1
35 265 1 48799 249 10.847 6.4
36 370 1 14067 195 3.146 6.7
37 312 1 12693 288 2.842 6.0
38 222 1 62184 229 11.882 6.9
39 280 1 9153 287 1.003 8.5
40 759 1 14250 224 3.487 6.2
41 114 1 3680 161 0.696 3.4
42 419 1 18063 221 4.877 6.6
43 435 1 65112 237 16.987 6.6
44 186 1 11340 220 1.723 4.9
45 87 1 4553 185 0.563 6.4
46 188 1 28960 260 6.187 5.8
47 303 1 19201 261 4.867 6.3
48 102 1 7533 118 1.793 10.5
49 127 1 26343 268 4.892 5.4
50 251 1 1641 300 0.454 5.1
51 205 0 145360 237 10.379 6.8
52 453 0 9066420 240 82.422 5.6
53 320 0 1038933 185 16.491 3.8
54 405 0 2739420 201 60.876 8.2
55 89 0 61620 193 0.474 4.1
56 74 0 827530 254 7.523 2.8
57 101 0 534100 230 5.450 6.3
58 321 0 328755 197 10.605 11.4
59 315 0 1413895 248 40.397 19.4
60 229 0 2909136 258 60.607 5.8
61 302 0 3604246 206 58.133 6.9
62 216 0 917504 199 8.192 3.5
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Pop BachDgr PoliceExp TotlPopul Unempl
-4.620e+01 1.575e+02 -2.243e-05 5.531e-01 5.526e+00 7.693e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-233.49 -98.33 -14.25 91.08 457.19
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.620e+01 8.925e+01 -0.518 0.6068
Pop 1.575e+02 6.121e+01 2.573 0.0128 *
BachDgr -2.243e-05 3.214e-05 -0.698 0.4880
PoliceExp 5.531e-01 2.260e-01 2.447 0.0176 *
TotlPopul 5.526e+00 2.716e+00 2.035 0.0466 *
Unempl 7.693e+00 8.574e+00 0.897 0.3735
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 145.6 on 56 degrees of freedom
Multiple R-squared: 0.2976, Adjusted R-squared: 0.2349
F-statistic: 4.746 on 5 and 56 DF, p-value: 0.001102
> 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.6431847 0.71363062 0.35681531
[2,] 0.5402855 0.91942905 0.45971452
[3,] 0.5698462 0.86030763 0.43015381
[4,] 0.5750154 0.84996923 0.42498461
[5,] 0.4746535 0.94930701 0.52534650
[6,] 0.3857126 0.77142517 0.61428742
[7,] 0.2961673 0.59233452 0.70383274
[8,] 0.2089918 0.41798366 0.79100817
[9,] 0.2077083 0.41541661 0.79229169
[10,] 0.2663813 0.53276268 0.73361866
[11,] 0.4810368 0.96207365 0.51896317
[12,] 0.4501722 0.90034437 0.54982782
[13,] 0.3781883 0.75637657 0.62181172
[14,] 0.3170063 0.63401252 0.68299374
[15,] 0.3113154 0.62263084 0.68868458
[16,] 0.3022732 0.60454637 0.69772682
[17,] 0.3037802 0.60756035 0.69621982
[18,] 0.3352439 0.67048775 0.66475613
[19,] 0.2648993 0.52979855 0.73510072
[20,] 0.2150143 0.43002869 0.78498566
[21,] 0.3127974 0.62559481 0.68720260
[22,] 0.2703099 0.54061989 0.72969005
[23,] 0.5008784 0.99824316 0.49912158
[24,] 0.4444368 0.88887356 0.55556322
[25,] 0.4170074 0.83401473 0.58299263
[26,] 0.4416263 0.88325251 0.55837375
[27,] 0.4016552 0.80331046 0.59834477
[28,] 0.3529792 0.70595847 0.64702077
[29,] 0.2811009 0.56220185 0.71889908
[30,] 0.2614509 0.52290182 0.73854909
[31,] 0.2028817 0.40576333 0.79711834
[32,] 0.9390080 0.12198397 0.06099199
[33,] 0.9179559 0.16408812 0.08204406
[34,] 0.9497269 0.10054627 0.05027313
[35,] 0.9782009 0.04359813 0.02179906
[36,] 0.9623978 0.07520450 0.03760225
[37,] 0.9547176 0.09056475 0.04528237
[38,] 0.9273935 0.14521299 0.07260650
[39,] 0.9245687 0.15086268 0.07543134
[40,] 0.9739610 0.05207807 0.02603904
[41,] 0.9886867 0.02262659 0.01131330
[42,] 0.9714620 0.05707597 0.02853798
[43,] 0.9565607 0.08687860 0.04343930
[44,] 0.9132591 0.17348190 0.08674095
[45,] 0.8725730 0.25485392 0.12742696
> postscript(file="/var/fisher/rcomp/tmp/1r06a1354819135.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/fisher/rcomp/tmp/24yms1354819135.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/fisher/rcomp/tmp/3fv6m1354819135.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/fisher/rcomp/tmp/4rn951354819135.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/fisher/rcomp/tmp/5dyl71354819135.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
230.101727 -144.529850 62.283211 106.218012 43.889999 57.042973
7 8 9 10 11 12
-86.017209 53.031250 330.153031 94.847249 -162.837226 -41.646879
13 14 15 16 17 18
148.952666 45.671068 -26.706366 -4.706363 -8.354698 246.863157
19 20 21 22 23 24
-184.166559 117.346918 94.023792 67.995603 -123.253217 -67.553918
25 26 27 28 29 30
163.707043 -169.082144 4.227350 -99.808706 -233.494975 -86.258770
31 32 33 34 35 36
276.318628 -27.991864 106.362274 -178.329423 -92.070499 82.268048
37 38 39 40 41 42
-20.136386 -133.274337 -60.732183 457.194144 -116.233340 108.180276
43 44 45 46 47 48
49.462690 -93.909532 -178.832054 -145.231416 -27.555206 -165.045612
49 50 51 52 53 54
-200.481235 -67.901140 13.701477 71.266824 166.810267 1.973891
55 56 57 58 59 60
-4.331654 -64.843546 -46.620756 119.303241 -116.745097 -181.796224
61 62
-59.231809 100.513385
> postscript(file="/var/fisher/rcomp/tmp/6zwjs1354819135.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 230.101727 NA
1 -144.529850 230.101727
2 62.283211 -144.529850
3 106.218012 62.283211
4 43.889999 106.218012
5 57.042973 43.889999
6 -86.017209 57.042973
7 53.031250 -86.017209
8 330.153031 53.031250
9 94.847249 330.153031
10 -162.837226 94.847249
11 -41.646879 -162.837226
12 148.952666 -41.646879
13 45.671068 148.952666
14 -26.706366 45.671068
15 -4.706363 -26.706366
16 -8.354698 -4.706363
17 246.863157 -8.354698
18 -184.166559 246.863157
19 117.346918 -184.166559
20 94.023792 117.346918
21 67.995603 94.023792
22 -123.253217 67.995603
23 -67.553918 -123.253217
24 163.707043 -67.553918
25 -169.082144 163.707043
26 4.227350 -169.082144
27 -99.808706 4.227350
28 -233.494975 -99.808706
29 -86.258770 -233.494975
30 276.318628 -86.258770
31 -27.991864 276.318628
32 106.362274 -27.991864
33 -178.329423 106.362274
34 -92.070499 -178.329423
35 82.268048 -92.070499
36 -20.136386 82.268048
37 -133.274337 -20.136386
38 -60.732183 -133.274337
39 457.194144 -60.732183
40 -116.233340 457.194144
41 108.180276 -116.233340
42 49.462690 108.180276
43 -93.909532 49.462690
44 -178.832054 -93.909532
45 -145.231416 -178.832054
46 -27.555206 -145.231416
47 -165.045612 -27.555206
48 -200.481235 -165.045612
49 -67.901140 -200.481235
50 13.701477 -67.901140
51 71.266824 13.701477
52 166.810267 71.266824
53 1.973891 166.810267
54 -4.331654 1.973891
55 -64.843546 -4.331654
56 -46.620756 -64.843546
57 119.303241 -46.620756
58 -116.745097 119.303241
59 -181.796224 -116.745097
60 -59.231809 -181.796224
61 100.513385 -59.231809
62 NA 100.513385
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -144.529850 230.101727
[2,] 62.283211 -144.529850
[3,] 106.218012 62.283211
[4,] 43.889999 106.218012
[5,] 57.042973 43.889999
[6,] -86.017209 57.042973
[7,] 53.031250 -86.017209
[8,] 330.153031 53.031250
[9,] 94.847249 330.153031
[10,] -162.837226 94.847249
[11,] -41.646879 -162.837226
[12,] 148.952666 -41.646879
[13,] 45.671068 148.952666
[14,] -26.706366 45.671068
[15,] -4.706363 -26.706366
[16,] -8.354698 -4.706363
[17,] 246.863157 -8.354698
[18,] -184.166559 246.863157
[19,] 117.346918 -184.166559
[20,] 94.023792 117.346918
[21,] 67.995603 94.023792
[22,] -123.253217 67.995603
[23,] -67.553918 -123.253217
[24,] 163.707043 -67.553918
[25,] -169.082144 163.707043
[26,] 4.227350 -169.082144
[27,] -99.808706 4.227350
[28,] -233.494975 -99.808706
[29,] -86.258770 -233.494975
[30,] 276.318628 -86.258770
[31,] -27.991864 276.318628
[32,] 106.362274 -27.991864
[33,] -178.329423 106.362274
[34,] -92.070499 -178.329423
[35,] 82.268048 -92.070499
[36,] -20.136386 82.268048
[37,] -133.274337 -20.136386
[38,] -60.732183 -133.274337
[39,] 457.194144 -60.732183
[40,] -116.233340 457.194144
[41,] 108.180276 -116.233340
[42,] 49.462690 108.180276
[43,] -93.909532 49.462690
[44,] -178.832054 -93.909532
[45,] -145.231416 -178.832054
[46,] -27.555206 -145.231416
[47,] -165.045612 -27.555206
[48,] -200.481235 -165.045612
[49,] -67.901140 -200.481235
[50,] 13.701477 -67.901140
[51,] 71.266824 13.701477
[52,] 166.810267 71.266824
[53,] 1.973891 166.810267
[54,] -4.331654 1.973891
[55,] -64.843546 -4.331654
[56,] -46.620756 -64.843546
[57,] 119.303241 -46.620756
[58,] -116.745097 119.303241
[59,] -181.796224 -116.745097
[60,] -59.231809 -181.796224
[61,] 100.513385 -59.231809
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -144.529850 230.101727
2 62.283211 -144.529850
3 106.218012 62.283211
4 43.889999 106.218012
5 57.042973 43.889999
6 -86.017209 57.042973
7 53.031250 -86.017209
8 330.153031 53.031250
9 94.847249 330.153031
10 -162.837226 94.847249
11 -41.646879 -162.837226
12 148.952666 -41.646879
13 45.671068 148.952666
14 -26.706366 45.671068
15 -4.706363 -26.706366
16 -8.354698 -4.706363
17 246.863157 -8.354698
18 -184.166559 246.863157
19 117.346918 -184.166559
20 94.023792 117.346918
21 67.995603 94.023792
22 -123.253217 67.995603
23 -67.553918 -123.253217
24 163.707043 -67.553918
25 -169.082144 163.707043
26 4.227350 -169.082144
27 -99.808706 4.227350
28 -233.494975 -99.808706
29 -86.258770 -233.494975
30 276.318628 -86.258770
31 -27.991864 276.318628
32 106.362274 -27.991864
33 -178.329423 106.362274
34 -92.070499 -178.329423
35 82.268048 -92.070499
36 -20.136386 82.268048
37 -133.274337 -20.136386
38 -60.732183 -133.274337
39 457.194144 -60.732183
40 -116.233340 457.194144
41 108.180276 -116.233340
42 49.462690 108.180276
43 -93.909532 49.462690
44 -178.832054 -93.909532
45 -145.231416 -178.832054
46 -27.555206 -145.231416
47 -165.045612 -27.555206
48 -200.481235 -165.045612
49 -67.901140 -200.481235
50 13.701477 -67.901140
51 71.266824 13.701477
52 166.810267 71.266824
53 1.973891 166.810267
54 -4.331654 1.973891
55 -64.843546 -4.331654
56 -46.620756 -64.843546
57 119.303241 -46.620756
58 -116.745097 119.303241
59 -181.796224 -116.745097
60 -59.231809 -181.796224
61 100.513385 -59.231809
> 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/fisher/rcomp/tmp/7zjms1354819135.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/fisher/rcomp/tmp/8k7ew1354819135.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/fisher/rcomp/tmp/9zzx61354819135.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/fisher/rcomp/tmp/10iqka1354819135.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/11p7su1354819135.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/fisher/rcomp/tmp/12ux0q1354819135.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/fisher/rcomp/tmp/134aak1354819135.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/fisher/rcomp/tmp/14rfko1354819135.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/fisher/rcomp/tmp/15kjz61354819135.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/fisher/rcomp/tmp/16uue01354819135.tab")
+ }
>
> try(system("convert tmp/1r06a1354819135.ps tmp/1r06a1354819135.png",intern=TRUE))
character(0)
> try(system("convert tmp/24yms1354819135.ps tmp/24yms1354819135.png",intern=TRUE))
character(0)
> try(system("convert tmp/3fv6m1354819135.ps tmp/3fv6m1354819135.png",intern=TRUE))
character(0)
> try(system("convert tmp/4rn951354819135.ps tmp/4rn951354819135.png",intern=TRUE))
character(0)
> try(system("convert tmp/5dyl71354819135.ps tmp/5dyl71354819135.png",intern=TRUE))
character(0)
> try(system("convert tmp/6zwjs1354819135.ps tmp/6zwjs1354819135.png",intern=TRUE))
character(0)
> try(system("convert tmp/7zjms1354819135.ps tmp/7zjms1354819135.png",intern=TRUE))
character(0)
> try(system("convert tmp/8k7ew1354819135.ps tmp/8k7ew1354819135.png",intern=TRUE))
character(0)
> try(system("convert tmp/9zzx61354819135.ps tmp/9zzx61354819135.png",intern=TRUE))
character(0)
> try(system("convert tmp/10iqka1354819135.ps tmp/10iqka1354819135.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.139 1.571 7.751