R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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(116222
+ ,344744
+ ,31899
+ ,492865
+ ,110924
+ ,338653
+ ,31384
+ ,480961
+ ,103753
+ ,327532
+ ,30650
+ ,461935
+ ,99983
+ ,326225
+ ,30400
+ ,456608
+ ,93302
+ ,318672
+ ,30003
+ ,441977
+ ,91496
+ ,317756
+ ,29896
+ ,439148
+ ,119321
+ ,337302
+ ,31557
+ ,488180
+ ,139261
+ ,349420
+ ,31883
+ ,520564
+ ,133739
+ ,336923
+ ,30830
+ ,501492
+ ,123913
+ ,330758
+ ,30354
+ ,485025
+ ,113438
+ ,321002
+ ,29756
+ ,464196
+ ,109416
+ ,320820
+ ,29934
+ ,460170
+ ,109406
+ ,327032
+ ,30599
+ ,467037
+ ,105645
+ ,324047
+ ,30378
+ ,460070
+ ,101328
+ ,316735
+ ,29925
+ ,447988
+ ,97686
+ ,315710
+ ,29471
+ ,442867
+ ,93093
+ ,313427
+ ,29567
+ ,436087
+ ,91382
+ ,310527
+ ,29419
+ ,431328
+ ,122257
+ ,330962
+ ,30796
+ ,484015
+ ,139183
+ ,339015
+ ,31475
+ ,509673
+ ,139887
+ ,341332
+ ,31708
+ ,512927
+ ,131822
+ ,339092
+ ,31917
+ ,502831
+ ,116805
+ ,323308
+ ,30871
+ ,470984
+ ,113706
+ ,325849
+ ,31512
+ ,471067
+ ,113012
+ ,330675
+ ,32362
+ ,476049
+ ,110452
+ ,332225
+ ,31928
+ ,474605
+ ,107005
+ ,331735
+ ,31699
+ ,470439
+ ,102841
+ ,328047
+ ,30363
+ ,461251
+ ,98173
+ ,326165
+ ,30386
+ ,454724
+ ,98181
+ ,327081
+ ,30364
+ ,455626
+ ,137277
+ ,346764
+ ,32806
+ ,516847
+ ,147579
+ ,344190
+ ,33423
+ ,525192
+ ,146571
+ ,343333
+ ,33071
+ ,522975
+ ,138920
+ ,345777
+ ,33888
+ ,518585
+ ,130340
+ ,344094
+ ,34805
+ ,509239
+ ,128140
+ ,348609
+ ,35489
+ ,512238
+ ,127059
+ ,354846
+ ,37259
+ ,519164
+ ,122860
+ ,356427
+ ,37722
+ ,517009
+ ,117702
+ ,353467
+ ,38764
+ ,509933
+ ,113537
+ ,355996
+ ,39594
+ ,509127
+ ,108366
+ ,352487
+ ,40004
+ ,500857
+ ,111078
+ ,355178
+ ,40715
+ ,506971
+ ,150739
+ ,374556
+ ,44028
+ ,569323
+ ,159129
+ ,375021
+ ,45564
+ ,579714
+ ,157928
+ ,375787
+ ,44277
+ ,577992
+ ,147768
+ ,372720
+ ,44976
+ ,565464
+ ,137507
+ ,364431
+ ,45406
+ ,547344
+ ,136919
+ ,370490
+ ,47379
+ ,554788
+ ,136151
+ ,376974
+ ,49200
+ ,562325
+ ,133001
+ ,377632
+ ,50221
+ ,560854
+ ,125554
+ ,378205
+ ,51573
+ ,555332
+ ,119647
+ ,370861
+ ,53091
+ ,543599
+ ,114158
+ ,369167
+ ,53337
+ ,536662
+ ,116193
+ ,371551
+ ,54978
+ ,542722
+ ,152803
+ ,382842
+ ,57885
+ ,593530
+ ,161761
+ ,381903
+ ,67099
+ ,610763
+ ,160942
+ ,384502
+ ,67169
+ ,612613
+ ,149470
+ ,392058
+ ,69796
+ ,611324
+ ,139208
+ ,384359
+ ,70600
+ ,594167
+ ,134588
+ ,388884
+ ,71982
+ ,595454
+ ,130322
+ ,386586
+ ,73957
+ ,590865
+ ,126611
+ ,387495
+ ,75273
+ ,589379
+ ,122401
+ ,385705
+ ,76322
+ ,584428
+ ,117352
+ ,378670
+ ,77078
+ ,573100
+ ,112135
+ ,377367
+ ,77954
+ ,567456
+ ,112879
+ ,376911
+ ,79238
+ ,569028
+ ,148729
+ ,389827
+ ,82179
+ ,620735
+ ,157230
+ ,387820
+ ,83834
+ ,628884
+ ,157221
+ ,387267
+ ,83744
+ ,628232
+ ,146681
+ ,380575
+ ,84861
+ ,612117
+ ,136524
+ ,372402
+ ,86478
+ ,595404
+ ,132111
+ ,376740
+ ,88290
+ ,597141)
+ ,dim=c(4
+ ,72)
+ ,dimnames=list(c('-25'
+ ,'25-50'
+ ,'50+'
+ ,'totaal')
+ ,1:72))
> y <- array(NA,dim=c(4,72),dimnames=list(c('-25','25-50','50+','totaal'),1:72))
> 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 = 'Include Monthly Dummies'
> par1 = '4'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
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
totaal -25 25-50 50+ M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 492865 116222 344744 31899 1 0 0 0 0 0 0 0 0 0 0
2 480961 110924 338653 31384 0 1 0 0 0 0 0 0 0 0 0
3 461935 103753 327532 30650 0 0 1 0 0 0 0 0 0 0 0
4 456608 99983 326225 30400 0 0 0 1 0 0 0 0 0 0 0
5 441977 93302 318672 30003 0 0 0 0 1 0 0 0 0 0 0
6 439148 91496 317756 29896 0 0 0 0 0 1 0 0 0 0 0
7 488180 119321 337302 31557 0 0 0 0 0 0 1 0 0 0 0
8 520564 139261 349420 31883 0 0 0 0 0 0 0 1 0 0 0
9 501492 133739 336923 30830 0 0 0 0 0 0 0 0 1 0 0
10 485025 123913 330758 30354 0 0 0 0 0 0 0 0 0 1 0
11 464196 113438 321002 29756 0 0 0 0 0 0 0 0 0 0 1
12 460170 109416 320820 29934 0 0 0 0 0 0 0 0 0 0 0
13 467037 109406 327032 30599 1 0 0 0 0 0 0 0 0 0 0
14 460070 105645 324047 30378 0 1 0 0 0 0 0 0 0 0 0
15 447988 101328 316735 29925 0 0 1 0 0 0 0 0 0 0 0
16 442867 97686 315710 29471 0 0 0 1 0 0 0 0 0 0 0
17 436087 93093 313427 29567 0 0 0 0 1 0 0 0 0 0 0
18 431328 91382 310527 29419 0 0 0 0 0 1 0 0 0 0 0
19 484015 122257 330962 30796 0 0 0 0 0 0 1 0 0 0 0
20 509673 139183 339015 31475 0 0 0 0 0 0 0 1 0 0 0
21 512927 139887 341332 31708 0 0 0 0 0 0 0 0 1 0 0
22 502831 131822 339092 31917 0 0 0 0 0 0 0 0 0 1 0
23 470984 116805 323308 30871 0 0 0 0 0 0 0 0 0 0 1
24 471067 113706 325849 31512 0 0 0 0 0 0 0 0 0 0 0
25 476049 113012 330675 32362 1 0 0 0 0 0 0 0 0 0 0
26 474605 110452 332225 31928 0 1 0 0 0 0 0 0 0 0 0
27 470439 107005 331735 31699 0 0 1 0 0 0 0 0 0 0 0
28 461251 102841 328047 30363 0 0 0 1 0 0 0 0 0 0 0
29 454724 98173 326165 30386 0 0 0 0 1 0 0 0 0 0 0
30 455626 98181 327081 30364 0 0 0 0 0 1 0 0 0 0 0
31 516847 137277 346764 32806 0 0 0 0 0 0 1 0 0 0 0
32 525192 147579 344190 33423 0 0 0 0 0 0 0 1 0 0 0
33 522975 146571 343333 33071 0 0 0 0 0 0 0 0 1 0 0
34 518585 138920 345777 33888 0 0 0 0 0 0 0 0 0 1 0
35 509239 130340 344094 34805 0 0 0 0 0 0 0 0 0 0 1
36 512238 128140 348609 35489 0 0 0 0 0 0 0 0 0 0 0
37 519164 127059 354846 37259 1 0 0 0 0 0 0 0 0 0 0
38 517009 122860 356427 37722 0 1 0 0 0 0 0 0 0 0 0
39 509933 117702 353467 38764 0 0 1 0 0 0 0 0 0 0 0
40 509127 113537 355996 39594 0 0 0 1 0 0 0 0 0 0 0
41 500857 108366 352487 40004 0 0 0 0 1 0 0 0 0 0 0
42 506971 111078 355178 40715 0 0 0 0 0 1 0 0 0 0 0
43 569323 150739 374556 44028 0 0 0 0 0 0 1 0 0 0 0
44 579714 159129 375021 45564 0 0 0 0 0 0 0 1 0 0 0
45 577992 157928 375787 44277 0 0 0 0 0 0 0 0 1 0 0
46 565464 147768 372720 44976 0 0 0 0 0 0 0 0 0 1 0
47 547344 137507 364431 45406 0 0 0 0 0 0 0 0 0 0 1
48 554788 136919 370490 47379 0 0 0 0 0 0 0 0 0 0 0
49 562325 136151 376974 49200 1 0 0 0 0 0 0 0 0 0 0
50 560854 133001 377632 50221 0 1 0 0 0 0 0 0 0 0 0
51 555332 125554 378205 51573 0 0 1 0 0 0 0 0 0 0 0
52 543599 119647 370861 53091 0 0 0 1 0 0 0 0 0 0 0
53 536662 114158 369167 53337 0 0 0 0 1 0 0 0 0 0 0
54 542722 116193 371551 54978 0 0 0 0 0 1 0 0 0 0 0
55 593530 152803 382842 57885 0 0 0 0 0 0 1 0 0 0 0
56 610763 161761 381903 67099 0 0 0 0 0 0 0 1 0 0 0
57 612613 160942 384502 67169 0 0 0 0 0 0 0 0 1 0 0
58 611324 149470 392058 69796 0 0 0 0 0 0 0 0 0 1 0
59 594167 139208 384359 70600 0 0 0 0 0 0 0 0 0 0 1
60 595454 134588 388884 71982 0 0 0 0 0 0 0 0 0 0 0
61 590865 130322 386586 73957 1 0 0 0 0 0 0 0 0 0 0
62 589379 126611 387495 75273 0 1 0 0 0 0 0 0 0 0 0
63 584428 122401 385705 76322 0 0 1 0 0 0 0 0 0 0 0
64 573100 117352 378670 77078 0 0 0 1 0 0 0 0 0 0 0
65 567456 112135 377367 77954 0 0 0 0 1 0 0 0 0 0 0
66 569028 112879 376911 79238 0 0 0 0 0 1 0 0 0 0 0
67 620735 148729 389827 82179 0 0 0 0 0 0 1 0 0 0 0
68 628884 157230 387820 83834 0 0 0 0 0 0 0 1 0 0 0
69 628232 157221 387267 83744 0 0 0 0 0 0 0 0 1 0 0
70 612117 146681 380575 84861 0 0 0 0 0 0 0 0 0 1 0
71 595404 136524 372402 86478 0 0 0 0 0 0 0 0 0 0 1
72 597141 132111 376740 88290 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `-25` `25-50` `50+` M1 M2
-4.939e-10 1.000e+00 1.000e+00 1.000e+00 1.313e-10 -8.132e-11
M3 M4 M5 M6 M7 M8
-3.910e-11 -7.167e-11 -8.074e-11 -7.998e-11 3.559e-11 8.867e-11
M9 M10 M11
8.843e-11 5.061e-11 2.511e-11
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.175e-10 -1.269e-11 -1.132e-12 1.562e-11 7.753e-10
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.939e-10 4.619e-10 -1.069e+00 0.2894
`-25` 1.000e+00 4.810e-15 2.079e+14 <2e-16 ***
`25-50` 1.000e+00 2.883e-15 3.469e+14 <2e-16 ***
`50+` 1.000e+00 1.835e-15 5.450e+14 <2e-16 ***
M1 1.313e-10 7.256e-11 1.810e+00 0.0756 .
M2 -8.132e-11 7.870e-11 -1.033e+00 0.3058
M3 -3.910e-11 8.662e-11 -4.510e-01 0.6534
M4 -7.167e-11 9.540e-11 -7.510e-01 0.4556
M5 -8.074e-11 1.083e-10 -7.450e-01 0.4591
M6 -7.998e-11 1.073e-10 -7.460e-01 0.4589
M7 3.559e-11 8.204e-11 4.340e-01 0.6660
M8 8.867e-11 1.200e-10 7.390e-01 0.4629
M9 8.843e-11 1.173e-10 7.540e-01 0.4542
M10 5.061e-11 8.678e-11 5.830e-01 0.5621
M11 2.511e-11 7.239e-11 3.470e-01 0.7299
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.184e-10 on 57 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 1.145e+30 on 14 and 57 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,] 9.999803e-01 3.943826e-05 1.971913e-05
[2,] 9.696424e-01 6.071524e-02 3.035762e-02
[3,] 9.474395e-03 1.894879e-02 9.905256e-01
[4,] 1.000000e+00 8.024867e-12 4.012434e-12
[5,] 9.993714e-01 1.257234e-03 6.286172e-04
[6,] 9.999886e-01 2.280694e-05 1.140347e-05
[7,] 7.973048e-01 4.053904e-01 2.026952e-01
[8,] 1.724260e-04 3.448520e-04 9.998276e-01
[9,] 6.851314e-01 6.297371e-01 3.148686e-01
[10,] 2.956627e-03 5.913253e-03 9.970434e-01
[11,] 9.999378e-01 1.244417e-04 6.222083e-05
[12,] 1.000000e+00 3.337972e-10 1.668986e-10
[13,] 1.305142e-05 2.610284e-05 9.999869e-01
[14,] 1.624453e-01 3.248905e-01 8.375547e-01
[15,] 1.007283e-01 2.014567e-01 8.992717e-01
[16,] 4.980808e-01 9.961616e-01 5.019192e-01
[17,] 5.381047e-08 1.076209e-07 9.999999e-01
[18,] 9.973314e-01 5.337164e-03 2.668582e-03
[19,] 4.270281e-03 8.540561e-03 9.957297e-01
[20,] 6.899439e-03 1.379888e-02 9.931006e-01
[21,] 3.283646e-06 6.567293e-06 9.999967e-01
[22,] 7.372554e-13 1.474511e-12 1.000000e+00
[23,] 9.871266e-01 2.574671e-02 1.287336e-02
[24,] 9.999645e-01 7.102440e-05 3.551220e-05
[25,] 6.422908e-01 7.154183e-01 3.577092e-01
[26,] 3.208019e-07 6.416039e-07 9.999997e-01
[27,] 9.999911e-01 1.777282e-05 8.886412e-06
[28,] 1.013329e-01 2.026658e-01 8.986671e-01
[29,] 3.011982e-01 6.023965e-01 6.988018e-01
[30,] 5.171312e-01 9.657375e-01 4.828688e-01
[31,] 9.999719e-01 5.620241e-05 2.810120e-05
[32,] 9.991891e-01 1.621886e-03 8.109429e-04
[33,] 2.511771e-38 5.023541e-38 1.000000e+00
[34,] 6.415104e-01 7.169792e-01 3.584896e-01
[35,] 2.377690e-01 4.755381e-01 7.622310e-01
[36,] 9.989709e-01 2.058219e-03 1.029109e-03
[37,] 9.777259e-01 4.454811e-02 2.227405e-02
> postscript(file="/var/www/html/rcomp/tmp/1h4ke1258294510.ps",horizontal=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/www/html/rcomp/tmp/2chle1258294510.ps",horizontal=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/www/html/rcomp/tmp/3uwrm1258294510.ps",horizontal=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/www/html/rcomp/tmp/4cga61258294510.ps",horizontal=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/www/html/rcomp/tmp/59pkf1258294510.ps",horizontal=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 = 72
Frequency = 1
1 2 3 4 5
7.752777e-10 -2.174881e-10 6.673452e-11 -3.684580e-11 1.248354e-12
6 7 8 9 10
-8.513859e-12 -4.420815e-11 -3.601828e-11 -3.043102e-11 -2.369481e-11
11 12 13 14 15
-1.874504e-11 -1.131868e-11 -1.649848e-10 4.702056e-11 1.998117e-13
16 17 18 19 20
1.231545e-11 1.578788e-11 1.650221e-11 -1.326645e-11 -1.475311e-11
21 22 23 24 25
-1.249583e-11 -2.301949e-12 1.543314e-13 -6.035131e-12 -1.536978e-10
26 27 28 29 30
4.291508e-11 -9.374128e-12 8.594361e-12 2.481750e-12 -2.453898e-12
31 32 33 34 35
3.833265e-12 1.731311e-11 9.247733e-12 1.524665e-11 2.326931e-13
36 37 38 39 40
4.292784e-12 -1.468653e-10 4.615940e-11 -1.174367e-11 -4.289610e-12
41 42 43 44 45
-1.081246e-11 -3.797290e-12 1.067812e-11 -1.427482e-12 -1.503617e-11
46 47 48 49 50
-6.905315e-12 -6.104573e-12 3.888277e-12 -1.465679e-10 5.296613e-11
51 52 53 54 55
-2.282514e-11 4.664194e-12 -7.868106e-12 -5.542843e-12 1.614625e-11
56 57 58 59 60
2.116723e-11 2.333708e-11 -8.976100e-12 -8.704678e-12 -1.365511e-11
61 62 63 64 65
-1.631619e-10 2.842691e-11 -2.299139e-11 1.556140e-11 -8.374186e-13
66 67 68 69 70
3.805683e-12 2.681696e-11 1.371853e-11 2.537821e-11 2.663152e-11
71 72
3.316727e-11 2.282786e-11
> postscript(file="/var/www/html/rcomp/tmp/6hws91258294510.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 7.752777e-10 NA
1 -2.174881e-10 7.752777e-10
2 6.673452e-11 -2.174881e-10
3 -3.684580e-11 6.673452e-11
4 1.248354e-12 -3.684580e-11
5 -8.513859e-12 1.248354e-12
6 -4.420815e-11 -8.513859e-12
7 -3.601828e-11 -4.420815e-11
8 -3.043102e-11 -3.601828e-11
9 -2.369481e-11 -3.043102e-11
10 -1.874504e-11 -2.369481e-11
11 -1.131868e-11 -1.874504e-11
12 -1.649848e-10 -1.131868e-11
13 4.702056e-11 -1.649848e-10
14 1.998117e-13 4.702056e-11
15 1.231545e-11 1.998117e-13
16 1.578788e-11 1.231545e-11
17 1.650221e-11 1.578788e-11
18 -1.326645e-11 1.650221e-11
19 -1.475311e-11 -1.326645e-11
20 -1.249583e-11 -1.475311e-11
21 -2.301949e-12 -1.249583e-11
22 1.543314e-13 -2.301949e-12
23 -6.035131e-12 1.543314e-13
24 -1.536978e-10 -6.035131e-12
25 4.291508e-11 -1.536978e-10
26 -9.374128e-12 4.291508e-11
27 8.594361e-12 -9.374128e-12
28 2.481750e-12 8.594361e-12
29 -2.453898e-12 2.481750e-12
30 3.833265e-12 -2.453898e-12
31 1.731311e-11 3.833265e-12
32 9.247733e-12 1.731311e-11
33 1.524665e-11 9.247733e-12
34 2.326931e-13 1.524665e-11
35 4.292784e-12 2.326931e-13
36 -1.468653e-10 4.292784e-12
37 4.615940e-11 -1.468653e-10
38 -1.174367e-11 4.615940e-11
39 -4.289610e-12 -1.174367e-11
40 -1.081246e-11 -4.289610e-12
41 -3.797290e-12 -1.081246e-11
42 1.067812e-11 -3.797290e-12
43 -1.427482e-12 1.067812e-11
44 -1.503617e-11 -1.427482e-12
45 -6.905315e-12 -1.503617e-11
46 -6.104573e-12 -6.905315e-12
47 3.888277e-12 -6.104573e-12
48 -1.465679e-10 3.888277e-12
49 5.296613e-11 -1.465679e-10
50 -2.282514e-11 5.296613e-11
51 4.664194e-12 -2.282514e-11
52 -7.868106e-12 4.664194e-12
53 -5.542843e-12 -7.868106e-12
54 1.614625e-11 -5.542843e-12
55 2.116723e-11 1.614625e-11
56 2.333708e-11 2.116723e-11
57 -8.976100e-12 2.333708e-11
58 -8.704678e-12 -8.976100e-12
59 -1.365511e-11 -8.704678e-12
60 -1.631619e-10 -1.365511e-11
61 2.842691e-11 -1.631619e-10
62 -2.299139e-11 2.842691e-11
63 1.556140e-11 -2.299139e-11
64 -8.374186e-13 1.556140e-11
65 3.805683e-12 -8.374186e-13
66 2.681696e-11 3.805683e-12
67 1.371853e-11 2.681696e-11
68 2.537821e-11 1.371853e-11
69 2.663152e-11 2.537821e-11
70 3.316727e-11 2.663152e-11
71 2.282786e-11 3.316727e-11
72 NA 2.282786e-11
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.174881e-10 7.752777e-10
[2,] 6.673452e-11 -2.174881e-10
[3,] -3.684580e-11 6.673452e-11
[4,] 1.248354e-12 -3.684580e-11
[5,] -8.513859e-12 1.248354e-12
[6,] -4.420815e-11 -8.513859e-12
[7,] -3.601828e-11 -4.420815e-11
[8,] -3.043102e-11 -3.601828e-11
[9,] -2.369481e-11 -3.043102e-11
[10,] -1.874504e-11 -2.369481e-11
[11,] -1.131868e-11 -1.874504e-11
[12,] -1.649848e-10 -1.131868e-11
[13,] 4.702056e-11 -1.649848e-10
[14,] 1.998117e-13 4.702056e-11
[15,] 1.231545e-11 1.998117e-13
[16,] 1.578788e-11 1.231545e-11
[17,] 1.650221e-11 1.578788e-11
[18,] -1.326645e-11 1.650221e-11
[19,] -1.475311e-11 -1.326645e-11
[20,] -1.249583e-11 -1.475311e-11
[21,] -2.301949e-12 -1.249583e-11
[22,] 1.543314e-13 -2.301949e-12
[23,] -6.035131e-12 1.543314e-13
[24,] -1.536978e-10 -6.035131e-12
[25,] 4.291508e-11 -1.536978e-10
[26,] -9.374128e-12 4.291508e-11
[27,] 8.594361e-12 -9.374128e-12
[28,] 2.481750e-12 8.594361e-12
[29,] -2.453898e-12 2.481750e-12
[30,] 3.833265e-12 -2.453898e-12
[31,] 1.731311e-11 3.833265e-12
[32,] 9.247733e-12 1.731311e-11
[33,] 1.524665e-11 9.247733e-12
[34,] 2.326931e-13 1.524665e-11
[35,] 4.292784e-12 2.326931e-13
[36,] -1.468653e-10 4.292784e-12
[37,] 4.615940e-11 -1.468653e-10
[38,] -1.174367e-11 4.615940e-11
[39,] -4.289610e-12 -1.174367e-11
[40,] -1.081246e-11 -4.289610e-12
[41,] -3.797290e-12 -1.081246e-11
[42,] 1.067812e-11 -3.797290e-12
[43,] -1.427482e-12 1.067812e-11
[44,] -1.503617e-11 -1.427482e-12
[45,] -6.905315e-12 -1.503617e-11
[46,] -6.104573e-12 -6.905315e-12
[47,] 3.888277e-12 -6.104573e-12
[48,] -1.465679e-10 3.888277e-12
[49,] 5.296613e-11 -1.465679e-10
[50,] -2.282514e-11 5.296613e-11
[51,] 4.664194e-12 -2.282514e-11
[52,] -7.868106e-12 4.664194e-12
[53,] -5.542843e-12 -7.868106e-12
[54,] 1.614625e-11 -5.542843e-12
[55,] 2.116723e-11 1.614625e-11
[56,] 2.333708e-11 2.116723e-11
[57,] -8.976100e-12 2.333708e-11
[58,] -8.704678e-12 -8.976100e-12
[59,] -1.365511e-11 -8.704678e-12
[60,] -1.631619e-10 -1.365511e-11
[61,] 2.842691e-11 -1.631619e-10
[62,] -2.299139e-11 2.842691e-11
[63,] 1.556140e-11 -2.299139e-11
[64,] -8.374186e-13 1.556140e-11
[65,] 3.805683e-12 -8.374186e-13
[66,] 2.681696e-11 3.805683e-12
[67,] 1.371853e-11 2.681696e-11
[68,] 2.537821e-11 1.371853e-11
[69,] 2.663152e-11 2.537821e-11
[70,] 3.316727e-11 2.663152e-11
[71,] 2.282786e-11 3.316727e-11
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.174881e-10 7.752777e-10
2 6.673452e-11 -2.174881e-10
3 -3.684580e-11 6.673452e-11
4 1.248354e-12 -3.684580e-11
5 -8.513859e-12 1.248354e-12
6 -4.420815e-11 -8.513859e-12
7 -3.601828e-11 -4.420815e-11
8 -3.043102e-11 -3.601828e-11
9 -2.369481e-11 -3.043102e-11
10 -1.874504e-11 -2.369481e-11
11 -1.131868e-11 -1.874504e-11
12 -1.649848e-10 -1.131868e-11
13 4.702056e-11 -1.649848e-10
14 1.998117e-13 4.702056e-11
15 1.231545e-11 1.998117e-13
16 1.578788e-11 1.231545e-11
17 1.650221e-11 1.578788e-11
18 -1.326645e-11 1.650221e-11
19 -1.475311e-11 -1.326645e-11
20 -1.249583e-11 -1.475311e-11
21 -2.301949e-12 -1.249583e-11
22 1.543314e-13 -2.301949e-12
23 -6.035131e-12 1.543314e-13
24 -1.536978e-10 -6.035131e-12
25 4.291508e-11 -1.536978e-10
26 -9.374128e-12 4.291508e-11
27 8.594361e-12 -9.374128e-12
28 2.481750e-12 8.594361e-12
29 -2.453898e-12 2.481750e-12
30 3.833265e-12 -2.453898e-12
31 1.731311e-11 3.833265e-12
32 9.247733e-12 1.731311e-11
33 1.524665e-11 9.247733e-12
34 2.326931e-13 1.524665e-11
35 4.292784e-12 2.326931e-13
36 -1.468653e-10 4.292784e-12
37 4.615940e-11 -1.468653e-10
38 -1.174367e-11 4.615940e-11
39 -4.289610e-12 -1.174367e-11
40 -1.081246e-11 -4.289610e-12
41 -3.797290e-12 -1.081246e-11
42 1.067812e-11 -3.797290e-12
43 -1.427482e-12 1.067812e-11
44 -1.503617e-11 -1.427482e-12
45 -6.905315e-12 -1.503617e-11
46 -6.104573e-12 -6.905315e-12
47 3.888277e-12 -6.104573e-12
48 -1.465679e-10 3.888277e-12
49 5.296613e-11 -1.465679e-10
50 -2.282514e-11 5.296613e-11
51 4.664194e-12 -2.282514e-11
52 -7.868106e-12 4.664194e-12
53 -5.542843e-12 -7.868106e-12
54 1.614625e-11 -5.542843e-12
55 2.116723e-11 1.614625e-11
56 2.333708e-11 2.116723e-11
57 -8.976100e-12 2.333708e-11
58 -8.704678e-12 -8.976100e-12
59 -1.365511e-11 -8.704678e-12
60 -1.631619e-10 -1.365511e-11
61 2.842691e-11 -1.631619e-10
62 -2.299139e-11 2.842691e-11
63 1.556140e-11 -2.299139e-11
64 -8.374186e-13 1.556140e-11
65 3.805683e-12 -8.374186e-13
66 2.681696e-11 3.805683e-12
67 1.371853e-11 2.681696e-11
68 2.537821e-11 1.371853e-11
69 2.663152e-11 2.537821e-11
70 3.316727e-11 2.663152e-11
71 2.282786e-11 3.316727e-11
> 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/www/html/rcomp/tmp/7vb0d1258294510.ps",horizontal=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/www/html/rcomp/tmp/8b5b91258294510.ps",horizontal=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/www/html/rcomp/tmp/9j3yw1258294510.ps",horizontal=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/www/html/rcomp/tmp/10byl71258294510.ps",horizontal=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/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/www/html/rcomp/tmp/11mpej1258294510.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/www/html/rcomp/tmp/127xxh1258294510.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/www/html/rcomp/tmp/130eaj1258294510.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/www/html/rcomp/tmp/14x8lf1258294510.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/www/html/rcomp/tmp/15bqry1258294510.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/www/html/rcomp/tmp/163g1r1258294510.tab")
+ }
>
> system("convert tmp/1h4ke1258294510.ps tmp/1h4ke1258294510.png")
> system("convert tmp/2chle1258294510.ps tmp/2chle1258294510.png")
> system("convert tmp/3uwrm1258294510.ps tmp/3uwrm1258294510.png")
> system("convert tmp/4cga61258294510.ps tmp/4cga61258294510.png")
> system("convert tmp/59pkf1258294510.ps tmp/59pkf1258294510.png")
> system("convert tmp/6hws91258294510.ps tmp/6hws91258294510.png")
> system("convert tmp/7vb0d1258294510.ps tmp/7vb0d1258294510.png")
> system("convert tmp/8b5b91258294510.ps tmp/8b5b91258294510.png")
> system("convert tmp/9j3yw1258294510.ps tmp/9j3yw1258294510.png")
> system("convert tmp/10byl71258294510.ps tmp/10byl71258294510.png")
>
>
> proc.time()
user system elapsed
2.573 1.608 3.082