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(344744
+ ,492865
+ ,338653
+ ,480961
+ ,327532
+ ,461935
+ ,326225
+ ,456608
+ ,318672
+ ,441977
+ ,317756
+ ,439148
+ ,337302
+ ,488180
+ ,349420
+ ,520564
+ ,336923
+ ,501492
+ ,330758
+ ,485025
+ ,321002
+ ,464196
+ ,320820
+ ,460170
+ ,327032
+ ,467037
+ ,324047
+ ,460070
+ ,316735
+ ,447988
+ ,315710
+ ,442867
+ ,313427
+ ,436087
+ ,310527
+ ,431328
+ ,330962
+ ,484015
+ ,339015
+ ,509673
+ ,341332
+ ,512927
+ ,339092
+ ,502831
+ ,323308
+ ,470984
+ ,325849
+ ,471067
+ ,330675
+ ,476049
+ ,332225
+ ,474605
+ ,331735
+ ,470439
+ ,328047
+ ,461251
+ ,326165
+ ,454724
+ ,327081
+ ,455626
+ ,346764
+ ,516847
+ ,344190
+ ,525192
+ ,343333
+ ,522975
+ ,345777
+ ,518585
+ ,344094
+ ,509239
+ ,348609
+ ,512238
+ ,354846
+ ,519164
+ ,356427
+ ,517009
+ ,353467
+ ,509933
+ ,355996
+ ,509127
+ ,352487
+ ,500857
+ ,355178
+ ,506971
+ ,374556
+ ,569323
+ ,375021
+ ,579714
+ ,375787
+ ,577992
+ ,372720
+ ,565464
+ ,364431
+ ,547344
+ ,370490
+ ,554788
+ ,376974
+ ,562325
+ ,377632
+ ,560854
+ ,378205
+ ,555332
+ ,370861
+ ,543599
+ ,369167
+ ,536662
+ ,371551
+ ,542722
+ ,382842
+ ,593530
+ ,381903
+ ,610763
+ ,384502
+ ,612613
+ ,392058
+ ,611324
+ ,384359
+ ,594167
+ ,388884
+ ,595454
+ ,386586
+ ,590865
+ ,387495
+ ,589379
+ ,385705
+ ,584428
+ ,378670
+ ,573100
+ ,377367
+ ,567456
+ ,376911
+ ,569028
+ ,389827
+ ,620735
+ ,387820
+ ,628884
+ ,387267
+ ,628232
+ ,380575
+ ,612117
+ ,372402
+ ,595404
+ ,376740
+ ,597141)
+ ,dim=c(2
+ ,72)
+ ,dimnames=list(c('Y'
+ ,'X')
+ ,1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Y','X'),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 = 'Linear Trend'
> par2 = 'Include Monthly 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.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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 344744 492865 1 0 0 0 0 0 0 0 0 0 0 1
2 338653 480961 0 1 0 0 0 0 0 0 0 0 0 2
3 327532 461935 0 0 1 0 0 0 0 0 0 0 0 3
4 326225 456608 0 0 0 1 0 0 0 0 0 0 0 4
5 318672 441977 0 0 0 0 1 0 0 0 0 0 0 5
6 317756 439148 0 0 0 0 0 1 0 0 0 0 0 6
7 337302 488180 0 0 0 0 0 0 1 0 0 0 0 7
8 349420 520564 0 0 0 0 0 0 0 1 0 0 0 8
9 336923 501492 0 0 0 0 0 0 0 0 1 0 0 9
10 330758 485025 0 0 0 0 0 0 0 0 0 1 0 10
11 321002 464196 0 0 0 0 0 0 0 0 0 0 1 11
12 320820 460170 0 0 0 0 0 0 0 0 0 0 0 12
13 327032 467037 1 0 0 0 0 0 0 0 0 0 0 13
14 324047 460070 0 1 0 0 0 0 0 0 0 0 0 14
15 316735 447988 0 0 1 0 0 0 0 0 0 0 0 15
16 315710 442867 0 0 0 1 0 0 0 0 0 0 0 16
17 313427 436087 0 0 0 0 1 0 0 0 0 0 0 17
18 310527 431328 0 0 0 0 0 1 0 0 0 0 0 18
19 330962 484015 0 0 0 0 0 0 1 0 0 0 0 19
20 339015 509673 0 0 0 0 0 0 0 1 0 0 0 20
21 341332 512927 0 0 0 0 0 0 0 0 1 0 0 21
22 339092 502831 0 0 0 0 0 0 0 0 0 1 0 22
23 323308 470984 0 0 0 0 0 0 0 0 0 0 1 23
24 325849 471067 0 0 0 0 0 0 0 0 0 0 0 24
25 330675 476049 1 0 0 0 0 0 0 0 0 0 0 25
26 332225 474605 0 1 0 0 0 0 0 0 0 0 0 26
27 331735 470439 0 0 1 0 0 0 0 0 0 0 0 27
28 328047 461251 0 0 0 1 0 0 0 0 0 0 0 28
29 326165 454724 0 0 0 0 1 0 0 0 0 0 0 29
30 327081 455626 0 0 0 0 0 1 0 0 0 0 0 30
31 346764 516847 0 0 0 0 0 0 1 0 0 0 0 31
32 344190 525192 0 0 0 0 0 0 0 1 0 0 0 32
33 343333 522975 0 0 0 0 0 0 0 0 1 0 0 33
34 345777 518585 0 0 0 0 0 0 0 0 0 1 0 34
35 344094 509239 0 0 0 0 0 0 0 0 0 0 1 35
36 348609 512238 0 0 0 0 0 0 0 0 0 0 0 36
37 354846 519164 1 0 0 0 0 0 0 0 0 0 0 37
38 356427 517009 0 1 0 0 0 0 0 0 0 0 0 38
39 353467 509933 0 0 1 0 0 0 0 0 0 0 0 39
40 355996 509127 0 0 0 1 0 0 0 0 0 0 0 40
41 352487 500857 0 0 0 0 1 0 0 0 0 0 0 41
42 355178 506971 0 0 0 0 0 1 0 0 0 0 0 42
43 374556 569323 0 0 0 0 0 0 1 0 0 0 0 43
44 375021 579714 0 0 0 0 0 0 0 1 0 0 0 44
45 375787 577992 0 0 0 0 0 0 0 0 1 0 0 45
46 372720 565464 0 0 0 0 0 0 0 0 0 1 0 46
47 364431 547344 0 0 0 0 0 0 0 0 0 0 1 47
48 370490 554788 0 0 0 0 0 0 0 0 0 0 0 48
49 376974 562325 1 0 0 0 0 0 0 0 0 0 0 49
50 377632 560854 0 1 0 0 0 0 0 0 0 0 0 50
51 378205 555332 0 0 1 0 0 0 0 0 0 0 0 51
52 370861 543599 0 0 0 1 0 0 0 0 0 0 0 52
53 369167 536662 0 0 0 0 1 0 0 0 0 0 0 53
54 371551 542722 0 0 0 0 0 1 0 0 0 0 0 54
55 382842 593530 0 0 0 0 0 0 1 0 0 0 0 55
56 381903 610763 0 0 0 0 0 0 0 1 0 0 0 56
57 384502 612613 0 0 0 0 0 0 0 0 1 0 0 57
58 392058 611324 0 0 0 0 0 0 0 0 0 1 0 58
59 384359 594167 0 0 0 0 0 0 0 0 0 0 1 59
60 388884 595454 0 0 0 0 0 0 0 0 0 0 0 60
61 386586 590865 1 0 0 0 0 0 0 0 0 0 0 61
62 387495 589379 0 1 0 0 0 0 0 0 0 0 0 62
63 385705 584428 0 0 1 0 0 0 0 0 0 0 0 63
64 378670 573100 0 0 0 1 0 0 0 0 0 0 0 64
65 377367 567456 0 0 0 0 1 0 0 0 0 0 0 65
66 376911 569028 0 0 0 0 0 1 0 0 0 0 0 66
67 389827 620735 0 0 0 0 0 0 1 0 0 0 0 67
68 387820 628884 0 0 0 0 0 0 0 1 0 0 0 68
69 387267 628232 0 0 0 0 0 0 0 0 1 0 0 69
70 380575 612117 0 0 0 0 0 0 0 0 0 1 0 70
71 372402 595404 0 0 0 0 0 0 0 0 0 0 1 71
72 376740 597141 0 0 0 0 0 0 0 0 0 0 0 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
7.252e+04 5.485e-01 3.445e+03 5.252e+03 6.444e+03 7.656e+03
M5 M6 M7 M8 M9 M10
9.291e+03 9.146e+03 -3.397e+03 -1.000e+04 -9.464e+03 -5.046e+03
M11 t
-2.975e+03 -2.132e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8549.3 -2406.2 474.9 2750.6 5547.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.252e+04 1.276e+04 5.684 4.49e-07 ***
X 5.484e-01 2.910e-02 18.850 < 2e-16 ***
M1 3.445e+03 2.220e+03 1.552 0.126137
M2 5.252e+03 2.196e+03 2.392 0.020023 *
M3 6.444e+03 2.195e+03 2.936 0.004765 **
M4 7.656e+03 2.232e+03 3.429 0.001120 **
M5 9.291e+03 2.311e+03 4.020 0.000170 ***
M6 9.146e+03 2.321e+03 3.940 0.000221 ***
M7 -3.397e+03 2.303e+03 -1.475 0.145553
M8 -1.000e+04 2.472e+03 -4.047 0.000156 ***
M9 -9.464e+03 2.401e+03 -3.941 0.000220 ***
M10 -5.046e+03 2.273e+03 -2.220 0.030320 *
M11 -2.975e+03 2.180e+03 -1.365 0.177576
t -2.132e+02 7.109e+01 -3.000 0.003979 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3775 on 58 degrees of freedom
Multiple R-squared: 0.9809, Adjusted R-squared: 0.9766
F-statistic: 229.4 on 13 and 58 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,] 0.0007152446 0.0014304893 0.99928476
[2,] 0.0002564466 0.0005128933 0.99974355
[3,] 0.0013660457 0.0027320913 0.99863395
[4,] 0.0013525958 0.0027051916 0.99864740
[5,] 0.0007330761 0.0014661522 0.99926692
[6,] 0.0006105437 0.0012210874 0.99938946
[7,] 0.0003516634 0.0007033269 0.99964834
[8,] 0.0002275665 0.0004551330 0.99977243
[9,] 0.0001414151 0.0002828301 0.99985858
[10,] 0.0005285929 0.0010571858 0.99947141
[11,] 0.0092599456 0.0185198912 0.99074005
[12,] 0.0205767321 0.0411534643 0.97942327
[13,] 0.0392326677 0.0784653354 0.96076733
[14,] 0.0323537397 0.0647074794 0.96764626
[15,] 0.3635500138 0.7271000276 0.63644999
[16,] 0.5085635353 0.9828729293 0.49143646
[17,] 0.5325835983 0.9348328034 0.46741640
[18,] 0.4997291167 0.9994582333 0.50027088
[19,] 0.4475678029 0.8951356059 0.55243220
[20,] 0.4078607737 0.8157215473 0.59213923
[21,] 0.3919106589 0.7838213178 0.60808934
[22,] 0.3496232235 0.6992464471 0.65037678
[23,] 0.3935777434 0.7871554868 0.60642226
[24,] 0.3799867708 0.7599735416 0.62001323
[25,] 0.3854149095 0.7708298191 0.61458509
[26,] 0.3916225359 0.7832450718 0.60837746
[27,] 0.7093903124 0.5812193753 0.29060969
[28,] 0.6789772542 0.6420454916 0.32102275
[29,] 0.6420277787 0.7159444426 0.35797222
[30,] 0.5685484090 0.8629031821 0.43145159
[31,] 0.5481721971 0.9036556059 0.45182780
[32,] 0.4686975274 0.9373950547 0.53130247
[33,] 0.4040391070 0.8080782140 0.59596089
[34,] 0.3582171828 0.7164343656 0.64178282
[35,] 0.2799746791 0.5599493582 0.72002532
[36,] 0.2156278741 0.4312557481 0.78437213
[37,] 0.2013749105 0.4027498209 0.79862509
[38,] 0.2568576226 0.5137152452 0.74314238
[39,] 0.9737384576 0.0525230848 0.02626154
> postscript(file="/var/www/html/rcomp/tmp/1nhnj1258485883.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/27tc61258485883.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/36kfd1258485883.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/4b3c61258485883.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/5yph51258485883.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 6
-1316.06650 -2472.86366 -4137.48025 -3521.09729 -4472.09023 -3477.66820
7 8 9 10 11 12
1932.87604 3108.75901 746.23991 -592.14133 -782.11673 -1518.13091
13 14 15 16 17 18
-2303.74508 -3062.25731 -4726.33523 -3940.93369 -3927.83531 -3858.89796
19 20 21 22 23 24
436.04862 1235.83004 1442.53741 534.95872 359.84432 93.23459
25 26 27 28 29 30
-1044.54467 -297.16576 519.19817 872.16022 1147.49986 1927.64176
31 32 33 34 35 36
790.08591 458.24332 491.53995 1138.48541 2723.61812 2831.71788
37 38 39 40 41 42
2038.74496 3207.07433 3149.43803 5122.26259 4726.55673 4423.15883
43 44 45 46 47 48
2360.30204 3945.32352 5330.13565 4929.39595 4720.65994 3934.88376
49 50 51 52 53 54
3053.80572 2923.99288 5547.05980 3639.83609 4328.04168 3747.26027
55 56 57 58 59 60
-71.24896 -3642.74654 -2384.01042 1674.18064 1527.28387 2584.33608
61 62 63 64 65 66
-428.19443 -298.78047 -351.88052 -2172.22791 -1802.17274 -2761.49470
67 68 69 70 71 72
-5448.06365 -5105.40935 -5626.44250 -7684.87939 -8549.28952 -7926.04140
> postscript(file="/var/www/html/rcomp/tmp/6x9401258485883.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 -1316.06650 NA
1 -2472.86366 -1316.06650
2 -4137.48025 -2472.86366
3 -3521.09729 -4137.48025
4 -4472.09023 -3521.09729
5 -3477.66820 -4472.09023
6 1932.87604 -3477.66820
7 3108.75901 1932.87604
8 746.23991 3108.75901
9 -592.14133 746.23991
10 -782.11673 -592.14133
11 -1518.13091 -782.11673
12 -2303.74508 -1518.13091
13 -3062.25731 -2303.74508
14 -4726.33523 -3062.25731
15 -3940.93369 -4726.33523
16 -3927.83531 -3940.93369
17 -3858.89796 -3927.83531
18 436.04862 -3858.89796
19 1235.83004 436.04862
20 1442.53741 1235.83004
21 534.95872 1442.53741
22 359.84432 534.95872
23 93.23459 359.84432
24 -1044.54467 93.23459
25 -297.16576 -1044.54467
26 519.19817 -297.16576
27 872.16022 519.19817
28 1147.49986 872.16022
29 1927.64176 1147.49986
30 790.08591 1927.64176
31 458.24332 790.08591
32 491.53995 458.24332
33 1138.48541 491.53995
34 2723.61812 1138.48541
35 2831.71788 2723.61812
36 2038.74496 2831.71788
37 3207.07433 2038.74496
38 3149.43803 3207.07433
39 5122.26259 3149.43803
40 4726.55673 5122.26259
41 4423.15883 4726.55673
42 2360.30204 4423.15883
43 3945.32352 2360.30204
44 5330.13565 3945.32352
45 4929.39595 5330.13565
46 4720.65994 4929.39595
47 3934.88376 4720.65994
48 3053.80572 3934.88376
49 2923.99288 3053.80572
50 5547.05980 2923.99288
51 3639.83609 5547.05980
52 4328.04168 3639.83609
53 3747.26027 4328.04168
54 -71.24896 3747.26027
55 -3642.74654 -71.24896
56 -2384.01042 -3642.74654
57 1674.18064 -2384.01042
58 1527.28387 1674.18064
59 2584.33608 1527.28387
60 -428.19443 2584.33608
61 -298.78047 -428.19443
62 -351.88052 -298.78047
63 -2172.22791 -351.88052
64 -1802.17274 -2172.22791
65 -2761.49470 -1802.17274
66 -5448.06365 -2761.49470
67 -5105.40935 -5448.06365
68 -5626.44250 -5105.40935
69 -7684.87939 -5626.44250
70 -8549.28952 -7684.87939
71 -7926.04140 -8549.28952
72 NA -7926.04140
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2472.86366 -1316.06650
[2,] -4137.48025 -2472.86366
[3,] -3521.09729 -4137.48025
[4,] -4472.09023 -3521.09729
[5,] -3477.66820 -4472.09023
[6,] 1932.87604 -3477.66820
[7,] 3108.75901 1932.87604
[8,] 746.23991 3108.75901
[9,] -592.14133 746.23991
[10,] -782.11673 -592.14133
[11,] -1518.13091 -782.11673
[12,] -2303.74508 -1518.13091
[13,] -3062.25731 -2303.74508
[14,] -4726.33523 -3062.25731
[15,] -3940.93369 -4726.33523
[16,] -3927.83531 -3940.93369
[17,] -3858.89796 -3927.83531
[18,] 436.04862 -3858.89796
[19,] 1235.83004 436.04862
[20,] 1442.53741 1235.83004
[21,] 534.95872 1442.53741
[22,] 359.84432 534.95872
[23,] 93.23459 359.84432
[24,] -1044.54467 93.23459
[25,] -297.16576 -1044.54467
[26,] 519.19817 -297.16576
[27,] 872.16022 519.19817
[28,] 1147.49986 872.16022
[29,] 1927.64176 1147.49986
[30,] 790.08591 1927.64176
[31,] 458.24332 790.08591
[32,] 491.53995 458.24332
[33,] 1138.48541 491.53995
[34,] 2723.61812 1138.48541
[35,] 2831.71788 2723.61812
[36,] 2038.74496 2831.71788
[37,] 3207.07433 2038.74496
[38,] 3149.43803 3207.07433
[39,] 5122.26259 3149.43803
[40,] 4726.55673 5122.26259
[41,] 4423.15883 4726.55673
[42,] 2360.30204 4423.15883
[43,] 3945.32352 2360.30204
[44,] 5330.13565 3945.32352
[45,] 4929.39595 5330.13565
[46,] 4720.65994 4929.39595
[47,] 3934.88376 4720.65994
[48,] 3053.80572 3934.88376
[49,] 2923.99288 3053.80572
[50,] 5547.05980 2923.99288
[51,] 3639.83609 5547.05980
[52,] 4328.04168 3639.83609
[53,] 3747.26027 4328.04168
[54,] -71.24896 3747.26027
[55,] -3642.74654 -71.24896
[56,] -2384.01042 -3642.74654
[57,] 1674.18064 -2384.01042
[58,] 1527.28387 1674.18064
[59,] 2584.33608 1527.28387
[60,] -428.19443 2584.33608
[61,] -298.78047 -428.19443
[62,] -351.88052 -298.78047
[63,] -2172.22791 -351.88052
[64,] -1802.17274 -2172.22791
[65,] -2761.49470 -1802.17274
[66,] -5448.06365 -2761.49470
[67,] -5105.40935 -5448.06365
[68,] -5626.44250 -5105.40935
[69,] -7684.87939 -5626.44250
[70,] -8549.28952 -7684.87939
[71,] -7926.04140 -8549.28952
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2472.86366 -1316.06650
2 -4137.48025 -2472.86366
3 -3521.09729 -4137.48025
4 -4472.09023 -3521.09729
5 -3477.66820 -4472.09023
6 1932.87604 -3477.66820
7 3108.75901 1932.87604
8 746.23991 3108.75901
9 -592.14133 746.23991
10 -782.11673 -592.14133
11 -1518.13091 -782.11673
12 -2303.74508 -1518.13091
13 -3062.25731 -2303.74508
14 -4726.33523 -3062.25731
15 -3940.93369 -4726.33523
16 -3927.83531 -3940.93369
17 -3858.89796 -3927.83531
18 436.04862 -3858.89796
19 1235.83004 436.04862
20 1442.53741 1235.83004
21 534.95872 1442.53741
22 359.84432 534.95872
23 93.23459 359.84432
24 -1044.54467 93.23459
25 -297.16576 -1044.54467
26 519.19817 -297.16576
27 872.16022 519.19817
28 1147.49986 872.16022
29 1927.64176 1147.49986
30 790.08591 1927.64176
31 458.24332 790.08591
32 491.53995 458.24332
33 1138.48541 491.53995
34 2723.61812 1138.48541
35 2831.71788 2723.61812
36 2038.74496 2831.71788
37 3207.07433 2038.74496
38 3149.43803 3207.07433
39 5122.26259 3149.43803
40 4726.55673 5122.26259
41 4423.15883 4726.55673
42 2360.30204 4423.15883
43 3945.32352 2360.30204
44 5330.13565 3945.32352
45 4929.39595 5330.13565
46 4720.65994 4929.39595
47 3934.88376 4720.65994
48 3053.80572 3934.88376
49 2923.99288 3053.80572
50 5547.05980 2923.99288
51 3639.83609 5547.05980
52 4328.04168 3639.83609
53 3747.26027 4328.04168
54 -71.24896 3747.26027
55 -3642.74654 -71.24896
56 -2384.01042 -3642.74654
57 1674.18064 -2384.01042
58 1527.28387 1674.18064
59 2584.33608 1527.28387
60 -428.19443 2584.33608
61 -298.78047 -428.19443
62 -351.88052 -298.78047
63 -2172.22791 -351.88052
64 -1802.17274 -2172.22791
65 -2761.49470 -1802.17274
66 -5448.06365 -2761.49470
67 -5105.40935 -5448.06365
68 -5626.44250 -5105.40935
69 -7684.87939 -5626.44250
70 -8549.28952 -7684.87939
71 -7926.04140 -8549.28952
> 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/787ra1258485883.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/8le3i1258485883.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/9uj2w1258485883.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/10ypdc1258485883.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/11ze1q1258485883.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/12bzkn1258485883.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/13w2o21258485883.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/14140o1258485883.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/1522c41258485883.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/1694r51258485883.tab")
+ }
>
> system("convert tmp/1nhnj1258485883.ps tmp/1nhnj1258485883.png")
> system("convert tmp/27tc61258485883.ps tmp/27tc61258485883.png")
> system("convert tmp/36kfd1258485883.ps tmp/36kfd1258485883.png")
> system("convert tmp/4b3c61258485883.ps tmp/4b3c61258485883.png")
> system("convert tmp/5yph51258485883.ps tmp/5yph51258485883.png")
> system("convert tmp/6x9401258485883.ps tmp/6x9401258485883.png")
> system("convert tmp/787ra1258485883.ps tmp/787ra1258485883.png")
> system("convert tmp/8le3i1258485883.ps tmp/8le3i1258485883.png")
> system("convert tmp/9uj2w1258485883.ps tmp/9uj2w1258485883.png")
> system("convert tmp/10ypdc1258485883.ps tmp/10ypdc1258485883.png")
>
>
> proc.time()
user system elapsed
2.511 1.637 2.995