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(18897
+ ,22424
+ ,19364
+ ,19434
+ ,22831
+ ,23072
+ ,37471
+ ,14690
+ ,17518
+ ,22125
+ ,18586
+ ,18389
+ ,22727
+ ,22551
+ ,36160
+ ,13824
+ ,8632
+ ,7653
+ ,8225
+ ,8405
+ ,8344
+ ,8695
+ ,9197
+ ,9477
+ ,832
+ ,554
+ ,822
+ ,854
+ ,830
+ ,935
+ ,1051
+ ,1150
+ ,3351
+ ,3357
+ ,3270
+ ,3346
+ ,3235
+ ,3329
+ ,3480
+ ,3447
+ ,8
+ ,8
+ ,3
+ ,4
+ ,5
+ ,5
+ ,4
+ ,4
+ ,1
+ ,1
+ ,1
+ ,1
+ ,1
+ ,1
+ ,1
+ ,2
+ ,7
+ ,10
+ ,11
+ ,9
+ ,10
+ ,9
+ ,10
+ ,9
+ ,217
+ ,222
+ ,204
+ ,205
+ ,191
+ ,197
+ ,196
+ ,191
+ ,911
+ ,947
+ ,918
+ ,939
+ ,937
+ ,967
+ ,1007
+ ,962
+ ,1932
+ ,1901
+ ,1862
+ ,1921
+ ,1823
+ ,1879
+ ,1982
+ ,2003
+ ,274
+ ,267
+ ,270
+ ,267
+ ,269
+ ,271
+ ,281
+ ,276
+ ,131
+ ,109
+ ,87
+ ,66
+ ,68
+ ,64
+ ,76
+ ,81
+ ,1708
+ ,1668
+ ,1738
+ ,1715
+ ,1726
+ ,1771
+ ,1861
+ ,2079
+ ,2609
+ ,1965
+ ,2308
+ ,2424
+ ,2486
+ ,2594
+ ,2729
+ ,2720
+ ,133
+ ,32
+ ,119
+ ,89
+ ,93
+ ,107
+ ,102
+ ,23
+ ,2476
+ ,1933
+ ,2189
+ ,2335
+ ,2393
+ ,2487
+ ,2627
+ ,2697
+ ,10
+ ,37
+ ,23
+ ,21
+ ,22
+ ,27
+ ,21
+ ,18
+ ,1510
+ ,1616
+ ,1378
+ ,1605
+ ,1534
+ ,1654
+ ,1421
+ ,1650
+ ,6427
+ ,7719
+ ,8279
+ ,6133
+ ,11706
+ ,9235
+ ,24339
+ ,1324
+ ,3812
+ ,5127
+ ,5890
+ ,4487
+ ,7888
+ ,6772
+ ,9522
+ ,3656
+ ,724
+ ,99
+ ,154
+ ,157
+ ,367
+ ,153
+ ,171
+ ,-211
+ ,1560
+ ,1996
+ ,1917
+ ,1223
+ ,2860
+ ,1964
+ ,13508
+ ,-2076
+ ,156
+ ,113
+ ,166
+ ,116
+ ,435
+ ,166
+ ,975
+ ,-178
+ ,3
+ ,3
+ ,2
+ ,2
+ ,15
+ ,13
+ ,27
+ ,13
+ ,172
+ ,380
+ ,151
+ ,148
+ ,141
+ ,167
+ ,137
+ ,120
+ ,65
+ ,1700
+ ,163
+ ,568
+ ,80
+ ,2043
+ ,768
+ ,338
+ ,593
+ ,2931
+ ,292
+ ,1348
+ ,774
+ ,631
+ ,122
+ ,698
+ ,281
+ ,469
+ ,227
+ ,309
+ ,266
+ ,267
+ ,292
+ ,321
+ ,1191
+ ,145
+ ,747
+ ,874
+ ,38
+ ,275
+ ,423
+ ,218
+ ,72
+ ,57
+ ,2
+ ,32
+ ,32
+ ,57
+ ,16
+ ,21
+ ,113
+ ,-6
+ ,27
+ ,120
+ ,32
+ ,184
+ ,860
+ ,604
+ ,19
+ ,86
+ ,2
+ ,18
+ ,3
+ ,3
+ ,12
+ ,246
+ ,97
+ ,11
+ ,27
+ ,120
+ ,32
+ ,185
+ ,860
+ ,381
+ ,18897
+ ,22424
+ ,19364
+ ,19434
+ ,22831
+ ,23072
+ ,37471
+ ,14690
+ ,16770
+ ,18775
+ ,17704
+ ,16289
+ ,21687
+ ,20252
+ ,35933
+ ,13873
+ ,6132
+ ,5145
+ ,5705
+ ,5818
+ ,5817
+ ,6171
+ ,6504
+ ,6749
+ ,648
+ ,299
+ ,484
+ ,535
+ ,511
+ ,548
+ ,638
+ ,758
+ ,1739
+ ,1710
+ ,1776
+ ,1910
+ ,1843
+ ,1990
+ ,2141
+ ,2097
+ ,160
+ ,167
+ ,176
+ ,193
+ ,183
+ ,202
+ ,163
+ ,179
+ ,621
+ ,570
+ ,592
+ ,743
+ ,655
+ ,735
+ ,851
+ ,835
+ ,804
+ ,821
+ ,842
+ ,831
+ ,847
+ ,869
+ ,886
+ ,881
+ ,3
+ ,3
+ ,3
+ ,3
+ ,3
+ ,3
+ ,3
+ ,3
+ ,150
+ ,149
+ ,163
+ ,140
+ ,156
+ ,182
+ ,238
+ ,199
+ ,549
+ ,528
+ ,558
+ ,354
+ ,470
+ ,438
+ ,450
+ ,453
+ ,95
+ ,80
+ ,132
+ ,-46
+ ,88
+ ,49
+ ,57
+ ,62
+ ,354
+ ,353
+ ,339
+ ,323
+ ,308
+ ,314
+ ,325
+ ,322
+ ,100
+ ,95
+ ,87
+ ,77
+ ,73
+ ,75
+ ,68
+ ,70
+ ,342
+ ,343
+ ,357
+ ,354
+ ,339
+ ,343
+ ,343
+ ,305
+ ,2854
+ ,2265
+ ,2530
+ ,2664
+ ,2653
+ ,2852
+ ,2932
+ ,3135
+ ,167
+ ,99
+ ,168
+ ,132
+ ,137
+ ,147
+ ,132
+ ,35
+ ,2687
+ ,2166
+ ,2362
+ ,2533
+ ,2516
+ ,2705
+ ,2799
+ ,3100
+ ,645
+ ,770
+ ,634
+ ,680
+ ,581
+ ,675
+ ,814
+ ,677
+ ,6113
+ ,7729
+ ,8065
+ ,5931
+ ,11602
+ ,9279
+ ,24726
+ ,1473
+ ,3567
+ ,4697
+ ,5792
+ ,4959
+ ,8473
+ ,6753
+ ,9199
+ ,3926
+ ,472
+ ,241
+ ,87
+ ,262
+ ,330
+ ,64
+ ,172
+ ,-142
+ ,1665
+ ,2360
+ ,1934
+ ,584
+ ,2229
+ ,1972
+ ,13856
+ ,-2338
+ ,328
+ ,318
+ ,154
+ ,6
+ ,256
+ ,181
+ ,1301
+ ,-241
+ ,0
+ ,1
+ ,0
+ ,0
+ ,12
+ ,10
+ ,21
+ ,10
+ ,81
+ ,112
+ ,99
+ ,120
+ ,302
+ ,300
+ ,177
+ ,259
+ ,1322
+ ,1286
+ ,1317
+ ,1325
+ ,1314
+ ,1322
+ ,1308
+ ,1370
+ ,154
+ ,143
+ ,156
+ ,152
+ ,144
+ ,151
+ ,151
+ ,171
+ ,1277
+ ,1448
+ ,1340
+ ,1689
+ ,1529
+ ,1544
+ ,1264
+ ,1656
+ ,1127
+ ,2253
+ ,486
+ ,694
+ ,699
+ ,1110
+ ,1165
+ ,1776
+ ,456
+ ,1356
+ ,63
+ ,2861
+ ,89
+ ,82
+ ,1019
+ ,926
+ ,224
+ ,200
+ ,149
+ ,91
+ ,165
+ ,216
+ ,94
+ ,131
+ ,1444
+ ,1990
+ ,1445
+ ,176
+ ,888
+ ,2517
+ ,413
+ ,-264
+ ,3
+ ,8
+ ,4
+ ,7
+ ,40
+ ,38
+ ,-64
+ ,-10
+ ,1444
+ ,2084
+ ,1443
+ ,187
+ ,850
+ ,2483
+ ,489
+ ,-231)
+ ,dim=c(8
+ ,69)
+ ,dimnames=list(c('2010-I'
+ ,'2010-II'
+ ,'2010-III'
+ ,'2010-IV'
+ ,'2011-I'
+ ,'2011-II'
+ ,'2011-III'
+ ,'2011-IV')
+ ,1:69))
> y <- array(NA,dim=c(8,69),dimnames=list(c('2010-I','2010-II','2010-III','2010-IV','2011-I','2011-II','2011-III','2011-IV'),1:69))
> 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 = '7'
> 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
2011-III 2010-I 2010-II 2010-III 2010-IV 2011-I 2011-II 2011-IV
1 37471 18897 22424 19364 19434 22831 23072 14690
2 36160 17518 22125 18586 18389 22727 22551 13824
3 9197 8632 7653 8225 8405 8344 8695 9477
4 1051 832 554 822 854 830 935 1150
5 3480 3351 3357 3270 3346 3235 3329 3447
6 4 8 8 3 4 5 5 4
7 1 1 1 1 1 1 1 2
8 10 7 10 11 9 10 9 9
9 196 217 222 204 205 191 197 191
10 1007 911 947 918 939 937 967 962
11 1982 1932 1901 1862 1921 1823 1879 2003
12 281 274 267 270 267 269 271 276
13 76 131 109 87 66 68 64 81
14 1861 1708 1668 1738 1715 1726 1771 2079
15 2729 2609 1965 2308 2424 2486 2594 2720
16 102 133 32 119 89 93 107 23
17 2627 2476 1933 2189 2335 2393 2487 2697
18 21 10 37 23 21 22 27 18
19 1421 1510 1616 1378 1605 1534 1654 1650
20 24339 6427 7719 8279 6133 11706 9235 1324
21 9522 3812 5127 5890 4487 7888 6772 3656
22 171 724 99 154 157 367 153 -211
23 13508 1560 1996 1917 1223 2860 1964 -2076
24 975 156 113 166 116 435 166 -178
25 27 3 3 2 2 15 13 13
26 137 172 380 151 148 141 167 120
27 768 65 1700 163 568 80 2043 338
28 122 593 2931 292 1348 774 631 698
29 292 281 469 227 309 266 267 321
30 423 1191 145 747 874 38 275 218
31 16 72 57 2 32 32 57 21
32 860 113 -6 27 120 32 184 604
33 12 19 86 2 18 3 3 246
34 860 97 11 27 120 32 185 381
35 37471 18897 22424 19364 19434 22831 23072 14690
36 35933 16770 18775 17704 16289 21687 20252 13873
37 6504 6132 5145 5705 5818 5817 6171 6749
38 638 648 299 484 535 511 548 758
39 2141 1739 1710 1776 1910 1843 1990 2097
40 163 160 167 176 193 183 202 179
41 851 621 570 592 743 655 735 835
42 886 804 821 842 831 847 869 881
43 3 3 3 3 3 3 3 3
44 238 150 149 163 140 156 182 199
45 450 549 528 558 354 470 438 453
46 57 95 80 132 -46 88 49 62
47 325 354 353 339 323 308 314 322
48 68 100 95 87 77 73 75 70
49 343 342 343 357 354 339 343 305
50 2932 2854 2265 2530 2664 2653 2852 3135
51 132 167 99 168 132 137 147 35
52 2799 2687 2166 2362 2533 2516 2705 3100
53 814 645 770 634 680 581 675 677
54 24726 6113 7729 8065 5931 11602 9279 1473
55 9199 3567 4697 5792 4959 8473 6753 3926
56 172 472 241 87 262 330 64 -142
57 13856 1665 2360 1934 584 2229 1972 -2338
58 1301 328 318 154 6 256 181 -241
59 21 0 1 0 0 12 10 10
60 177 81 112 99 120 302 300 259
61 1308 1322 1286 1317 1325 1314 1322 1370
62 151 154 143 156 152 144 151 171
63 1264 1277 1448 1340 1689 1529 1544 1656
64 1165 1127 2253 486 694 699 1110 1776
65 1019 456 1356 63 2861 89 82 926
66 94 224 200 149 91 165 216 131
67 413 1444 1990 1445 176 888 2517 -264
68 -64 3 8 4 7 40 38 -10
69 489 1444 2084 1443 187 850 2483 -231
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `2010-I` `2010-II` `2010-III` `2010-IV` `2011-I`
-18.9779 2.4812 0.1358 -1.2037 0.4509 2.2482
`2011-II` `2011-IV`
-0.8466 -1.9952
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2621.8 -253.0 -47.5 200.6 3481.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -18.9779 141.2897 -0.134 0.893592
`2010-I` 2.4812 0.6032 4.114 0.000119 ***
`2010-II` 0.1358 0.3884 0.350 0.727915
`2010-III` -1.2037 1.1601 -1.038 0.303550
`2010-IV` 0.4509 0.3956 1.140 0.258885
`2011-I` 2.2482 0.4714 4.770 1.19e-05 ***
`2011-II` -0.8466 0.5681 -1.490 0.141332
`2011-IV` -1.9952 0.1902 -10.491 2.76e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1039 on 61 degrees of freedom
Multiple R-squared: 0.9891, Adjusted R-squared: 0.9879
F-statistic: 792.6 on 7 and 61 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,] 7.257962e-05 1.451592e-04 0.9999274204
[2,] 2.841364e-06 5.682728e-06 0.9999971586
[3,] 6.069166e-07 1.213833e-06 0.9999993931
[4,] 4.630459e-07 9.260918e-07 0.9999995370
[5,] 5.442379e-08 1.088476e-07 0.9999999456
[6,] 4.000362e-09 8.000724e-09 0.9999999960
[7,] 2.559524e-10 5.119047e-10 0.9999999997
[8,] 2.251841e-11 4.503681e-11 1.0000000000
[9,] 2.650626e-09 5.301252e-09 0.9999999973
[10,] 2.358313e-09 4.716626e-09 0.9999999976
[11,] 3.870017e-07 7.740034e-07 0.9999996130
[12,] 7.432918e-04 1.486584e-03 0.9992567082
[13,] 4.738069e-03 9.476138e-03 0.9952619311
[14,] 2.562785e-03 5.125570e-03 0.9974372152
[15,] 1.211507e-03 2.423014e-03 0.9987884929
[16,] 6.017080e-04 1.203416e-03 0.9993982920
[17,] 1.989666e-03 3.979332e-03 0.9980103340
[18,] 5.465485e-03 1.093097e-02 0.9945345147
[19,] 3.340300e-03 6.680600e-03 0.9966596998
[20,] 3.667171e-01 7.334341e-01 0.6332829329
[21,] 2.953001e-01 5.906002e-01 0.7046999020
[22,] 6.965330e-01 6.069339e-01 0.3034669671
[23,] 6.421548e-01 7.156903e-01 0.3578451671
[24,] 8.587508e-01 2.824984e-01 0.1412491885
[25,] 8.315182e-01 3.369636e-01 0.1684817906
[26,] 8.933202e-01 2.133597e-01 0.1066798354
[27,] 8.652856e-01 2.694288e-01 0.1347144001
[28,] 8.462044e-01 3.075911e-01 0.1537955600
[29,] 8.250158e-01 3.499684e-01 0.1749841875
[30,] 7.715035e-01 4.569930e-01 0.2284964897
[31,] 7.916860e-01 4.166281e-01 0.2083140283
[32,] 7.338487e-01 5.323026e-01 0.2661512971
[33,] 6.701290e-01 6.597420e-01 0.3298709806
[34,] 6.099681e-01 7.800638e-01 0.3900319075
[35,] 6.805323e-01 6.389354e-01 0.3194677008
[36,] 6.205352e-01 7.589295e-01 0.3794647508
[37,] 5.493099e-01 9.013803e-01 0.4506901450
[38,] 4.590169e-01 9.180339e-01 0.5409830555
[39,] 3.792307e-01 7.584614e-01 0.6207693176
[40,] 3.190959e-01 6.381918e-01 0.6809041201
[41,] 2.418627e-01 4.837253e-01 0.7581373343
[42,] 8.486553e-01 3.026893e-01 0.1513446606
[43,] 7.748394e-01 4.503213e-01 0.2251606428
[44,] 7.746897e-01 4.506207e-01 0.2253103334
[45,] 9.984944e-01 3.011267e-03 0.0015056334
[46,] 9.993032e-01 1.393670e-03 0.0006968350
[47,] 9.995481e-01 9.038924e-04 0.0004519462
[48,] 9.978952e-01 4.209679e-03 0.0021048397
> postscript(file="/var/wessaorg/rcomp/tmp/1v1d21355151226.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/2y1wm1355151226.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/3poxy1355151226.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/4e2x81355151226.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/5f6ev1355151226.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 = 69
Frequency = 1
1 2 3 4 5 6
-381.996700 -631.297216 380.980558 754.872113 -420.894473 4.822662
7 8 9 10 11 12
20.702637 22.528968 -81.996508 -49.843028 -186.914513 -36.162178
13 14 15 16 17 18
-106.975327 501.411993 -272.982300 -182.862680 -71.141708 37.670686
19 20 21 22 23 24
-347.338587 -1294.232626 -252.996168 -2621.831784 2232.026863 -453.516859
25 26 27 28 29 30
42.853132 -143.538792 2559.117391 -1797.832710 -47.505455 -1445.364329
31 32 33 34 35 36
-145.214467 1866.772236 453.077688 1460.069331 -381.996700 1827.234775
37 38 39 40 41 42
465.995968 177.428640 614.888120 3.881202 445.047800 26.723747
43 44 45 46 47 48
18.166340 198.060230 -234.841449 -23.625158 -104.023297 -64.993694
49 50 51 52 53 54
-126.265109 111.401259 -247.841400 376.769705 200.630312 272.346197
55 56 57 58 59 60
-1033.076622 -1997.327925 3481.298222 -257.531783 41.281764 136.637626
61 62 63 64 65 66
-241.380807 32.981457 -56.937370 1265.651544 225.289966 -258.358197
67 68 69
-2153.370994 -129.604427 -1975.007793
> postscript(file="/var/wessaorg/rcomp/tmp/6344p1355151226.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 -381.996700 NA
1 -631.297216 -381.996700
2 380.980558 -631.297216
3 754.872113 380.980558
4 -420.894473 754.872113
5 4.822662 -420.894473
6 20.702637 4.822662
7 22.528968 20.702637
8 -81.996508 22.528968
9 -49.843028 -81.996508
10 -186.914513 -49.843028
11 -36.162178 -186.914513
12 -106.975327 -36.162178
13 501.411993 -106.975327
14 -272.982300 501.411993
15 -182.862680 -272.982300
16 -71.141708 -182.862680
17 37.670686 -71.141708
18 -347.338587 37.670686
19 -1294.232626 -347.338587
20 -252.996168 -1294.232626
21 -2621.831784 -252.996168
22 2232.026863 -2621.831784
23 -453.516859 2232.026863
24 42.853132 -453.516859
25 -143.538792 42.853132
26 2559.117391 -143.538792
27 -1797.832710 2559.117391
28 -47.505455 -1797.832710
29 -1445.364329 -47.505455
30 -145.214467 -1445.364329
31 1866.772236 -145.214467
32 453.077688 1866.772236
33 1460.069331 453.077688
34 -381.996700 1460.069331
35 1827.234775 -381.996700
36 465.995968 1827.234775
37 177.428640 465.995968
38 614.888120 177.428640
39 3.881202 614.888120
40 445.047800 3.881202
41 26.723747 445.047800
42 18.166340 26.723747
43 198.060230 18.166340
44 -234.841449 198.060230
45 -23.625158 -234.841449
46 -104.023297 -23.625158
47 -64.993694 -104.023297
48 -126.265109 -64.993694
49 111.401259 -126.265109
50 -247.841400 111.401259
51 376.769705 -247.841400
52 200.630312 376.769705
53 272.346197 200.630312
54 -1033.076622 272.346197
55 -1997.327925 -1033.076622
56 3481.298222 -1997.327925
57 -257.531783 3481.298222
58 41.281764 -257.531783
59 136.637626 41.281764
60 -241.380807 136.637626
61 32.981457 -241.380807
62 -56.937370 32.981457
63 1265.651544 -56.937370
64 225.289966 1265.651544
65 -258.358197 225.289966
66 -2153.370994 -258.358197
67 -129.604427 -2153.370994
68 -1975.007793 -129.604427
69 NA -1975.007793
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -631.297216 -381.996700
[2,] 380.980558 -631.297216
[3,] 754.872113 380.980558
[4,] -420.894473 754.872113
[5,] 4.822662 -420.894473
[6,] 20.702637 4.822662
[7,] 22.528968 20.702637
[8,] -81.996508 22.528968
[9,] -49.843028 -81.996508
[10,] -186.914513 -49.843028
[11,] -36.162178 -186.914513
[12,] -106.975327 -36.162178
[13,] 501.411993 -106.975327
[14,] -272.982300 501.411993
[15,] -182.862680 -272.982300
[16,] -71.141708 -182.862680
[17,] 37.670686 -71.141708
[18,] -347.338587 37.670686
[19,] -1294.232626 -347.338587
[20,] -252.996168 -1294.232626
[21,] -2621.831784 -252.996168
[22,] 2232.026863 -2621.831784
[23,] -453.516859 2232.026863
[24,] 42.853132 -453.516859
[25,] -143.538792 42.853132
[26,] 2559.117391 -143.538792
[27,] -1797.832710 2559.117391
[28,] -47.505455 -1797.832710
[29,] -1445.364329 -47.505455
[30,] -145.214467 -1445.364329
[31,] 1866.772236 -145.214467
[32,] 453.077688 1866.772236
[33,] 1460.069331 453.077688
[34,] -381.996700 1460.069331
[35,] 1827.234775 -381.996700
[36,] 465.995968 1827.234775
[37,] 177.428640 465.995968
[38,] 614.888120 177.428640
[39,] 3.881202 614.888120
[40,] 445.047800 3.881202
[41,] 26.723747 445.047800
[42,] 18.166340 26.723747
[43,] 198.060230 18.166340
[44,] -234.841449 198.060230
[45,] -23.625158 -234.841449
[46,] -104.023297 -23.625158
[47,] -64.993694 -104.023297
[48,] -126.265109 -64.993694
[49,] 111.401259 -126.265109
[50,] -247.841400 111.401259
[51,] 376.769705 -247.841400
[52,] 200.630312 376.769705
[53,] 272.346197 200.630312
[54,] -1033.076622 272.346197
[55,] -1997.327925 -1033.076622
[56,] 3481.298222 -1997.327925
[57,] -257.531783 3481.298222
[58,] 41.281764 -257.531783
[59,] 136.637626 41.281764
[60,] -241.380807 136.637626
[61,] 32.981457 -241.380807
[62,] -56.937370 32.981457
[63,] 1265.651544 -56.937370
[64,] 225.289966 1265.651544
[65,] -258.358197 225.289966
[66,] -2153.370994 -258.358197
[67,] -129.604427 -2153.370994
[68,] -1975.007793 -129.604427
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -631.297216 -381.996700
2 380.980558 -631.297216
3 754.872113 380.980558
4 -420.894473 754.872113
5 4.822662 -420.894473
6 20.702637 4.822662
7 22.528968 20.702637
8 -81.996508 22.528968
9 -49.843028 -81.996508
10 -186.914513 -49.843028
11 -36.162178 -186.914513
12 -106.975327 -36.162178
13 501.411993 -106.975327
14 -272.982300 501.411993
15 -182.862680 -272.982300
16 -71.141708 -182.862680
17 37.670686 -71.141708
18 -347.338587 37.670686
19 -1294.232626 -347.338587
20 -252.996168 -1294.232626
21 -2621.831784 -252.996168
22 2232.026863 -2621.831784
23 -453.516859 2232.026863
24 42.853132 -453.516859
25 -143.538792 42.853132
26 2559.117391 -143.538792
27 -1797.832710 2559.117391
28 -47.505455 -1797.832710
29 -1445.364329 -47.505455
30 -145.214467 -1445.364329
31 1866.772236 -145.214467
32 453.077688 1866.772236
33 1460.069331 453.077688
34 -381.996700 1460.069331
35 1827.234775 -381.996700
36 465.995968 1827.234775
37 177.428640 465.995968
38 614.888120 177.428640
39 3.881202 614.888120
40 445.047800 3.881202
41 26.723747 445.047800
42 18.166340 26.723747
43 198.060230 18.166340
44 -234.841449 198.060230
45 -23.625158 -234.841449
46 -104.023297 -23.625158
47 -64.993694 -104.023297
48 -126.265109 -64.993694
49 111.401259 -126.265109
50 -247.841400 111.401259
51 376.769705 -247.841400
52 200.630312 376.769705
53 272.346197 200.630312
54 -1033.076622 272.346197
55 -1997.327925 -1033.076622
56 3481.298222 -1997.327925
57 -257.531783 3481.298222
58 41.281764 -257.531783
59 136.637626 41.281764
60 -241.380807 136.637626
61 32.981457 -241.380807
62 -56.937370 32.981457
63 1265.651544 -56.937370
64 225.289966 1265.651544
65 -258.358197 225.289966
66 -2153.370994 -258.358197
67 -129.604427 -2153.370994
68 -1975.007793 -129.604427
> 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/7jdkw1355151226.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/8qzv41355151226.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/9p2841355151226.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/10f0pe1355151226.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/11quto1355151226.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/12bdfz1355151226.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/139kp71355151226.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/14kyfo1355151226.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/15wc0p1355151226.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/168yfq1355151226.tab")
+ }
>
> try(system("convert tmp/1v1d21355151226.ps tmp/1v1d21355151226.png",intern=TRUE))
character(0)
> try(system("convert tmp/2y1wm1355151226.ps tmp/2y1wm1355151226.png",intern=TRUE))
character(0)
> try(system("convert tmp/3poxy1355151226.ps tmp/3poxy1355151226.png",intern=TRUE))
character(0)
> try(system("convert tmp/4e2x81355151226.ps tmp/4e2x81355151226.png",intern=TRUE))
character(0)
> try(system("convert tmp/5f6ev1355151226.ps tmp/5f6ev1355151226.png",intern=TRUE))
character(0)
> try(system("convert tmp/6344p1355151226.ps tmp/6344p1355151226.png",intern=TRUE))
character(0)
> try(system("convert tmp/7jdkw1355151226.ps tmp/7jdkw1355151226.png",intern=TRUE))
character(0)
> try(system("convert tmp/8qzv41355151226.ps tmp/8qzv41355151226.png",intern=TRUE))
character(0)
> try(system("convert tmp/9p2841355151226.ps tmp/9p2841355151226.png",intern=TRUE))
character(0)
> try(system("convert tmp/10f0pe1355151226.ps tmp/10f0pe1355151226.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
5.875 0.786 6.663