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 = 'No 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
1 344744 492865 1 0 0 0 0 0 0 0 0 0 0
2 338653 480961 0 1 0 0 0 0 0 0 0 0 0
3 327532 461935 0 0 1 0 0 0 0 0 0 0 0
4 326225 456608 0 0 0 1 0 0 0 0 0 0 0
5 318672 441977 0 0 0 0 1 0 0 0 0 0 0
6 317756 439148 0 0 0 0 0 1 0 0 0 0 0
7 337302 488180 0 0 0 0 0 0 1 0 0 0 0
8 349420 520564 0 0 0 0 0 0 0 1 0 0 0
9 336923 501492 0 0 0 0 0 0 0 0 1 0 0
10 330758 485025 0 0 0 0 0 0 0 0 0 1 0
11 321002 464196 0 0 0 0 0 0 0 0 0 0 1
12 320820 460170 0 0 0 0 0 0 0 0 0 0 0
13 327032 467037 1 0 0 0 0 0 0 0 0 0 0
14 324047 460070 0 1 0 0 0 0 0 0 0 0 0
15 316735 447988 0 0 1 0 0 0 0 0 0 0 0
16 315710 442867 0 0 0 1 0 0 0 0 0 0 0
17 313427 436087 0 0 0 0 1 0 0 0 0 0 0
18 310527 431328 0 0 0 0 0 1 0 0 0 0 0
19 330962 484015 0 0 0 0 0 0 1 0 0 0 0
20 339015 509673 0 0 0 0 0 0 0 1 0 0 0
21 341332 512927 0 0 0 0 0 0 0 0 1 0 0
22 339092 502831 0 0 0 0 0 0 0 0 0 1 0
23 323308 470984 0 0 0 0 0 0 0 0 0 0 1
24 325849 471067 0 0 0 0 0 0 0 0 0 0 0
25 330675 476049 1 0 0 0 0 0 0 0 0 0 0
26 332225 474605 0 1 0 0 0 0 0 0 0 0 0
27 331735 470439 0 0 1 0 0 0 0 0 0 0 0
28 328047 461251 0 0 0 1 0 0 0 0 0 0 0
29 326165 454724 0 0 0 0 1 0 0 0 0 0 0
30 327081 455626 0 0 0 0 0 1 0 0 0 0 0
31 346764 516847 0 0 0 0 0 0 1 0 0 0 0
32 344190 525192 0 0 0 0 0 0 0 1 0 0 0
33 343333 522975 0 0 0 0 0 0 0 0 1 0 0
34 345777 518585 0 0 0 0 0 0 0 0 0 1 0
35 344094 509239 0 0 0 0 0 0 0 0 0 0 1
36 348609 512238 0 0 0 0 0 0 0 0 0 0 0
37 354846 519164 1 0 0 0 0 0 0 0 0 0 0
38 356427 517009 0 1 0 0 0 0 0 0 0 0 0
39 353467 509933 0 0 1 0 0 0 0 0 0 0 0
40 355996 509127 0 0 0 1 0 0 0 0 0 0 0
41 352487 500857 0 0 0 0 1 0 0 0 0 0 0
42 355178 506971 0 0 0 0 0 1 0 0 0 0 0
43 374556 569323 0 0 0 0 0 0 1 0 0 0 0
44 375021 579714 0 0 0 0 0 0 0 1 0 0 0
45 375787 577992 0 0 0 0 0 0 0 0 1 0 0
46 372720 565464 0 0 0 0 0 0 0 0 0 1 0
47 364431 547344 0 0 0 0 0 0 0 0 0 0 1
48 370490 554788 0 0 0 0 0 0 0 0 0 0 0
49 376974 562325 1 0 0 0 0 0 0 0 0 0 0
50 377632 560854 0 1 0 0 0 0 0 0 0 0 0
51 378205 555332 0 0 1 0 0 0 0 0 0 0 0
52 370861 543599 0 0 0 1 0 0 0 0 0 0 0
53 369167 536662 0 0 0 0 1 0 0 0 0 0 0
54 371551 542722 0 0 0 0 0 1 0 0 0 0 0
55 382842 593530 0 0 0 0 0 0 1 0 0 0 0
56 381903 610763 0 0 0 0 0 0 0 1 0 0 0
57 384502 612613 0 0 0 0 0 0 0 0 1 0 0
58 392058 611324 0 0 0 0 0 0 0 0 0 1 0
59 384359 594167 0 0 0 0 0 0 0 0 0 0 1
60 388884 595454 0 0 0 0 0 0 0 0 0 0 0
61 386586 590865 1 0 0 0 0 0 0 0 0 0 0
62 387495 589379 0 1 0 0 0 0 0 0 0 0 0
63 385705 584428 0 0 1 0 0 0 0 0 0 0 0
64 378670 573100 0 0 0 1 0 0 0 0 0 0 0
65 377367 567456 0 0 0 0 1 0 0 0 0 0 0
66 376911 569028 0 0 0 0 0 1 0 0 0 0 0
67 389827 620735 0 0 0 0 0 0 1 0 0 0 0
68 387820 628884 0 0 0 0 0 0 0 1 0 0 0
69 387267 628232 0 0 0 0 0 0 0 0 1 0 0
70 380575 612117 0 0 0 0 0 0 0 0 0 1 0
71 372402 595404 0 0 0 0 0 0 0 0 0 0 1
72 376740 597141 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) X M1 M2 M3 M4
1.078e+05 4.653e-01 4.647e+03 5.889e+03 6.136e+03 6.532e+03
M5 M6 M7 M8 M9 M10
7.278e+03 7.017e+03 -1.199e+03 -6.603e+03 -6.534e+03 -3.173e+03
M11
-2.894e+03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-9529.4 -1856.9 -138.6 2995.6 6023.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.078e+05 5.297e+03 20.343 < 2e-16 ***
X 4.653e-01 9.469e-03 49.142 < 2e-16 ***
M1 4.647e+03 2.326e+03 1.997 0.05040 .
M2 5.889e+03 2.329e+03 2.529 0.01415 *
M3 6.136e+03 2.337e+03 2.626 0.01099 *
M4 6.532e+03 2.345e+03 2.785 0.00718 **
M5 7.278e+03 2.357e+03 3.088 0.00307 **
M6 7.017e+03 2.355e+03 2.980 0.00419 **
M7 -1.199e+03 2.326e+03 -0.515 0.60834
M8 -6.603e+03 2.341e+03 -2.821 0.00652 **
M9 -6.534e+03 2.337e+03 -2.795 0.00698 **
M10 -3.173e+03 2.329e+03 -1.362 0.17826
M11 -2.894e+03 2.323e+03 -1.246 0.21773
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4023 on 59 degrees of freedom
Multiple R-squared: 0.978, Adjusted R-squared: 0.9735
F-statistic: 218.2 on 12 and 59 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,] 3.048194e-03 6.096388e-03 0.9969518
[2,] 8.828138e-04 1.765628e-03 0.9991172
[3,] 5.898687e-04 1.179737e-03 0.9994101
[4,] 2.172193e-03 4.344385e-03 0.9978278
[5,] 1.377704e-03 2.755409e-03 0.9986223
[6,] 2.642682e-03 5.285364e-03 0.9973573
[7,] 3.736312e-03 7.472624e-03 0.9962637
[8,] 2.067868e-03 4.135737e-03 0.9979321
[9,] 1.167036e-03 2.334072e-03 0.9988330
[10,] 9.765610e-04 1.953122e-03 0.9990234
[11,] 6.273424e-04 1.254685e-03 0.9993727
[12,] 3.223250e-04 6.446501e-04 0.9996777
[13,] 1.599204e-04 3.198408e-04 0.9998401
[14,] 8.212407e-05 1.642481e-04 0.9999179
[15,] 4.219512e-05 8.439024e-05 0.9999578
[16,] 6.853166e-04 1.370633e-03 0.9993147
[17,] 2.558438e-03 5.116876e-03 0.9974416
[18,] 4.910404e-03 9.820807e-03 0.9950896
[19,] 5.150983e-03 1.030197e-02 0.9948490
[20,] 3.149216e-03 6.298432e-03 0.9968508
[21,] 2.044954e-03 4.089908e-03 0.9979550
[22,] 1.834876e-03 3.669753e-03 0.9981651
[23,] 1.582430e-03 3.164861e-03 0.9984176
[24,] 2.986613e-03 5.973227e-03 0.9970134
[25,] 2.413491e-03 4.826983e-03 0.9975865
[26,] 2.439827e-03 4.879655e-03 0.9975602
[27,] 2.129386e-03 4.258772e-03 0.9978706
[28,] 2.294876e-03 4.589751e-03 0.9977051
[29,] 1.535626e-03 3.071253e-03 0.9984644
[30,] 7.972863e-04 1.594573e-03 0.9992027
[31,] 3.688931e-04 7.377861e-04 0.9996311
[32,] 1.660643e-04 3.321286e-04 0.9998339
[33,] 7.372729e-05 1.474546e-04 0.9999263
[34,] 3.809175e-05 7.618351e-05 0.9999619
[35,] 1.958328e-05 3.916656e-05 0.9999804
[36,] 7.898054e-06 1.579611e-05 0.9999921
[37,] 2.879324e-06 5.758648e-06 0.9999971
[38,] 8.858298e-07 1.771660e-06 0.9999991
[39,] 2.545821e-07 5.091641e-07 0.9999997
[40,] 6.227498e-07 1.245500e-06 0.9999994
[41,] 1.058326e-05 2.116651e-05 0.9999894
> postscript(file="/var/www/html/rcomp/tmp/1buhs1258649687.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/2vahx1258649687.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/320ff1258649687.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/472sw1258649687.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/5ilvq1258649687.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 7
2987.9672 1194.0560 -1320.0968 -543.8570 -2035.0055 -1373.4834 3571.4524
8 9 10 11 12 13 14
6023.6412 2333.1651 469.5876 127.7753 -1074.7476 -2705.0691 -3690.3949
15 16 17 18 19 20 21
-5626.9122 -4664.5338 -4539.1159 -4963.4755 -830.3803 686.7275 1420.9304
22 23 24 25 26 27 28
517.6318 -724.9953 -1116.6260 -3255.7698 -2276.2030 -1074.4015 -882.4598
29 30 31 32 33 34 35
-473.7746 283.5406 -306.6294 -1359.9814 -1253.8690 -128.4335 2259.1828
36 37 38 39 40 41 42
2484.6023 851.8254 2193.2533 2279.2125 4787.6217 4380.4055 4487.3358
43 44 45 46 47 48 49
3065.8592 4099.4087 5598.1747 4999.5982 4864.1628 4565.1181 2895.0146
50 51 52 53 54 55 56
2995.1456 5890.9567 3611.2047 4398.6819 4223.7409 87.2218 -3467.1280
57 58 59 60 61 62 63
-1797.5789 2996.8178 3003.2539 4035.3462 -773.9683 -415.8570 -148.7587
64 65 66 67 68 69 70
-2307.9759 -1731.1914 -2657.6583 -5587.5237 -5982.6679 -6300.8224 -8855.2018
71 72
-9529.3795 -8893.6929
> postscript(file="/var/www/html/rcomp/tmp/6h5ld1258649687.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 2987.9672 NA
1 1194.0560 2987.9672
2 -1320.0968 1194.0560
3 -543.8570 -1320.0968
4 -2035.0055 -543.8570
5 -1373.4834 -2035.0055
6 3571.4524 -1373.4834
7 6023.6412 3571.4524
8 2333.1651 6023.6412
9 469.5876 2333.1651
10 127.7753 469.5876
11 -1074.7476 127.7753
12 -2705.0691 -1074.7476
13 -3690.3949 -2705.0691
14 -5626.9122 -3690.3949
15 -4664.5338 -5626.9122
16 -4539.1159 -4664.5338
17 -4963.4755 -4539.1159
18 -830.3803 -4963.4755
19 686.7275 -830.3803
20 1420.9304 686.7275
21 517.6318 1420.9304
22 -724.9953 517.6318
23 -1116.6260 -724.9953
24 -3255.7698 -1116.6260
25 -2276.2030 -3255.7698
26 -1074.4015 -2276.2030
27 -882.4598 -1074.4015
28 -473.7746 -882.4598
29 283.5406 -473.7746
30 -306.6294 283.5406
31 -1359.9814 -306.6294
32 -1253.8690 -1359.9814
33 -128.4335 -1253.8690
34 2259.1828 -128.4335
35 2484.6023 2259.1828
36 851.8254 2484.6023
37 2193.2533 851.8254
38 2279.2125 2193.2533
39 4787.6217 2279.2125
40 4380.4055 4787.6217
41 4487.3358 4380.4055
42 3065.8592 4487.3358
43 4099.4087 3065.8592
44 5598.1747 4099.4087
45 4999.5982 5598.1747
46 4864.1628 4999.5982
47 4565.1181 4864.1628
48 2895.0146 4565.1181
49 2995.1456 2895.0146
50 5890.9567 2995.1456
51 3611.2047 5890.9567
52 4398.6819 3611.2047
53 4223.7409 4398.6819
54 87.2218 4223.7409
55 -3467.1280 87.2218
56 -1797.5789 -3467.1280
57 2996.8178 -1797.5789
58 3003.2539 2996.8178
59 4035.3462 3003.2539
60 -773.9683 4035.3462
61 -415.8570 -773.9683
62 -148.7587 -415.8570
63 -2307.9759 -148.7587
64 -1731.1914 -2307.9759
65 -2657.6583 -1731.1914
66 -5587.5237 -2657.6583
67 -5982.6679 -5587.5237
68 -6300.8224 -5982.6679
69 -8855.2018 -6300.8224
70 -9529.3795 -8855.2018
71 -8893.6929 -9529.3795
72 NA -8893.6929
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1194.0560 2987.9672
[2,] -1320.0968 1194.0560
[3,] -543.8570 -1320.0968
[4,] -2035.0055 -543.8570
[5,] -1373.4834 -2035.0055
[6,] 3571.4524 -1373.4834
[7,] 6023.6412 3571.4524
[8,] 2333.1651 6023.6412
[9,] 469.5876 2333.1651
[10,] 127.7753 469.5876
[11,] -1074.7476 127.7753
[12,] -2705.0691 -1074.7476
[13,] -3690.3949 -2705.0691
[14,] -5626.9122 -3690.3949
[15,] -4664.5338 -5626.9122
[16,] -4539.1159 -4664.5338
[17,] -4963.4755 -4539.1159
[18,] -830.3803 -4963.4755
[19,] 686.7275 -830.3803
[20,] 1420.9304 686.7275
[21,] 517.6318 1420.9304
[22,] -724.9953 517.6318
[23,] -1116.6260 -724.9953
[24,] -3255.7698 -1116.6260
[25,] -2276.2030 -3255.7698
[26,] -1074.4015 -2276.2030
[27,] -882.4598 -1074.4015
[28,] -473.7746 -882.4598
[29,] 283.5406 -473.7746
[30,] -306.6294 283.5406
[31,] -1359.9814 -306.6294
[32,] -1253.8690 -1359.9814
[33,] -128.4335 -1253.8690
[34,] 2259.1828 -128.4335
[35,] 2484.6023 2259.1828
[36,] 851.8254 2484.6023
[37,] 2193.2533 851.8254
[38,] 2279.2125 2193.2533
[39,] 4787.6217 2279.2125
[40,] 4380.4055 4787.6217
[41,] 4487.3358 4380.4055
[42,] 3065.8592 4487.3358
[43,] 4099.4087 3065.8592
[44,] 5598.1747 4099.4087
[45,] 4999.5982 5598.1747
[46,] 4864.1628 4999.5982
[47,] 4565.1181 4864.1628
[48,] 2895.0146 4565.1181
[49,] 2995.1456 2895.0146
[50,] 5890.9567 2995.1456
[51,] 3611.2047 5890.9567
[52,] 4398.6819 3611.2047
[53,] 4223.7409 4398.6819
[54,] 87.2218 4223.7409
[55,] -3467.1280 87.2218
[56,] -1797.5789 -3467.1280
[57,] 2996.8178 -1797.5789
[58,] 3003.2539 2996.8178
[59,] 4035.3462 3003.2539
[60,] -773.9683 4035.3462
[61,] -415.8570 -773.9683
[62,] -148.7587 -415.8570
[63,] -2307.9759 -148.7587
[64,] -1731.1914 -2307.9759
[65,] -2657.6583 -1731.1914
[66,] -5587.5237 -2657.6583
[67,] -5982.6679 -5587.5237
[68,] -6300.8224 -5982.6679
[69,] -8855.2018 -6300.8224
[70,] -9529.3795 -8855.2018
[71,] -8893.6929 -9529.3795
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1194.0560 2987.9672
2 -1320.0968 1194.0560
3 -543.8570 -1320.0968
4 -2035.0055 -543.8570
5 -1373.4834 -2035.0055
6 3571.4524 -1373.4834
7 6023.6412 3571.4524
8 2333.1651 6023.6412
9 469.5876 2333.1651
10 127.7753 469.5876
11 -1074.7476 127.7753
12 -2705.0691 -1074.7476
13 -3690.3949 -2705.0691
14 -5626.9122 -3690.3949
15 -4664.5338 -5626.9122
16 -4539.1159 -4664.5338
17 -4963.4755 -4539.1159
18 -830.3803 -4963.4755
19 686.7275 -830.3803
20 1420.9304 686.7275
21 517.6318 1420.9304
22 -724.9953 517.6318
23 -1116.6260 -724.9953
24 -3255.7698 -1116.6260
25 -2276.2030 -3255.7698
26 -1074.4015 -2276.2030
27 -882.4598 -1074.4015
28 -473.7746 -882.4598
29 283.5406 -473.7746
30 -306.6294 283.5406
31 -1359.9814 -306.6294
32 -1253.8690 -1359.9814
33 -128.4335 -1253.8690
34 2259.1828 -128.4335
35 2484.6023 2259.1828
36 851.8254 2484.6023
37 2193.2533 851.8254
38 2279.2125 2193.2533
39 4787.6217 2279.2125
40 4380.4055 4787.6217
41 4487.3358 4380.4055
42 3065.8592 4487.3358
43 4099.4087 3065.8592
44 5598.1747 4099.4087
45 4999.5982 5598.1747
46 4864.1628 4999.5982
47 4565.1181 4864.1628
48 2895.0146 4565.1181
49 2995.1456 2895.0146
50 5890.9567 2995.1456
51 3611.2047 5890.9567
52 4398.6819 3611.2047
53 4223.7409 4398.6819
54 87.2218 4223.7409
55 -3467.1280 87.2218
56 -1797.5789 -3467.1280
57 2996.8178 -1797.5789
58 3003.2539 2996.8178
59 4035.3462 3003.2539
60 -773.9683 4035.3462
61 -415.8570 -773.9683
62 -148.7587 -415.8570
63 -2307.9759 -148.7587
64 -1731.1914 -2307.9759
65 -2657.6583 -1731.1914
66 -5587.5237 -2657.6583
67 -5982.6679 -5587.5237
68 -6300.8224 -5982.6679
69 -8855.2018 -6300.8224
70 -9529.3795 -8855.2018
71 -8893.6929 -9529.3795
> 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/7ilkk1258649687.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/889ma1258649687.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/946z01258649687.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/10jnuk1258649687.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/11nili1258649687.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/12wivu1258649687.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/13gacs1258649687.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/14lvcm1258649687.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/15hhw11258649687.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/16pi0n1258649687.tab")
+ }
>
> system("convert tmp/1buhs1258649687.ps tmp/1buhs1258649687.png")
> system("convert tmp/2vahx1258649687.ps tmp/2vahx1258649687.png")
> system("convert tmp/320ff1258649687.ps tmp/320ff1258649687.png")
> system("convert tmp/472sw1258649687.ps tmp/472sw1258649687.png")
> system("convert tmp/5ilvq1258649687.ps tmp/5ilvq1258649687.png")
> system("convert tmp/6h5ld1258649687.ps tmp/6h5ld1258649687.png")
> system("convert tmp/7ilkk1258649687.ps tmp/7ilkk1258649687.png")
> system("convert tmp/889ma1258649687.ps tmp/889ma1258649687.png")
> system("convert tmp/946z01258649687.ps tmp/946z01258649687.png")
> system("convert tmp/10jnuk1258649687.ps tmp/10jnuk1258649687.png")
>
>
> proc.time()
user system elapsed
2.577 1.585 2.976