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])
+ }
+ }
> 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
2010-I 2010-II 2010-III 2010-IV 2011-I 2011-II 2011-III 2011-IV\r
1 18897 22424 19364 19434 22831 23072 37471 14690
2 17518 22125 18586 18389 22727 22551 36160 13824
3 8632 7653 8225 8405 8344 8695 9197 9477
4 832 554 822 854 830 935 1051 1150
5 3351 3357 3270 3346 3235 3329 3480 3447
6 8 8 3 4 5 5 4 4
7 1 1 1 1 1 1 1 2
8 7 10 11 9 10 9 10 9
9 217 222 204 205 191 197 196 191
10 911 947 918 939 937 967 1007 962
11 1932 1901 1862 1921 1823 1879 1982 2003
12 274 267 270 267 269 271 281 276
13 131 109 87 66 68 64 76 81
14 1708 1668 1738 1715 1726 1771 1861 2079
15 2609 1965 2308 2424 2486 2594 2729 2720
16 133 32 119 89 93 107 102 23
17 2476 1933 2189 2335 2393 2487 2627 2697
18 10 37 23 21 22 27 21 18
19 1510 1616 1378 1605 1534 1654 1421 1650
20 6427 7719 8279 6133 11706 9235 24339 1324
21 3812 5127 5890 4487 7888 6772 9522 3656
22 724 99 154 157 367 153 171 -211
23 1560 1996 1917 1223 2860 1964 13508 -2076
24 156 113 166 116 435 166 975 -178
25 3 3 2 2 15 13 27 13
26 172 380 151 148 141 167 137 120
27 65 1700 163 568 80 2043 768 338
28 593 2931 292 1348 774 631 122 698
29 281 469 227 309 266 267 292 321
30 1191 145 747 874 38 275 423 218
31 72 57 2 32 32 57 16 21
32 113 -6 27 120 32 184 860 604
33 19 86 2 18 3 3 12 246
34 97 11 27 120 32 185 860 381
35 18897 22424 19364 19434 22831 23072 37471 14690
36 16770 18775 17704 16289 21687 20252 35933 13873
37 6132 5145 5705 5818 5817 6171 6504 6749
38 648 299 484 535 511 548 638 758
39 1739 1710 1776 1910 1843 1990 2141 2097
40 160 167 176 193 183 202 163 179
41 621 570 592 743 655 735 851 835
42 804 821 842 831 847 869 886 881
43 3 3 3 3 3 3 3 3
44 150 149 163 140 156 182 238 199
45 549 528 558 354 470 438 450 453
46 95 80 132 -46 88 49 57 62
47 354 353 339 323 308 314 325 322
48 100 95 87 77 73 75 68 70
49 342 343 357 354 339 343 343 305
50 2854 2265 2530 2664 2653 2852 2932 3135
51 167 99 168 132 137 147 132 35
52 2687 2166 2362 2533 2516 2705 2799 3100
53 645 770 634 680 581 675 814 677
54 6113 7729 8065 5931 11602 9279 24726 1473
55 3567 4697 5792 4959 8473 6753 9199 3926
56 472 241 87 262 330 64 172 -142
57 1665 2360 1934 584 2229 1972 13856 -2338
58 328 318 154 6 256 181 1301 -241
59 0 1 0 0 12 10 21 10
60 81 112 99 120 302 300 177 259
61 1322 1286 1317 1325 1314 1322 1308 1370
62 154 143 156 152 144 151 151 171
63 1277 1448 1340 1689 1529 1544 1264 1656
64 1127 2253 486 694 699 1110 1165 1776
65 456 1356 63 2861 89 82 1019 926
66 224 200 149 91 165 216 94 131
67 1444 1990 1445 176 888 2517 413 -264
68 3 8 4 7 40 38 -64 -10
69 1444 2084 1443 187 850 2483 489 -231
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `2010-II` `2010-III` `2010-IV` `2011-I`
13.55398 0.23236 1.35519 0.02169 -0.63595
`2011-II` `2011-III` `2011-IV\\r`
-0.21129 0.08752 0.20833
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-620.04 -67.39 -19.97 75.61 770.05
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.55398 26.48350 0.512 0.610646
`2010-II` 0.23236 0.06670 3.484 0.000921 ***
`2010-III` 1.35519 0.13491 10.045 1.50e-14 ***
`2010-IV` 0.02169 0.07504 0.289 0.773569
`2011-I` -0.63595 0.06426 -9.896 2.66e-14 ***
`2011-II` -0.21129 0.10520 -2.008 0.049028 *
`2011-III` 0.08752 0.02128 4.114 0.000119 ***
`2011-IV\\r` 0.20833 0.05354 3.891 0.000250 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 195.2 on 61 degrees of freedom
Multiple R-squared: 0.9982, Adjusted R-squared: 0.998
F-statistic: 4745 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,] 2.985443e-04 5.970887e-04 0.9997014557
[2,] 1.040974e-04 2.081947e-04 0.9998959026
[3,] 1.441756e-04 2.883513e-04 0.9998558244
[4,] 5.383071e-05 1.076614e-04 0.9999461693
[5,] 7.906568e-06 1.581314e-05 0.9999920934
[6,] 8.413403e-06 1.682681e-05 0.9999915866
[7,] 1.561135e-06 3.122270e-06 0.9999984389
[8,] 2.745603e-07 5.491206e-07 0.9999997254
[9,] 1.869438e-07 3.738876e-07 0.9999998131
[10,] 6.734445e-07 1.346889e-06 0.9999993266
[11,] 5.403290e-06 1.080658e-05 0.9999945967
[12,] 5.369730e-02 1.073946e-01 0.9463026955
[13,] 3.779110e-02 7.558220e-02 0.9622088995
[14,] 8.724514e-02 1.744903e-01 0.9127548625
[15,] 5.738498e-02 1.147700e-01 0.9426150189
[16,] 5.911976e-02 1.182395e-01 0.9408802381
[17,] 2.856404e-01 5.712807e-01 0.7143596332
[18,] 2.847859e-01 5.695717e-01 0.7152141329
[19,] 2.295585e-01 4.591170e-01 0.7704415176
[20,] 4.313203e-01 8.626406e-01 0.5686797102
[21,] 3.632797e-01 7.265594e-01 0.6367202940
[22,] 3.146187e-01 6.292374e-01 0.6853812863
[23,] 2.608070e-01 5.216139e-01 0.7391930287
[24,] 2.263194e-01 4.526388e-01 0.7736805948
[25,] 1.756802e-01 3.513604e-01 0.8243198171
[26,] 5.824577e-01 8.350846e-01 0.4175422817
[27,] 5.363012e-01 9.273976e-01 0.4636988163
[28,] 4.690139e-01 9.380279e-01 0.5309860673
[29,] 4.753612e-01 9.507224e-01 0.5246388247
[30,] 3.989637e-01 7.979275e-01 0.6010362567
[31,] 3.524190e-01 7.048381e-01 0.6475809702
[32,] 2.873789e-01 5.747578e-01 0.7126211211
[33,] 2.233943e-01 4.467885e-01 0.7766057497
[34,] 1.721288e-01 3.442577e-01 0.8278711677
[35,] 1.422437e-01 2.844874e-01 0.8577562876
[36,] 1.090828e-01 2.181657e-01 0.8909171561
[37,] 7.551272e-02 1.510254e-01 0.9244872753
[38,] 4.997905e-02 9.995811e-02 0.9500209455
[39,] 3.217337e-02 6.434674e-02 0.9678266321
[40,] 3.433212e-02 6.866425e-02 0.9656678770
[41,] 2.036013e-02 4.072026e-02 0.9796398684
[42,] 4.731299e-01 9.462597e-01 0.5268701467
[43,] 3.777444e-01 7.554888e-01 0.6222556072
[44,] 8.048314e-01 3.903372e-01 0.1951686089
[45,] 9.993915e-01 1.216963e-03 0.0006084815
[46,] 9.982241e-01 3.551757e-03 0.0017758784
[47,] 9.963682e-01 7.263588e-03 0.0036317942
[48,] 9.857978e-01 2.840441e-02 0.0142022049
> postscript(file="/var/wessaorg/rcomp/tmp/1xluz1355777214.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/2vvmm1355777214.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/3oa5o1355777214.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/4jdci1355777214.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/5m91h1355777214.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
64.0555491 -49.5406444 -124.2185943 -48.9318163 -208.6110245 -8.5123742
7 8 9 10 11 12
-13.8201495 -18.4688883 -22.8962174 -75.3676776 -122.6866850 -27.0469432
13 14 15 16 17 18
6.0270103 -209.7847349 282.0794062 16.8456668 251.6797585 -29.6681805
19 20 21 22 23 24
75.6098252 256.8867428 -620.0378963 770.0491081 -57.7558118 152.1708255
25 26 27 28 29 30
-6.7901398 -49.7295967 -231.8647973 -57.0938438 -22.7121731 112.3093181
31 32 33 34 35 36
68.4159280 -80.2219025 -67.3943669 -53.5042712 64.0555491 84.1984292
37 38 39 40 41 42
93.3879266 124.4645707 -151.8468913 -27.5537334 -19.9716253 -98.2310926
43 44 45 46 47 48
-13.7275114 -46.7316112 -93.4279651 -66.6185815 -41.2988878 -13.4627293
49 50 51 52 53 54
-48.2329192 207.8092717 -0.7514868 195.0726154 -121.5797950 113.2005576
55 56 57 58 59 60
-302.5127236 516.7772065 -421.9858098 169.1111450 -7.9632763 90.6518917
61 62 63 64 65 66
-88.8090337 -32.8452113 -82.6086398 123.3761918 -228.2253349 75.1300810
67 68 69
121.3796670 23.1663380 57.1330089
> postscript(file="/var/wessaorg/rcomp/tmp/68nq51355777214.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 64.0555491 NA
1 -49.5406444 64.0555491
2 -124.2185943 -49.5406444
3 -48.9318163 -124.2185943
4 -208.6110245 -48.9318163
5 -8.5123742 -208.6110245
6 -13.8201495 -8.5123742
7 -18.4688883 -13.8201495
8 -22.8962174 -18.4688883
9 -75.3676776 -22.8962174
10 -122.6866850 -75.3676776
11 -27.0469432 -122.6866850
12 6.0270103 -27.0469432
13 -209.7847349 6.0270103
14 282.0794062 -209.7847349
15 16.8456668 282.0794062
16 251.6797585 16.8456668
17 -29.6681805 251.6797585
18 75.6098252 -29.6681805
19 256.8867428 75.6098252
20 -620.0378963 256.8867428
21 770.0491081 -620.0378963
22 -57.7558118 770.0491081
23 152.1708255 -57.7558118
24 -6.7901398 152.1708255
25 -49.7295967 -6.7901398
26 -231.8647973 -49.7295967
27 -57.0938438 -231.8647973
28 -22.7121731 -57.0938438
29 112.3093181 -22.7121731
30 68.4159280 112.3093181
31 -80.2219025 68.4159280
32 -67.3943669 -80.2219025
33 -53.5042712 -67.3943669
34 64.0555491 -53.5042712
35 84.1984292 64.0555491
36 93.3879266 84.1984292
37 124.4645707 93.3879266
38 -151.8468913 124.4645707
39 -27.5537334 -151.8468913
40 -19.9716253 -27.5537334
41 -98.2310926 -19.9716253
42 -13.7275114 -98.2310926
43 -46.7316112 -13.7275114
44 -93.4279651 -46.7316112
45 -66.6185815 -93.4279651
46 -41.2988878 -66.6185815
47 -13.4627293 -41.2988878
48 -48.2329192 -13.4627293
49 207.8092717 -48.2329192
50 -0.7514868 207.8092717
51 195.0726154 -0.7514868
52 -121.5797950 195.0726154
53 113.2005576 -121.5797950
54 -302.5127236 113.2005576
55 516.7772065 -302.5127236
56 -421.9858098 516.7772065
57 169.1111450 -421.9858098
58 -7.9632763 169.1111450
59 90.6518917 -7.9632763
60 -88.8090337 90.6518917
61 -32.8452113 -88.8090337
62 -82.6086398 -32.8452113
63 123.3761918 -82.6086398
64 -228.2253349 123.3761918
65 75.1300810 -228.2253349
66 121.3796670 75.1300810
67 23.1663380 121.3796670
68 57.1330089 23.1663380
69 NA 57.1330089
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -49.5406444 64.0555491
[2,] -124.2185943 -49.5406444
[3,] -48.9318163 -124.2185943
[4,] -208.6110245 -48.9318163
[5,] -8.5123742 -208.6110245
[6,] -13.8201495 -8.5123742
[7,] -18.4688883 -13.8201495
[8,] -22.8962174 -18.4688883
[9,] -75.3676776 -22.8962174
[10,] -122.6866850 -75.3676776
[11,] -27.0469432 -122.6866850
[12,] 6.0270103 -27.0469432
[13,] -209.7847349 6.0270103
[14,] 282.0794062 -209.7847349
[15,] 16.8456668 282.0794062
[16,] 251.6797585 16.8456668
[17,] -29.6681805 251.6797585
[18,] 75.6098252 -29.6681805
[19,] 256.8867428 75.6098252
[20,] -620.0378963 256.8867428
[21,] 770.0491081 -620.0378963
[22,] -57.7558118 770.0491081
[23,] 152.1708255 -57.7558118
[24,] -6.7901398 152.1708255
[25,] -49.7295967 -6.7901398
[26,] -231.8647973 -49.7295967
[27,] -57.0938438 -231.8647973
[28,] -22.7121731 -57.0938438
[29,] 112.3093181 -22.7121731
[30,] 68.4159280 112.3093181
[31,] -80.2219025 68.4159280
[32,] -67.3943669 -80.2219025
[33,] -53.5042712 -67.3943669
[34,] 64.0555491 -53.5042712
[35,] 84.1984292 64.0555491
[36,] 93.3879266 84.1984292
[37,] 124.4645707 93.3879266
[38,] -151.8468913 124.4645707
[39,] -27.5537334 -151.8468913
[40,] -19.9716253 -27.5537334
[41,] -98.2310926 -19.9716253
[42,] -13.7275114 -98.2310926
[43,] -46.7316112 -13.7275114
[44,] -93.4279651 -46.7316112
[45,] -66.6185815 -93.4279651
[46,] -41.2988878 -66.6185815
[47,] -13.4627293 -41.2988878
[48,] -48.2329192 -13.4627293
[49,] 207.8092717 -48.2329192
[50,] -0.7514868 207.8092717
[51,] 195.0726154 -0.7514868
[52,] -121.5797950 195.0726154
[53,] 113.2005576 -121.5797950
[54,] -302.5127236 113.2005576
[55,] 516.7772065 -302.5127236
[56,] -421.9858098 516.7772065
[57,] 169.1111450 -421.9858098
[58,] -7.9632763 169.1111450
[59,] 90.6518917 -7.9632763
[60,] -88.8090337 90.6518917
[61,] -32.8452113 -88.8090337
[62,] -82.6086398 -32.8452113
[63,] 123.3761918 -82.6086398
[64,] -228.2253349 123.3761918
[65,] 75.1300810 -228.2253349
[66,] 121.3796670 75.1300810
[67,] 23.1663380 121.3796670
[68,] 57.1330089 23.1663380
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -49.5406444 64.0555491
2 -124.2185943 -49.5406444
3 -48.9318163 -124.2185943
4 -208.6110245 -48.9318163
5 -8.5123742 -208.6110245
6 -13.8201495 -8.5123742
7 -18.4688883 -13.8201495
8 -22.8962174 -18.4688883
9 -75.3676776 -22.8962174
10 -122.6866850 -75.3676776
11 -27.0469432 -122.6866850
12 6.0270103 -27.0469432
13 -209.7847349 6.0270103
14 282.0794062 -209.7847349
15 16.8456668 282.0794062
16 251.6797585 16.8456668
17 -29.6681805 251.6797585
18 75.6098252 -29.6681805
19 256.8867428 75.6098252
20 -620.0378963 256.8867428
21 770.0491081 -620.0378963
22 -57.7558118 770.0491081
23 152.1708255 -57.7558118
24 -6.7901398 152.1708255
25 -49.7295967 -6.7901398
26 -231.8647973 -49.7295967
27 -57.0938438 -231.8647973
28 -22.7121731 -57.0938438
29 112.3093181 -22.7121731
30 68.4159280 112.3093181
31 -80.2219025 68.4159280
32 -67.3943669 -80.2219025
33 -53.5042712 -67.3943669
34 64.0555491 -53.5042712
35 84.1984292 64.0555491
36 93.3879266 84.1984292
37 124.4645707 93.3879266
38 -151.8468913 124.4645707
39 -27.5537334 -151.8468913
40 -19.9716253 -27.5537334
41 -98.2310926 -19.9716253
42 -13.7275114 -98.2310926
43 -46.7316112 -13.7275114
44 -93.4279651 -46.7316112
45 -66.6185815 -93.4279651
46 -41.2988878 -66.6185815
47 -13.4627293 -41.2988878
48 -48.2329192 -13.4627293
49 207.8092717 -48.2329192
50 -0.7514868 207.8092717
51 195.0726154 -0.7514868
52 -121.5797950 195.0726154
53 113.2005576 -121.5797950
54 -302.5127236 113.2005576
55 516.7772065 -302.5127236
56 -421.9858098 516.7772065
57 169.1111450 -421.9858098
58 -7.9632763 169.1111450
59 90.6518917 -7.9632763
60 -88.8090337 90.6518917
61 -32.8452113 -88.8090337
62 -82.6086398 -32.8452113
63 123.3761918 -82.6086398
64 -228.2253349 123.3761918
65 75.1300810 -228.2253349
66 121.3796670 75.1300810
67 23.1663380 121.3796670
68 57.1330089 23.1663380
> 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/76z941355777214.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/821w81355777214.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/9u3rc1355777214.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/10fzy21355777214.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/11ee511355777214.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/128gn41355777214.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/13w3iz1355777214.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/14i7nn1355777214.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/15ivfr1355777214.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/163kyc1355777214.tab")
+ }
>
> try(system("convert tmp/1xluz1355777214.ps tmp/1xluz1355777214.png",intern=TRUE))
character(0)
> try(system("convert tmp/2vvmm1355777214.ps tmp/2vvmm1355777214.png",intern=TRUE))
character(0)
> try(system("convert tmp/3oa5o1355777214.ps tmp/3oa5o1355777214.png",intern=TRUE))
character(0)
> try(system("convert tmp/4jdci1355777214.ps tmp/4jdci1355777214.png",intern=TRUE))
character(0)
> try(system("convert tmp/5m91h1355777214.ps tmp/5m91h1355777214.png",intern=TRUE))
character(0)
> try(system("convert tmp/68nq51355777214.ps tmp/68nq51355777214.png",intern=TRUE))
character(0)
> try(system("convert tmp/76z941355777214.ps tmp/76z941355777214.png",intern=TRUE))
character(0)
> try(system("convert tmp/821w81355777214.ps tmp/821w81355777214.png",intern=TRUE))
character(0)
> try(system("convert tmp/9u3rc1355777214.ps tmp/9u3rc1355777214.png",intern=TRUE))
character(0)
> try(system("convert tmp/10fzy21355777214.ps tmp/10fzy21355777214.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.249 1.185 7.427