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 = 'Do not include Seasonal 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+
1 492865 116222 344744 31899
2 480961 110924 338653 31384
3 461935 103753 327532 30650
4 456608 99983 326225 30400
5 441977 93302 318672 30003
6 439148 91496 317756 29896
7 488180 119321 337302 31557
8 520564 139261 349420 31883
9 501492 133739 336923 30830
10 485025 123913 330758 30354
11 464196 113438 321002 29756
12 460170 109416 320820 29934
13 467037 109406 327032 30599
14 460070 105645 324047 30378
15 447988 101328 316735 29925
16 442867 97686 315710 29471
17 436087 93093 313427 29567
18 431328 91382 310527 29419
19 484015 122257 330962 30796
20 509673 139183 339015 31475
21 512927 139887 341332 31708
22 502831 131822 339092 31917
23 470984 116805 323308 30871
24 471067 113706 325849 31512
25 476049 113012 330675 32362
26 474605 110452 332225 31928
27 470439 107005 331735 31699
28 461251 102841 328047 30363
29 454724 98173 326165 30386
30 455626 98181 327081 30364
31 516847 137277 346764 32806
32 525192 147579 344190 33423
33 522975 146571 343333 33071
34 518585 138920 345777 33888
35 509239 130340 344094 34805
36 512238 128140 348609 35489
37 519164 127059 354846 37259
38 517009 122860 356427 37722
39 509933 117702 353467 38764
40 509127 113537 355996 39594
41 500857 108366 352487 40004
42 506971 111078 355178 40715
43 569323 150739 374556 44028
44 579714 159129 375021 45564
45 577992 157928 375787 44277
46 565464 147768 372720 44976
47 547344 137507 364431 45406
48 554788 136919 370490 47379
49 562325 136151 376974 49200
50 560854 133001 377632 50221
51 555332 125554 378205 51573
52 543599 119647 370861 53091
53 536662 114158 369167 53337
54 542722 116193 371551 54978
55 593530 152803 382842 57885
56 610763 161761 381903 67099
57 612613 160942 384502 67169
58 611324 149470 392058 69796
59 594167 139208 384359 70600
60 595454 134588 388884 71982
61 590865 130322 386586 73957
62 589379 126611 387495 75273
63 584428 122401 385705 76322
64 573100 117352 378670 77078
65 567456 112135 377367 77954
66 569028 112879 376911 79238
67 620735 148729 389827 82179
68 628884 157230 387820 83834
69 628232 157221 387267 83744
70 612117 146681 380575 84861
71 595404 136524 372402 86478
72 597141 132111 376740 88290
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `-25` `25-50` `50+`
-1.921e-10 1.000e+00 1.000e+00 1.000e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.680e-10 -1.932e-11 -1.250e-11 -3.068e-12 9.267e-10
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.921e-10 3.830e-10 -5.020e-01 0.618
`-25` 1.000e+00 1.114e-15 8.975e+14 <2e-16 ***
`25-50` 1.000e+00 1.477e-15 6.770e+14 <2e-16 ***
`50+` 1.000e+00 1.540e-15 6.492e+14 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.186e-10 on 68 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 5.331e+30 on 3 and 68 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,] 8.973984e-03 1.794797e-02 9.910260e-01
[2,] 2.550680e-05 5.101361e-05 9.999745e-01
[3,] 4.257988e-04 8.515977e-04 9.995742e-01
[4,] 3.315931e-07 6.631861e-07 9.999997e-01
[5,] 4.461802e-02 8.923603e-02 9.553820e-01
[6,] 1.469152e-01 2.938303e-01 8.530848e-01
[7,] 6.847833e-02 1.369567e-01 9.315217e-01
[8,] 6.440636e-03 1.288127e-02 9.935594e-01
[9,] 9.097137e-01 1.805727e-01 9.028634e-02
[10,] 9.671877e-01 6.562457e-02 3.281229e-02
[11,] 9.978816e-01 4.236808e-03 2.118404e-03
[12,] 9.999983e-01 3.358702e-06 1.679351e-06
[13,] 9.785248e-01 4.295049e-02 2.147524e-02
[14,] 8.971769e-07 1.794354e-06 9.999991e-01
[15,] 1.000000e+00 1.263778e-16 6.318892e-17
[16,] 9.999341e-01 1.318678e-04 6.593390e-05
[17,] 9.999999e-01 2.196109e-07 1.098054e-07
[18,] 7.766445e-01 4.467111e-01 2.233555e-01
[19,] 3.968068e-08 7.936137e-08 1.000000e+00
[20,] 7.097629e-01 5.804742e-01 2.902371e-01
[21,] 2.282985e-04 4.565969e-04 9.997717e-01
[22,] 9.999997e-01 6.092635e-07 3.046318e-07
[23,] 1.000000e+00 5.709800e-15 2.854900e-15
[24,] 4.271938e-09 8.543875e-09 1.000000e+00
[25,] 1.120817e-01 2.241633e-01 8.879183e-01
[26,] 5.163442e-02 1.032688e-01 9.483656e-01
[27,] 5.450652e-01 9.098695e-01 4.549348e-01
[28,] 9.499879e-11 1.899976e-10 1.000000e+00
[29,] 9.999283e-01 1.433011e-04 7.165055e-05
[30,] 4.627348e-04 9.254696e-04 9.995373e-01
[31,] 4.649493e-04 9.298986e-04 9.995351e-01
[32,] 1.532589e-09 3.065177e-09 1.000000e+00
[33,] 7.750845e-19 1.550169e-18 1.000000e+00
[34,] 9.991612e-01 1.677613e-03 8.388064e-04
[35,] 9.999999e-01 1.470475e-07 7.352377e-08
[36,] 6.947868e-01 6.104264e-01 3.052132e-01
[37,] 2.085514e-09 4.171027e-09 1.000000e+00
[38,] 1.000000e+00 1.941960e-09 9.709798e-10
[39,] 4.018911e-02 8.037821e-02 9.598109e-01
[40,] 2.263258e-02 4.526515e-02 9.773674e-01
[41,] 6.300508e-01 7.398983e-01 3.699492e-01
[42,] 1.000000e+00 4.567560e-08 2.283780e-08
[43,] 9.999891e-01 2.183910e-05 1.091955e-05
[44,] 4.687507e-51 9.375014e-51 1.000000e+00
[45,] 9.003589e-01 1.992822e-01 9.964109e-02
[46,] 8.868937e-02 1.773787e-01 9.113106e-01
[47,] 1.000000e+00 3.686205e-12 1.843102e-12
[48,] 9.992614e-01 1.477147e-03 7.385735e-04
[49,] 1.000000e+00 7.291564e-09 3.645782e-09
[50,] 9.999993e-01 1.491320e-06 7.456599e-07
[51,] 2.322313e-05 4.644626e-05 9.999768e-01
[52,] 9.966693e-01 6.661399e-03 3.330699e-03
[53,] 1.000000e+00 1.416078e-11 7.080388e-12
[54,] 2.837588e-04 5.675176e-04 9.997162e-01
[55,] 8.842260e-01 2.315480e-01 1.157740e-01
[56,] 9.712455e-01 5.750900e-02 2.875450e-02
[57,] 3.062316e-02 6.124633e-02 9.693768e-01
[58,] 2.214409e-03 4.428818e-03 9.977856e-01
[59,] 9.994716e-01 1.056769e-03 5.283845e-04
> postscript(file="/var/www/html/rcomp/tmp/1i8rg1258293315.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/20x0o1258293315.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/3ba5g1258293315.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/4pl0h1258293315.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/58dza1258293315.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
9.266735e-10 -2.679697e-10 6.787060e-11 -5.521909e-11 -1.248143e-11
6 7 8 9 10
-1.575271e-11 -1.440767e-11 -1.252278e-11 -6.852625e-12 -8.917928e-12
11 12 13 14 15
-4.219761e-12 -5.941961e-12 -1.729876e-11 -7.931544e-12 -7.981304e-12
16 17 18 19 20
-1.532081e-11 -6.380059e-12 -3.115322e-12 -6.552318e-12 -9.516933e-12
21 22 23 24 25
-6.275062e-12 -5.210670e-12 4.768375e-12 -9.578683e-12 -1.470485e-11
26 27 28 29 30
-1.738710e-11 -1.424848e-11 -1.810989e-11 -1.770397e-11 -2.025745e-11
31 32 33 34 35
-2.276741e-11 -2.924588e-12 -8.559174e-12 -5.150837e-12 -1.410442e-11
36 37 38 39 40
-1.832134e-11 -2.313649e-11 -2.298438e-11 -2.336569e-11 -2.727412e-11
41 42 43 44 45
-2.828732e-11 -2.693416e-11 -2.438971e-11 -1.719779e-11 -2.436292e-11
46 47 48 49 50
-1.849331e-11 -1.651628e-11 -1.900429e-11 -2.396285e-11 -2.341076e-11
51 52 53 54 55
-2.612305e-11 -2.089622e-11 -2.330777e-11 -2.475189e-11 -1.699797e-11
56 57 58 59 60
-5.445473e-14 9.853098e-12 -9.216613e-13 1.467349e-12 -2.134128e-12
61 62 63 64 65
-7.765273e-12 -1.252382e-11 -8.172916e-12 5.436850e-12 -1.452428e-12
66 67 68 69 70
-3.189818e-13 1.475805e-11 1.602106e-11 2.650989e-11 2.038970e-11
71 72
2.749469e-11 1.718406e-11
> postscript(file="/var/www/html/rcomp/tmp/6p4gp1258293315.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 9.266735e-10 NA
1 -2.679697e-10 9.266735e-10
2 6.787060e-11 -2.679697e-10
3 -5.521909e-11 6.787060e-11
4 -1.248143e-11 -5.521909e-11
5 -1.575271e-11 -1.248143e-11
6 -1.440767e-11 -1.575271e-11
7 -1.252278e-11 -1.440767e-11
8 -6.852625e-12 -1.252278e-11
9 -8.917928e-12 -6.852625e-12
10 -4.219761e-12 -8.917928e-12
11 -5.941961e-12 -4.219761e-12
12 -1.729876e-11 -5.941961e-12
13 -7.931544e-12 -1.729876e-11
14 -7.981304e-12 -7.931544e-12
15 -1.532081e-11 -7.981304e-12
16 -6.380059e-12 -1.532081e-11
17 -3.115322e-12 -6.380059e-12
18 -6.552318e-12 -3.115322e-12
19 -9.516933e-12 -6.552318e-12
20 -6.275062e-12 -9.516933e-12
21 -5.210670e-12 -6.275062e-12
22 4.768375e-12 -5.210670e-12
23 -9.578683e-12 4.768375e-12
24 -1.470485e-11 -9.578683e-12
25 -1.738710e-11 -1.470485e-11
26 -1.424848e-11 -1.738710e-11
27 -1.810989e-11 -1.424848e-11
28 -1.770397e-11 -1.810989e-11
29 -2.025745e-11 -1.770397e-11
30 -2.276741e-11 -2.025745e-11
31 -2.924588e-12 -2.276741e-11
32 -8.559174e-12 -2.924588e-12
33 -5.150837e-12 -8.559174e-12
34 -1.410442e-11 -5.150837e-12
35 -1.832134e-11 -1.410442e-11
36 -2.313649e-11 -1.832134e-11
37 -2.298438e-11 -2.313649e-11
38 -2.336569e-11 -2.298438e-11
39 -2.727412e-11 -2.336569e-11
40 -2.828732e-11 -2.727412e-11
41 -2.693416e-11 -2.828732e-11
42 -2.438971e-11 -2.693416e-11
43 -1.719779e-11 -2.438971e-11
44 -2.436292e-11 -1.719779e-11
45 -1.849331e-11 -2.436292e-11
46 -1.651628e-11 -1.849331e-11
47 -1.900429e-11 -1.651628e-11
48 -2.396285e-11 -1.900429e-11
49 -2.341076e-11 -2.396285e-11
50 -2.612305e-11 -2.341076e-11
51 -2.089622e-11 -2.612305e-11
52 -2.330777e-11 -2.089622e-11
53 -2.475189e-11 -2.330777e-11
54 -1.699797e-11 -2.475189e-11
55 -5.445473e-14 -1.699797e-11
56 9.853098e-12 -5.445473e-14
57 -9.216613e-13 9.853098e-12
58 1.467349e-12 -9.216613e-13
59 -2.134128e-12 1.467349e-12
60 -7.765273e-12 -2.134128e-12
61 -1.252382e-11 -7.765273e-12
62 -8.172916e-12 -1.252382e-11
63 5.436850e-12 -8.172916e-12
64 -1.452428e-12 5.436850e-12
65 -3.189818e-13 -1.452428e-12
66 1.475805e-11 -3.189818e-13
67 1.602106e-11 1.475805e-11
68 2.650989e-11 1.602106e-11
69 2.038970e-11 2.650989e-11
70 2.749469e-11 2.038970e-11
71 1.718406e-11 2.749469e-11
72 NA 1.718406e-11
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.679697e-10 9.266735e-10
[2,] 6.787060e-11 -2.679697e-10
[3,] -5.521909e-11 6.787060e-11
[4,] -1.248143e-11 -5.521909e-11
[5,] -1.575271e-11 -1.248143e-11
[6,] -1.440767e-11 -1.575271e-11
[7,] -1.252278e-11 -1.440767e-11
[8,] -6.852625e-12 -1.252278e-11
[9,] -8.917928e-12 -6.852625e-12
[10,] -4.219761e-12 -8.917928e-12
[11,] -5.941961e-12 -4.219761e-12
[12,] -1.729876e-11 -5.941961e-12
[13,] -7.931544e-12 -1.729876e-11
[14,] -7.981304e-12 -7.931544e-12
[15,] -1.532081e-11 -7.981304e-12
[16,] -6.380059e-12 -1.532081e-11
[17,] -3.115322e-12 -6.380059e-12
[18,] -6.552318e-12 -3.115322e-12
[19,] -9.516933e-12 -6.552318e-12
[20,] -6.275062e-12 -9.516933e-12
[21,] -5.210670e-12 -6.275062e-12
[22,] 4.768375e-12 -5.210670e-12
[23,] -9.578683e-12 4.768375e-12
[24,] -1.470485e-11 -9.578683e-12
[25,] -1.738710e-11 -1.470485e-11
[26,] -1.424848e-11 -1.738710e-11
[27,] -1.810989e-11 -1.424848e-11
[28,] -1.770397e-11 -1.810989e-11
[29,] -2.025745e-11 -1.770397e-11
[30,] -2.276741e-11 -2.025745e-11
[31,] -2.924588e-12 -2.276741e-11
[32,] -8.559174e-12 -2.924588e-12
[33,] -5.150837e-12 -8.559174e-12
[34,] -1.410442e-11 -5.150837e-12
[35,] -1.832134e-11 -1.410442e-11
[36,] -2.313649e-11 -1.832134e-11
[37,] -2.298438e-11 -2.313649e-11
[38,] -2.336569e-11 -2.298438e-11
[39,] -2.727412e-11 -2.336569e-11
[40,] -2.828732e-11 -2.727412e-11
[41,] -2.693416e-11 -2.828732e-11
[42,] -2.438971e-11 -2.693416e-11
[43,] -1.719779e-11 -2.438971e-11
[44,] -2.436292e-11 -1.719779e-11
[45,] -1.849331e-11 -2.436292e-11
[46,] -1.651628e-11 -1.849331e-11
[47,] -1.900429e-11 -1.651628e-11
[48,] -2.396285e-11 -1.900429e-11
[49,] -2.341076e-11 -2.396285e-11
[50,] -2.612305e-11 -2.341076e-11
[51,] -2.089622e-11 -2.612305e-11
[52,] -2.330777e-11 -2.089622e-11
[53,] -2.475189e-11 -2.330777e-11
[54,] -1.699797e-11 -2.475189e-11
[55,] -5.445473e-14 -1.699797e-11
[56,] 9.853098e-12 -5.445473e-14
[57,] -9.216613e-13 9.853098e-12
[58,] 1.467349e-12 -9.216613e-13
[59,] -2.134128e-12 1.467349e-12
[60,] -7.765273e-12 -2.134128e-12
[61,] -1.252382e-11 -7.765273e-12
[62,] -8.172916e-12 -1.252382e-11
[63,] 5.436850e-12 -8.172916e-12
[64,] -1.452428e-12 5.436850e-12
[65,] -3.189818e-13 -1.452428e-12
[66,] 1.475805e-11 -3.189818e-13
[67,] 1.602106e-11 1.475805e-11
[68,] 2.650989e-11 1.602106e-11
[69,] 2.038970e-11 2.650989e-11
[70,] 2.749469e-11 2.038970e-11
[71,] 1.718406e-11 2.749469e-11
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.679697e-10 9.266735e-10
2 6.787060e-11 -2.679697e-10
3 -5.521909e-11 6.787060e-11
4 -1.248143e-11 -5.521909e-11
5 -1.575271e-11 -1.248143e-11
6 -1.440767e-11 -1.575271e-11
7 -1.252278e-11 -1.440767e-11
8 -6.852625e-12 -1.252278e-11
9 -8.917928e-12 -6.852625e-12
10 -4.219761e-12 -8.917928e-12
11 -5.941961e-12 -4.219761e-12
12 -1.729876e-11 -5.941961e-12
13 -7.931544e-12 -1.729876e-11
14 -7.981304e-12 -7.931544e-12
15 -1.532081e-11 -7.981304e-12
16 -6.380059e-12 -1.532081e-11
17 -3.115322e-12 -6.380059e-12
18 -6.552318e-12 -3.115322e-12
19 -9.516933e-12 -6.552318e-12
20 -6.275062e-12 -9.516933e-12
21 -5.210670e-12 -6.275062e-12
22 4.768375e-12 -5.210670e-12
23 -9.578683e-12 4.768375e-12
24 -1.470485e-11 -9.578683e-12
25 -1.738710e-11 -1.470485e-11
26 -1.424848e-11 -1.738710e-11
27 -1.810989e-11 -1.424848e-11
28 -1.770397e-11 -1.810989e-11
29 -2.025745e-11 -1.770397e-11
30 -2.276741e-11 -2.025745e-11
31 -2.924588e-12 -2.276741e-11
32 -8.559174e-12 -2.924588e-12
33 -5.150837e-12 -8.559174e-12
34 -1.410442e-11 -5.150837e-12
35 -1.832134e-11 -1.410442e-11
36 -2.313649e-11 -1.832134e-11
37 -2.298438e-11 -2.313649e-11
38 -2.336569e-11 -2.298438e-11
39 -2.727412e-11 -2.336569e-11
40 -2.828732e-11 -2.727412e-11
41 -2.693416e-11 -2.828732e-11
42 -2.438971e-11 -2.693416e-11
43 -1.719779e-11 -2.438971e-11
44 -2.436292e-11 -1.719779e-11
45 -1.849331e-11 -2.436292e-11
46 -1.651628e-11 -1.849331e-11
47 -1.900429e-11 -1.651628e-11
48 -2.396285e-11 -1.900429e-11
49 -2.341076e-11 -2.396285e-11
50 -2.612305e-11 -2.341076e-11
51 -2.089622e-11 -2.612305e-11
52 -2.330777e-11 -2.089622e-11
53 -2.475189e-11 -2.330777e-11
54 -1.699797e-11 -2.475189e-11
55 -5.445473e-14 -1.699797e-11
56 9.853098e-12 -5.445473e-14
57 -9.216613e-13 9.853098e-12
58 1.467349e-12 -9.216613e-13
59 -2.134128e-12 1.467349e-12
60 -7.765273e-12 -2.134128e-12
61 -1.252382e-11 -7.765273e-12
62 -8.172916e-12 -1.252382e-11
63 5.436850e-12 -8.172916e-12
64 -1.452428e-12 5.436850e-12
65 -3.189818e-13 -1.452428e-12
66 1.475805e-11 -3.189818e-13
67 1.602106e-11 1.475805e-11
68 2.650989e-11 1.602106e-11
69 2.038970e-11 2.650989e-11
70 2.749469e-11 2.038970e-11
71 1.718406e-11 2.749469e-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/7db7t1258293315.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/8f2zy1258293315.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/9r6cp1258293315.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/10amod1258293315.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/11qa4n1258293315.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/12aw1g1258293315.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/134x2r1258293315.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/14w5w81258293315.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/15121h1258293315.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/16tzfb1258293315.tab")
+ }
>
> system("convert tmp/1i8rg1258293315.ps tmp/1i8rg1258293315.png")
> system("convert tmp/20x0o1258293315.ps tmp/20x0o1258293315.png")
> system("convert tmp/3ba5g1258293315.ps tmp/3ba5g1258293315.png")
> system("convert tmp/4pl0h1258293315.ps tmp/4pl0h1258293315.png")
> system("convert tmp/58dza1258293315.ps tmp/58dza1258293315.png")
> system("convert tmp/6p4gp1258293315.ps tmp/6p4gp1258293315.png")
> system("convert tmp/7db7t1258293315.ps tmp/7db7t1258293315.png")
> system("convert tmp/8f2zy1258293315.ps tmp/8f2zy1258293315.png")
> system("convert tmp/9r6cp1258293315.ps tmp/9r6cp1258293315.png")
> system("convert tmp/10amod1258293315.ps tmp/10amod1258293315.png")
>
>
> proc.time()
user system elapsed
2.572 1.548 3.486