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
+ ,79.6
+ ,7.2
+ ,367
+ ,1148
+ ,600
+ ,0.55
+ ,1
+ ,8.5
+ ,443
+ ,18068
+ ,372
+ ,3.665
+ ,32.3
+ ,5.7
+ ,365
+ ,7729
+ ,142
+ ,2.351
+ ,45.1
+ ,7.3
+ ,614
+ ,100484
+ ,432
+ ,29.76
+ ,190.8
+ ,7.5
+ ,385
+ ,16728
+ ,290
+ ,3.294
+ ,31.8
+ ,5
+ ,286
+ ,14630
+ ,346
+ ,3.287
+ ,678.4
+ ,6.7
+ ,397
+ ,4008
+ ,328
+ ,0.666
+ ,340.8
+ ,6.2
+ ,764
+ ,38927
+ ,354
+ ,12.938
+ ,239.6
+ ,7.3
+ ,427
+ ,22322
+ ,266
+ ,6.478
+ ,111.9
+ ,5
+ ,153
+ ,3711
+ ,320
+ ,1.108
+ ,172.5
+ ,2.8
+ ,231
+ ,3136
+ ,197
+ ,1.007
+ ,12.2
+ ,6.1
+ ,524
+ ,50508
+ ,266
+ ,11.431
+ ,205.6
+ ,7.1
+ ,328
+ ,28886
+ ,173
+ ,5.544
+ ,154.6
+ ,5.9
+ ,240
+ ,16996
+ ,190
+ ,2.777
+ ,49.7
+ ,4.6
+ ,286
+ ,13035
+ ,239
+ ,2.478
+ ,30.3
+ ,4.4
+ ,285
+ ,12973
+ ,190
+ ,3.685
+ ,92.8
+ ,7.4
+ ,569
+ ,16309
+ ,241
+ ,4.22
+ ,96.9
+ ,7.1
+ ,96
+ ,5227
+ ,189
+ ,1.228
+ ,39.8
+ ,7.5
+ ,498
+ ,19235
+ ,358
+ ,4.781
+ ,489.2
+ ,5.9
+ ,481
+ ,44487
+ ,315
+ ,6.016
+ ,767.6
+ ,9
+ ,468
+ ,44213
+ ,303
+ ,9.295
+ ,163.6
+ ,9.2
+ ,177
+ ,23619
+ ,228
+ ,4.375
+ ,55
+ ,5.1
+ ,198
+ ,9106
+ ,134
+ ,2.573
+ ,54.9
+ ,8.6
+ ,458
+ ,24917
+ ,189
+ ,5.117
+ ,74.3
+ ,6.6
+ ,108
+ ,3872
+ ,196
+ ,0.799
+ ,5.5
+ ,6.9
+ ,246
+ ,8945
+ ,183
+ ,1.578
+ ,20.5
+ ,2.7
+ ,291
+ ,2373
+ ,417
+ ,1.202
+ ,10.9
+ ,5.5
+ ,68
+ ,7128
+ ,233
+ ,1.109
+ ,123.7
+ ,7.2
+ ,311
+ ,23624
+ ,349
+ ,7.73
+ ,1042
+ ,6.6
+ ,606
+ ,5242
+ ,284
+ ,1.515
+ ,12.5
+ ,6.9
+ ,512
+ ,92629
+ ,499
+ ,17.99
+ ,381
+ ,7.2
+ ,426
+ ,28795
+ ,231
+ ,6.629
+ ,136.1
+ ,5.8
+ ,47
+ ,4487
+ ,143
+ ,0.639
+ ,9.3
+ ,4.1
+ ,265
+ ,48799
+ ,249
+ ,10.847
+ ,264.9
+ ,6.4
+ ,370
+ ,14067
+ ,195
+ ,3.146
+ ,45.8
+ ,6.7
+ ,312
+ ,12693
+ ,288
+ ,2.842
+ ,29.6
+ ,6
+ ,222
+ ,62184
+ ,229
+ ,11.882
+ ,265.1
+ ,6.9
+ ,280
+ ,9153
+ ,287
+ ,1.003
+ ,960.3
+ ,8.5
+ ,759
+ ,14250
+ ,224
+ ,3.487
+ ,115.8
+ ,6.2
+ ,114
+ ,3680
+ ,161
+ ,0.696
+ ,9.2
+ ,3.4
+ ,419
+ ,18063
+ ,221
+ ,4.877
+ ,118.3
+ ,6.6
+ ,435
+ ,65112
+ ,237
+ ,16.987
+ ,64.9
+ ,6.6
+ ,186
+ ,11340
+ ,220
+ ,1.723
+ ,21
+ ,4.9
+ ,87
+ ,4553
+ ,185
+ ,0.563
+ ,60.8
+ ,6.4
+ ,188
+ ,28960
+ ,260
+ ,6.187
+ ,156.3
+ ,5.8
+ ,303
+ ,19201
+ ,261
+ ,4.867
+ ,73.1
+ ,6.3
+ ,102
+ ,7533
+ ,118
+ ,1.793
+ ,74.5
+ ,10.5
+ ,127
+ ,26343
+ ,268
+ ,4.892
+ ,90.1
+ ,5.4
+ ,251
+ ,1641
+ ,300
+ ,0.454
+ ,4.7
+ ,5.1
+ ,205
+ ,145360
+ ,237
+ ,10.379
+ ,889
+ ,6.8
+ ,453
+ ,9066420
+ ,240
+ ,82.422
+ ,609
+ ,5.6
+ ,320
+ ,1038933
+ ,185
+ ,16.491
+ ,1259
+ ,3.8
+ ,405
+ ,2739420
+ ,201
+ ,60.876
+ ,289
+ ,8.2
+ ,89
+ ,61620
+ ,193
+ ,0.474
+ ,475
+ ,4.1
+ ,74
+ ,827530
+ ,254
+ ,7.523
+ ,490
+ ,2.8
+ ,101
+ ,534100
+ ,230
+ ,5.45
+ ,333
+ ,6.3
+ ,321
+ ,328755
+ ,197
+ ,10.605
+ ,300
+ ,11.4
+ ,315
+ ,1413895
+ ,248
+ ,40.397
+ ,210
+ ,19.4
+ ,229
+ ,2909136
+ ,258
+ ,60.607
+ ,650
+ ,5.8
+ ,302
+ ,3604246
+ ,206
+ ,58.133
+ ,512
+ ,6.9
+ ,216
+ ,917504
+ ,199
+ ,8.192
+ ,256
+ ,3.5)
+ ,dim=c(6
+ ,62)
+ ,dimnames=list(c('Assaults'
+ ,'BachDegrees'
+ ,'PoliceExp'
+ ,'Population'
+ ,'Density'
+ ,'Unemployment')
+ ,1:62))
> y <- array(NA,dim=c(6,62),dimnames=list(c('Assaults','BachDegrees','PoliceExp','Population','Density','Unemployment'),1:62))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = ''
> par9 = ''
> par8 = ''
> par7 = ''
> par6 = ''
> par5 = ''
> par4 = ''
> 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 Population Density Unemployment
1 521 18308 185 4.041 79.6 7.2
2 367 1148 600 0.550 1.0 8.5
3 443 18068 372 3.665 32.3 5.7
4 365 7729 142 2.351 45.1 7.3
5 614 100484 432 29.760 190.8 7.5
6 385 16728 290 3.294 31.8 5.0
7 286 14630 346 3.287 678.4 6.7
8 397 4008 328 0.666 340.8 6.2
9 764 38927 354 12.938 239.6 7.3
10 427 22322 266 6.478 111.9 5.0
11 153 3711 320 1.108 172.5 2.8
12 231 3136 197 1.007 12.2 6.1
13 524 50508 266 11.431 205.6 7.1
14 328 28886 173 5.544 154.6 5.9
15 240 16996 190 2.777 49.7 4.6
16 286 13035 239 2.478 30.3 4.4
17 285 12973 190 3.685 92.8 7.4
18 569 16309 241 4.220 96.9 7.1
19 96 5227 189 1.228 39.8 7.5
20 498 19235 358 4.781 489.2 5.9
21 481 44487 315 6.016 767.6 9.0
22 468 44213 303 9.295 163.6 9.2
23 177 23619 228 4.375 55.0 5.1
24 198 9106 134 2.573 54.9 8.6
25 458 24917 189 5.117 74.3 6.6
26 108 3872 196 0.799 5.5 6.9
27 246 8945 183 1.578 20.5 2.7
28 291 2373 417 1.202 10.9 5.5
29 68 7128 233 1.109 123.7 7.2
30 311 23624 349 7.730 1042.0 6.6
31 606 5242 284 1.515 12.5 6.9
32 512 92629 499 17.990 381.0 7.2
33 426 28795 231 6.629 136.1 5.8
34 47 4487 143 0.639 9.3 4.1
35 265 48799 249 10.847 264.9 6.4
36 370 14067 195 3.146 45.8 6.7
37 312 12693 288 2.842 29.6 6.0
38 222 62184 229 11.882 265.1 6.9
39 280 9153 287 1.003 960.3 8.5
40 759 14250 224 3.487 115.8 6.2
41 114 3680 161 0.696 9.2 3.4
42 419 18063 221 4.877 118.3 6.6
43 435 65112 237 16.987 64.9 6.6
44 186 11340 220 1.723 21.0 4.9
45 87 4553 185 0.563 60.8 6.4
46 188 28960 260 6.187 156.3 5.8
47 303 19201 261 4.867 73.1 6.3
48 102 7533 118 1.793 74.5 10.5
49 127 26343 268 4.892 90.1 5.4
50 251 1641 300 0.454 4.7 5.1
51 205 145360 237 10.379 889.0 6.8
52 453 9066420 240 82.422 609.0 5.6
53 320 1038933 185 16.491 1259.0 3.8
54 405 2739420 201 60.876 289.0 8.2
55 89 61620 193 0.474 475.0 4.1
56 74 827530 254 7.523 490.0 2.8
57 101 534100 230 5.450 333.0 6.3
58 321 328755 197 10.605 300.0 11.4
59 315 1413895 248 40.397 210.0 19.4
60 229 2909136 258 60.607 650.0 5.8
61 302 3604246 206 58.133 512.0 6.9
62 216 917504 199 8.192 256.0 3.5
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) BachDegrees PoliceExp Population Density
7.522e+01 -2.917e-05 7.091e-01 3.988e+00 -5.059e-02
Unemployment
7.031e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-221.01 -119.49 -17.66 96.13 473.73
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.522e+01 8.148e+01 0.923 0.3599
BachDegrees -2.917e-05 3.372e-05 -0.865 0.3908
PoliceExp 7.091e-01 2.341e-01 3.029 0.0037 **
Population 3.988e+00 2.788e+00 1.430 0.1582
Density -5.059e-02 7.411e-02 -0.683 0.4976
Unemployment 7.031e+00 9.039e+00 0.778 0.4400
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 153.4 on 56 degrees of freedom
Multiple R-squared: 0.2211, Adjusted R-squared: 0.1515
F-statistic: 3.179 on 5 and 56 DF, p-value: 0.0135
> 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.5469562 0.90608759 0.45304380
[2,] 0.4586445 0.91728892 0.54135554
[3,] 0.4900150 0.98002999 0.50998500
[4,] 0.5161491 0.96770180 0.48385090
[5,] 0.4181489 0.83629773 0.58185114
[6,] 0.3345650 0.66913003 0.66543499
[7,] 0.2524743 0.50494861 0.74752570
[8,] 0.1728362 0.34567241 0.82716379
[9,] 0.1723690 0.34473805 0.82763097
[10,] 0.2232072 0.44641446 0.77679277
[11,] 0.4201366 0.84027315 0.57986342
[12,] 0.3829030 0.76580596 0.61709702
[13,] 0.3234033 0.64680650 0.67659675
[14,] 0.2649142 0.52982850 0.73508575
[15,] 0.2484453 0.49689063 0.75155468
[16,] 0.2394210 0.47884206 0.76057897
[17,] 0.2513756 0.50275118 0.74862441
[18,] 0.2705722 0.54114439 0.72942781
[19,] 0.2088703 0.41774067 0.79112967
[20,] 0.1697475 0.33949500 0.83025250
[21,] 0.2489350 0.49787008 0.75106496
[22,] 0.2180240 0.43604803 0.78197598
[23,] 0.4209374 0.84187481 0.57906259
[24,] 0.3685981 0.73719617 0.63140191
[25,] 0.3516229 0.70324574 0.64837713
[26,] 0.3583963 0.71679251 0.64160374
[27,] 0.3186221 0.63724418 0.68137791
[28,] 0.2834917 0.56698330 0.71650835
[29,] 0.2202074 0.44041490 0.77979255
[30,] 0.1969094 0.39381885 0.80309058
[31,] 0.1474336 0.29486729 0.85256636
[32,] 0.9197309 0.16053823 0.08026911
[33,] 0.8887268 0.22254647 0.11127323
[34,] 0.9399162 0.12016760 0.06008380
[35,] 0.9795246 0.04095070 0.02047535
[36,] 0.9658146 0.06837078 0.03418539
[37,] 0.9548424 0.09031524 0.04515762
[38,] 0.9274115 0.14517708 0.07258854
[39,] 0.9339043 0.13219136 0.06609568
[40,] 0.9368565 0.12628706 0.06314353
[41,] 0.8958373 0.20832546 0.10416273
[42,] 0.9861404 0.02771920 0.01385960
[43,] 0.9686313 0.06273748 0.03136874
[44,] 0.9324254 0.13514914 0.06757457
[45,] 0.9195535 0.16089297 0.08044648
> postscript(file="/var/wessaorg/rcomp/tmp/195va1355044117.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/2d3e01355044117.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/3tccl1355044117.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/4v0bg1355044117.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/54dlp1355044117.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
252.428937 -195.537238 51.475661 130.901402 73.642456 57.957033
7 8 9 10 11 12
-60.023123 60.316853 348.109687 108.494318 -164.392118 -30.100510
13 14 15 16 17 18
176.542078 75.187650 -10.347293 2.408989 13.409177 261.526524
19 20 21 22 23 24
-168.694742 133.696959 135.289531 85.750360 -109.718659 -39.915143
25 26 27 28 29 30
186.445032 -157.504108 17.043619 -122.746045 -221.010390 -35.505827
31 32 33 34 35 36
275.635553 -17.429214 127.498710 -160.388416 -60.202867 99.588576
37 38 39 40 41 42
-19.082150 -96.263171 -13.632201 473.727671 -101.486295 127.737592
43 44 45 46 47 48
82.773265 -85.143131 -163.430043 -128.275389 -16.730305 -133.872129
49 50 51 52 53 54
-190.397497 -74.323304 -78.247723 134.833950 115.128850 -18.615725
55 56 57 58 59 60
-127.956538 -182.079635 -170.905169 8.422598 -181.685941 -193.872246
61 62
-68.574936 -17.883810
> postscript(file="/var/wessaorg/rcomp/tmp/6tcd61355044117.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 252.428937 NA
1 -195.537238 252.428937
2 51.475661 -195.537238
3 130.901402 51.475661
4 73.642456 130.901402
5 57.957033 73.642456
6 -60.023123 57.957033
7 60.316853 -60.023123
8 348.109687 60.316853
9 108.494318 348.109687
10 -164.392118 108.494318
11 -30.100510 -164.392118
12 176.542078 -30.100510
13 75.187650 176.542078
14 -10.347293 75.187650
15 2.408989 -10.347293
16 13.409177 2.408989
17 261.526524 13.409177
18 -168.694742 261.526524
19 133.696959 -168.694742
20 135.289531 133.696959
21 85.750360 135.289531
22 -109.718659 85.750360
23 -39.915143 -109.718659
24 186.445032 -39.915143
25 -157.504108 186.445032
26 17.043619 -157.504108
27 -122.746045 17.043619
28 -221.010390 -122.746045
29 -35.505827 -221.010390
30 275.635553 -35.505827
31 -17.429214 275.635553
32 127.498710 -17.429214
33 -160.388416 127.498710
34 -60.202867 -160.388416
35 99.588576 -60.202867
36 -19.082150 99.588576
37 -96.263171 -19.082150
38 -13.632201 -96.263171
39 473.727671 -13.632201
40 -101.486295 473.727671
41 127.737592 -101.486295
42 82.773265 127.737592
43 -85.143131 82.773265
44 -163.430043 -85.143131
45 -128.275389 -163.430043
46 -16.730305 -128.275389
47 -133.872129 -16.730305
48 -190.397497 -133.872129
49 -74.323304 -190.397497
50 -78.247723 -74.323304
51 134.833950 -78.247723
52 115.128850 134.833950
53 -18.615725 115.128850
54 -127.956538 -18.615725
55 -182.079635 -127.956538
56 -170.905169 -182.079635
57 8.422598 -170.905169
58 -181.685941 8.422598
59 -193.872246 -181.685941
60 -68.574936 -193.872246
61 -17.883810 -68.574936
62 NA -17.883810
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -195.537238 252.428937
[2,] 51.475661 -195.537238
[3,] 130.901402 51.475661
[4,] 73.642456 130.901402
[5,] 57.957033 73.642456
[6,] -60.023123 57.957033
[7,] 60.316853 -60.023123
[8,] 348.109687 60.316853
[9,] 108.494318 348.109687
[10,] -164.392118 108.494318
[11,] -30.100510 -164.392118
[12,] 176.542078 -30.100510
[13,] 75.187650 176.542078
[14,] -10.347293 75.187650
[15,] 2.408989 -10.347293
[16,] 13.409177 2.408989
[17,] 261.526524 13.409177
[18,] -168.694742 261.526524
[19,] 133.696959 -168.694742
[20,] 135.289531 133.696959
[21,] 85.750360 135.289531
[22,] -109.718659 85.750360
[23,] -39.915143 -109.718659
[24,] 186.445032 -39.915143
[25,] -157.504108 186.445032
[26,] 17.043619 -157.504108
[27,] -122.746045 17.043619
[28,] -221.010390 -122.746045
[29,] -35.505827 -221.010390
[30,] 275.635553 -35.505827
[31,] -17.429214 275.635553
[32,] 127.498710 -17.429214
[33,] -160.388416 127.498710
[34,] -60.202867 -160.388416
[35,] 99.588576 -60.202867
[36,] -19.082150 99.588576
[37,] -96.263171 -19.082150
[38,] -13.632201 -96.263171
[39,] 473.727671 -13.632201
[40,] -101.486295 473.727671
[41,] 127.737592 -101.486295
[42,] 82.773265 127.737592
[43,] -85.143131 82.773265
[44,] -163.430043 -85.143131
[45,] -128.275389 -163.430043
[46,] -16.730305 -128.275389
[47,] -133.872129 -16.730305
[48,] -190.397497 -133.872129
[49,] -74.323304 -190.397497
[50,] -78.247723 -74.323304
[51,] 134.833950 -78.247723
[52,] 115.128850 134.833950
[53,] -18.615725 115.128850
[54,] -127.956538 -18.615725
[55,] -182.079635 -127.956538
[56,] -170.905169 -182.079635
[57,] 8.422598 -170.905169
[58,] -181.685941 8.422598
[59,] -193.872246 -181.685941
[60,] -68.574936 -193.872246
[61,] -17.883810 -68.574936
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -195.537238 252.428937
2 51.475661 -195.537238
3 130.901402 51.475661
4 73.642456 130.901402
5 57.957033 73.642456
6 -60.023123 57.957033
7 60.316853 -60.023123
8 348.109687 60.316853
9 108.494318 348.109687
10 -164.392118 108.494318
11 -30.100510 -164.392118
12 176.542078 -30.100510
13 75.187650 176.542078
14 -10.347293 75.187650
15 2.408989 -10.347293
16 13.409177 2.408989
17 261.526524 13.409177
18 -168.694742 261.526524
19 133.696959 -168.694742
20 135.289531 133.696959
21 85.750360 135.289531
22 -109.718659 85.750360
23 -39.915143 -109.718659
24 186.445032 -39.915143
25 -157.504108 186.445032
26 17.043619 -157.504108
27 -122.746045 17.043619
28 -221.010390 -122.746045
29 -35.505827 -221.010390
30 275.635553 -35.505827
31 -17.429214 275.635553
32 127.498710 -17.429214
33 -160.388416 127.498710
34 -60.202867 -160.388416
35 99.588576 -60.202867
36 -19.082150 99.588576
37 -96.263171 -19.082150
38 -13.632201 -96.263171
39 473.727671 -13.632201
40 -101.486295 473.727671
41 127.737592 -101.486295
42 82.773265 127.737592
43 -85.143131 82.773265
44 -163.430043 -85.143131
45 -128.275389 -163.430043
46 -16.730305 -128.275389
47 -133.872129 -16.730305
48 -190.397497 -133.872129
49 -74.323304 -190.397497
50 -78.247723 -74.323304
51 134.833950 -78.247723
52 115.128850 134.833950
53 -18.615725 115.128850
54 -127.956538 -18.615725
55 -182.079635 -127.956538
56 -170.905169 -182.079635
57 8.422598 -170.905169
58 -181.685941 8.422598
59 -193.872246 -181.685941
60 -68.574936 -193.872246
61 -17.883810 -68.574936
> 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/75z1x1355044117.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/89a4t1355044117.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/9ikl51355044117.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/10ov5p1355044117.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/110zph1355044117.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/12dcme1355044117.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/136qq11355044117.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/14q8yx1355044117.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/15c9311355044117.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/16cyg91355044117.tab")
+ }
>
> try(system("convert tmp/195va1355044117.ps tmp/195va1355044117.png",intern=TRUE))
character(0)
> try(system("convert tmp/2d3e01355044117.ps tmp/2d3e01355044117.png",intern=TRUE))
character(0)
> try(system("convert tmp/3tccl1355044117.ps tmp/3tccl1355044117.png",intern=TRUE))
character(0)
> try(system("convert tmp/4v0bg1355044117.ps tmp/4v0bg1355044117.png",intern=TRUE))
character(0)
> try(system("convert tmp/54dlp1355044117.ps tmp/54dlp1355044117.png",intern=TRUE))
character(0)
> try(system("convert tmp/6tcd61355044117.ps tmp/6tcd61355044117.png",intern=TRUE))
character(0)
> try(system("convert tmp/75z1x1355044117.ps tmp/75z1x1355044117.png",intern=TRUE))
character(0)
> try(system("convert tmp/89a4t1355044117.ps tmp/89a4t1355044117.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ikl51355044117.ps tmp/9ikl51355044117.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ov5p1355044117.ps tmp/10ov5p1355044117.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.435 1.157 7.639