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(7024
+ ,2735
+ ,6981
+ ,6962
+ ,6699
+ ,6539
+ ,6940
+ ,2659
+ ,7024
+ ,6981
+ ,6962
+ ,6699
+ ,6774
+ ,2654
+ ,6940
+ ,7024
+ ,6981
+ ,6962
+ ,6671
+ ,2670
+ ,6774
+ ,6940
+ ,7024
+ ,6981
+ ,6965
+ ,2785
+ ,6671
+ ,6774
+ ,6940
+ ,7024
+ ,6969
+ ,2845
+ ,6965
+ ,6671
+ ,6774
+ ,6940
+ ,6822
+ ,2723
+ ,6969
+ ,6965
+ ,6671
+ ,6774
+ ,6878
+ ,2746
+ ,6822
+ ,6969
+ ,6965
+ ,6671
+ ,6691
+ ,2767
+ ,6878
+ ,6822
+ ,6969
+ ,6965
+ ,6837
+ ,2940
+ ,6691
+ ,6878
+ ,6822
+ ,6969
+ ,7018
+ ,2977
+ ,6837
+ ,6691
+ ,6878
+ ,6822
+ ,7167
+ ,2993
+ ,7018
+ ,6837
+ ,6691
+ ,6878
+ ,7076
+ ,2892
+ ,7167
+ ,7018
+ ,6837
+ ,6691
+ ,7171
+ ,2824
+ ,7076
+ ,7167
+ ,7018
+ ,6837
+ ,7093
+ ,2771
+ ,7171
+ ,7076
+ ,7167
+ ,7018
+ ,6971
+ ,2686
+ ,7093
+ ,7171
+ ,7076
+ ,7167
+ ,7142
+ ,2738
+ ,6971
+ ,7093
+ ,7171
+ ,7076
+ ,7047
+ ,2723
+ ,7142
+ ,6971
+ ,7093
+ ,7171
+ ,6999
+ ,2731
+ ,7047
+ ,7142
+ ,6971
+ ,7093
+ ,6650
+ ,2632
+ ,6999
+ ,7047
+ ,7142
+ ,6971
+ ,6475
+ ,2606
+ ,6650
+ ,6999
+ ,7047
+ ,7142
+ ,6437
+ ,2605
+ ,6475
+ ,6650
+ ,6999
+ ,7047
+ ,6639
+ ,2646
+ ,6437
+ ,6475
+ ,6650
+ ,6999
+ ,6422
+ ,2627
+ ,6639
+ ,6437
+ ,6475
+ ,6650
+ ,6272
+ ,2535
+ ,6422
+ ,6639
+ ,6437
+ ,6475
+ ,6232
+ ,2456
+ ,6272
+ ,6422
+ ,6639
+ ,6437
+ ,6003
+ ,2404
+ ,6232
+ ,6272
+ ,6422
+ ,6639
+ ,5673
+ ,2319
+ ,6003
+ ,6232
+ ,6272
+ ,6422
+ ,6050
+ ,2519
+ ,5673
+ ,6003
+ ,6232
+ ,6272
+ ,5977
+ ,2504
+ ,6050
+ ,5673
+ ,6003
+ ,6232
+ ,5796
+ ,2382
+ ,5977
+ ,6050
+ ,5673
+ ,6003
+ ,5752
+ ,2394
+ ,5796
+ ,5977
+ ,6050
+ ,5673
+ ,5609
+ ,2381
+ ,5752
+ ,5796
+ ,5977
+ ,6050
+ ,5839
+ ,2501
+ ,5609
+ ,5752
+ ,5796
+ ,5977
+ ,6069
+ ,2532
+ ,5839
+ ,5609
+ ,5752
+ ,5796
+ ,6006
+ ,2515
+ ,6069
+ ,5839
+ ,5609
+ ,5752
+ ,5809
+ ,2429
+ ,6006
+ ,6069
+ ,5839
+ ,5609
+ ,5797
+ ,2389
+ ,5809
+ ,6006
+ ,6069
+ ,5839
+ ,5502
+ ,2261
+ ,5797
+ ,5809
+ ,6006
+ ,6069
+ ,5568
+ ,2272
+ ,5502
+ ,5797
+ ,5809
+ ,6006
+ ,5864
+ ,2439
+ ,5568
+ ,5502
+ ,5797
+ ,5809
+ ,5764
+ ,2373
+ ,5864
+ ,5568
+ ,5502
+ ,5797
+ ,5615
+ ,2327
+ ,5764
+ ,5864
+ ,5568
+ ,5502
+ ,5615
+ ,2364
+ ,5615
+ ,5764
+ ,5864
+ ,5568
+ ,5681
+ ,2388
+ ,5615
+ ,5615
+ ,5764
+ ,5864
+ ,5915
+ ,2553
+ ,5681
+ ,5615
+ ,5615
+ ,5764
+ ,6334
+ ,2663
+ ,5915
+ ,5681
+ ,5615
+ ,5615
+ ,6494
+ ,2694
+ ,6334
+ ,5915
+ ,5681
+ ,5615
+ ,6620
+ ,2679
+ ,6494
+ ,6334
+ ,5915
+ ,5681
+ ,6578
+ ,2611
+ ,6620
+ ,6494
+ ,6334
+ ,5915
+ ,6495
+ ,2580
+ ,6578
+ ,6620
+ ,6494
+ ,6334
+ ,6538
+ ,2627
+ ,6495
+ ,6578
+ ,6620
+ ,6494
+ ,6737
+ ,2732
+ ,6538
+ ,6495
+ ,6578
+ ,6620
+ ,6651
+ ,2707
+ ,6737
+ ,6538
+ ,6495
+ ,6578
+ ,6530
+ ,2633
+ ,6651
+ ,6737
+ ,6538
+ ,6495
+ ,6563
+ ,2683
+ ,6530
+ ,6651
+ ,6737
+ ,6538)
+ ,dim=c(6
+ ,56)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y-1'
+ ,'Y-2'
+ ,'Y-3'
+ ,'Y-4')
+ ,1:56))
> y <- array(NA,dim=c(6,56),dimnames=list(c('Y','X','Y-1','Y-2','Y-3','Y-4'),1:56))
> 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'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> 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 Y-1 Y-2 Y-3 Y-4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 7024 2735 6981 6962 6699 6539 1 0 0 0 0 0 0 0 0 0 0 1
2 6940 2659 7024 6981 6962 6699 0 1 0 0 0 0 0 0 0 0 0 2
3 6774 2654 6940 7024 6981 6962 0 0 1 0 0 0 0 0 0 0 0 3
4 6671 2670 6774 6940 7024 6981 0 0 0 1 0 0 0 0 0 0 0 4
5 6965 2785 6671 6774 6940 7024 0 0 0 0 1 0 0 0 0 0 0 5
6 6969 2845 6965 6671 6774 6940 0 0 0 0 0 1 0 0 0 0 0 6
7 6822 2723 6969 6965 6671 6774 0 0 0 0 0 0 1 0 0 0 0 7
8 6878 2746 6822 6969 6965 6671 0 0 0 0 0 0 0 1 0 0 0 8
9 6691 2767 6878 6822 6969 6965 0 0 0 0 0 0 0 0 1 0 0 9
10 6837 2940 6691 6878 6822 6969 0 0 0 0 0 0 0 0 0 1 0 10
11 7018 2977 6837 6691 6878 6822 0 0 0 0 0 0 0 0 0 0 1 11
12 7167 2993 7018 6837 6691 6878 0 0 0 0 0 0 0 0 0 0 0 12
13 7076 2892 7167 7018 6837 6691 1 0 0 0 0 0 0 0 0 0 0 13
14 7171 2824 7076 7167 7018 6837 0 1 0 0 0 0 0 0 0 0 0 14
15 7093 2771 7171 7076 7167 7018 0 0 1 0 0 0 0 0 0 0 0 15
16 6971 2686 7093 7171 7076 7167 0 0 0 1 0 0 0 0 0 0 0 16
17 7142 2738 6971 7093 7171 7076 0 0 0 0 1 0 0 0 0 0 0 17
18 7047 2723 7142 6971 7093 7171 0 0 0 0 0 1 0 0 0 0 0 18
19 6999 2731 7047 7142 6971 7093 0 0 0 0 0 0 1 0 0 0 0 19
20 6650 2632 6999 7047 7142 6971 0 0 0 0 0 0 0 1 0 0 0 20
21 6475 2606 6650 6999 7047 7142 0 0 0 0 0 0 0 0 1 0 0 21
22 6437 2605 6475 6650 6999 7047 0 0 0 0 0 0 0 0 0 1 0 22
23 6639 2646 6437 6475 6650 6999 0 0 0 0 0 0 0 0 0 0 1 23
24 6422 2627 6639 6437 6475 6650 0 0 0 0 0 0 0 0 0 0 0 24
25 6272 2535 6422 6639 6437 6475 1 0 0 0 0 0 0 0 0 0 0 25
26 6232 2456 6272 6422 6639 6437 0 1 0 0 0 0 0 0 0 0 0 26
27 6003 2404 6232 6272 6422 6639 0 0 1 0 0 0 0 0 0 0 0 27
28 5673 2319 6003 6232 6272 6422 0 0 0 1 0 0 0 0 0 0 0 28
29 6050 2519 5673 6003 6232 6272 0 0 0 0 1 0 0 0 0 0 0 29
30 5977 2504 6050 5673 6003 6232 0 0 0 0 0 1 0 0 0 0 0 30
31 5796 2382 5977 6050 5673 6003 0 0 0 0 0 0 1 0 0 0 0 31
32 5752 2394 5796 5977 6050 5673 0 0 0 0 0 0 0 1 0 0 0 32
33 5609 2381 5752 5796 5977 6050 0 0 0 0 0 0 0 0 1 0 0 33
34 5839 2501 5609 5752 5796 5977 0 0 0 0 0 0 0 0 0 1 0 34
35 6069 2532 5839 5609 5752 5796 0 0 0 0 0 0 0 0 0 0 1 35
36 6006 2515 6069 5839 5609 5752 0 0 0 0 0 0 0 0 0 0 0 36
37 5809 2429 6006 6069 5839 5609 1 0 0 0 0 0 0 0 0 0 0 37
38 5797 2389 5809 6006 6069 5839 0 1 0 0 0 0 0 0 0 0 0 38
39 5502 2261 5797 5809 6006 6069 0 0 1 0 0 0 0 0 0 0 0 39
40 5568 2272 5502 5797 5809 6006 0 0 0 1 0 0 0 0 0 0 0 40
41 5864 2439 5568 5502 5797 5809 0 0 0 0 1 0 0 0 0 0 0 41
42 5764 2373 5864 5568 5502 5797 0 0 0 0 0 1 0 0 0 0 0 42
43 5615 2327 5764 5864 5568 5502 0 0 0 0 0 0 1 0 0 0 0 43
44 5615 2364 5615 5764 5864 5568 0 0 0 0 0 0 0 1 0 0 0 44
45 5681 2388 5615 5615 5764 5864 0 0 0 0 0 0 0 0 1 0 0 45
46 5915 2553 5681 5615 5615 5764 0 0 0 0 0 0 0 0 0 1 0 46
47 6334 2663 5915 5681 5615 5615 0 0 0 0 0 0 0 0 0 0 1 47
48 6494 2694 6334 5915 5681 5615 0 0 0 0 0 0 0 0 0 0 0 48
49 6620 2679 6494 6334 5915 5681 1 0 0 0 0 0 0 0 0 0 0 49
50 6578 2611 6620 6494 6334 5915 0 1 0 0 0 0 0 0 0 0 0 50
51 6495 2580 6578 6620 6494 6334 0 0 1 0 0 0 0 0 0 0 0 51
52 6538 2627 6495 6578 6620 6494 0 0 0 1 0 0 0 0 0 0 0 52
53 6737 2732 6538 6495 6578 6620 0 0 0 0 1 0 0 0 0 0 0 53
54 6651 2707 6737 6538 6495 6578 0 0 0 0 0 1 0 0 0 0 0 54
55 6530 2633 6651 6737 6538 6495 0 0 0 0 0 0 1 0 0 0 0 55
56 6563 2683 6530 6651 6737 6538 0 0 0 0 0 0 0 1 0 0 0 56
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X `Y-1` `Y-2` `Y-3` `Y-4`
-620.64552 1.12560 0.58825 0.20224 -0.11945 -0.04458
M1 M2 M3 M4 M5 M6
-3.14892 121.75108 45.26329 74.96962 280.84583 64.06565
M7 M8 M9 M10 M11 t
-16.32052 35.63891 -1.96193 75.23377 193.58405 0.18437
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-157.220 -42.773 4.715 43.507 146.237
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -620.64552 275.31328 -2.254 0.030022 *
X 1.12560 0.18896 5.957 6.52e-07 ***
`Y-1` 0.58825 0.13735 4.283 0.000121 ***
`Y-2` 0.20224 0.16725 1.209 0.234054
`Y-3` -0.11945 0.16684 -0.716 0.478410
`Y-4` -0.04458 0.12344 -0.361 0.719963
M1 -3.14892 77.24258 -0.041 0.967695
M2 121.75108 89.33078 1.363 0.180929
M3 45.26329 75.02506 0.603 0.549888
M4 74.96962 78.71369 0.952 0.346896
M5 280.84583 78.65633 3.571 0.000986 ***
M6 64.06565 62.49997 1.025 0.311820
M7 -16.32052 73.53264 -0.222 0.825542
M8 35.63891 94.44294 0.377 0.708006
M9 -1.96193 76.09124 -0.026 0.979565
M10 75.23377 78.19029 0.962 0.342038
M11 193.58405 69.29045 2.794 0.008116 **
t 0.18437 0.99801 0.185 0.854421
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 82.39 on 38 degrees of freedom
Multiple R-squared: 0.9825, Adjusted R-squared: 0.9747
F-statistic: 125.6 on 17 and 38 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.9366068 0.12678648 0.06339324
[2,] 0.9503221 0.09935588 0.04967794
[3,] 0.9690222 0.06195569 0.03097785
[4,] 0.9729479 0.05410416 0.02705208
[5,] 0.9591191 0.08176170 0.04088085
[6,] 0.9898534 0.02029314 0.01014657
[7,] 0.9844110 0.03117794 0.01558897
[8,] 0.9766253 0.04674948 0.02337474
[9,] 0.9642354 0.07152912 0.03576456
[10,] 0.9410269 0.11794615 0.05897307
[11,] 0.9276224 0.14475511 0.07237755
[12,] 0.8739842 0.25203164 0.12601582
[13,] 0.9494581 0.10108374 0.05054187
[14,] 0.9079505 0.18409901 0.09204950
[15,] 0.9182036 0.16359281 0.08179641
> postscript(file="/var/www/html/rcomp/tmp/1ex0k1258726945.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/29ifh1258726945.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/3axl71258726945.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/4n1zv1258726945.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/51wri1258726945.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 = 56
Frequency = 1
1 2 3 4 5 6
146.237206 32.108480 2.751774 -27.527471 17.013756 -5.614146
7 8 9 10 11 12
-16.605076 77.551044 -85.297138 -130.108063 -157.220429 11.329793
13 14 15 16 17 18
-78.173454 19.808709 66.155920 32.385171 33.625013 91.105989
19 20 21 22 23 24
117.553850 -109.721253 -6.755740 42.546447 93.781344 -56.037040
25 26 27 28 29 30
-25.060984 53.334142 -3.879045 -152.886777 21.903556 -1.784316
31 32 33 34 35 36
-48.188818 -6.284848 -26.658357 59.034132 15.903332 -35.417047
37 38 39 40 41 42
-121.009765 -46.716959 -71.707619 101.641311 14.225981 -18.130260
43 44 45 46 47 48
-41.458030 10.922862 118.711236 28.527484 47.535754 80.124294
49 50 51 52 53 54
78.006997 -58.534371 6.678970 46.387767 -86.768305 -65.577266
55 56
-11.301925 27.532195
> postscript(file="/var/www/html/rcomp/tmp/6ihx21258726945.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 = 56
Frequency = 1
lag(myerror, k = 1) myerror
0 146.237206 NA
1 32.108480 146.237206
2 2.751774 32.108480
3 -27.527471 2.751774
4 17.013756 -27.527471
5 -5.614146 17.013756
6 -16.605076 -5.614146
7 77.551044 -16.605076
8 -85.297138 77.551044
9 -130.108063 -85.297138
10 -157.220429 -130.108063
11 11.329793 -157.220429
12 -78.173454 11.329793
13 19.808709 -78.173454
14 66.155920 19.808709
15 32.385171 66.155920
16 33.625013 32.385171
17 91.105989 33.625013
18 117.553850 91.105989
19 -109.721253 117.553850
20 -6.755740 -109.721253
21 42.546447 -6.755740
22 93.781344 42.546447
23 -56.037040 93.781344
24 -25.060984 -56.037040
25 53.334142 -25.060984
26 -3.879045 53.334142
27 -152.886777 -3.879045
28 21.903556 -152.886777
29 -1.784316 21.903556
30 -48.188818 -1.784316
31 -6.284848 -48.188818
32 -26.658357 -6.284848
33 59.034132 -26.658357
34 15.903332 59.034132
35 -35.417047 15.903332
36 -121.009765 -35.417047
37 -46.716959 -121.009765
38 -71.707619 -46.716959
39 101.641311 -71.707619
40 14.225981 101.641311
41 -18.130260 14.225981
42 -41.458030 -18.130260
43 10.922862 -41.458030
44 118.711236 10.922862
45 28.527484 118.711236
46 47.535754 28.527484
47 80.124294 47.535754
48 78.006997 80.124294
49 -58.534371 78.006997
50 6.678970 -58.534371
51 46.387767 6.678970
52 -86.768305 46.387767
53 -65.577266 -86.768305
54 -11.301925 -65.577266
55 27.532195 -11.301925
56 NA 27.532195
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 32.108480 146.237206
[2,] 2.751774 32.108480
[3,] -27.527471 2.751774
[4,] 17.013756 -27.527471
[5,] -5.614146 17.013756
[6,] -16.605076 -5.614146
[7,] 77.551044 -16.605076
[8,] -85.297138 77.551044
[9,] -130.108063 -85.297138
[10,] -157.220429 -130.108063
[11,] 11.329793 -157.220429
[12,] -78.173454 11.329793
[13,] 19.808709 -78.173454
[14,] 66.155920 19.808709
[15,] 32.385171 66.155920
[16,] 33.625013 32.385171
[17,] 91.105989 33.625013
[18,] 117.553850 91.105989
[19,] -109.721253 117.553850
[20,] -6.755740 -109.721253
[21,] 42.546447 -6.755740
[22,] 93.781344 42.546447
[23,] -56.037040 93.781344
[24,] -25.060984 -56.037040
[25,] 53.334142 -25.060984
[26,] -3.879045 53.334142
[27,] -152.886777 -3.879045
[28,] 21.903556 -152.886777
[29,] -1.784316 21.903556
[30,] -48.188818 -1.784316
[31,] -6.284848 -48.188818
[32,] -26.658357 -6.284848
[33,] 59.034132 -26.658357
[34,] 15.903332 59.034132
[35,] -35.417047 15.903332
[36,] -121.009765 -35.417047
[37,] -46.716959 -121.009765
[38,] -71.707619 -46.716959
[39,] 101.641311 -71.707619
[40,] 14.225981 101.641311
[41,] -18.130260 14.225981
[42,] -41.458030 -18.130260
[43,] 10.922862 -41.458030
[44,] 118.711236 10.922862
[45,] 28.527484 118.711236
[46,] 47.535754 28.527484
[47,] 80.124294 47.535754
[48,] 78.006997 80.124294
[49,] -58.534371 78.006997
[50,] 6.678970 -58.534371
[51,] 46.387767 6.678970
[52,] -86.768305 46.387767
[53,] -65.577266 -86.768305
[54,] -11.301925 -65.577266
[55,] 27.532195 -11.301925
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 32.108480 146.237206
2 2.751774 32.108480
3 -27.527471 2.751774
4 17.013756 -27.527471
5 -5.614146 17.013756
6 -16.605076 -5.614146
7 77.551044 -16.605076
8 -85.297138 77.551044
9 -130.108063 -85.297138
10 -157.220429 -130.108063
11 11.329793 -157.220429
12 -78.173454 11.329793
13 19.808709 -78.173454
14 66.155920 19.808709
15 32.385171 66.155920
16 33.625013 32.385171
17 91.105989 33.625013
18 117.553850 91.105989
19 -109.721253 117.553850
20 -6.755740 -109.721253
21 42.546447 -6.755740
22 93.781344 42.546447
23 -56.037040 93.781344
24 -25.060984 -56.037040
25 53.334142 -25.060984
26 -3.879045 53.334142
27 -152.886777 -3.879045
28 21.903556 -152.886777
29 -1.784316 21.903556
30 -48.188818 -1.784316
31 -6.284848 -48.188818
32 -26.658357 -6.284848
33 59.034132 -26.658357
34 15.903332 59.034132
35 -35.417047 15.903332
36 -121.009765 -35.417047
37 -46.716959 -121.009765
38 -71.707619 -46.716959
39 101.641311 -71.707619
40 14.225981 101.641311
41 -18.130260 14.225981
42 -41.458030 -18.130260
43 10.922862 -41.458030
44 118.711236 10.922862
45 28.527484 118.711236
46 47.535754 28.527484
47 80.124294 47.535754
48 78.006997 80.124294
49 -58.534371 78.006997
50 6.678970 -58.534371
51 46.387767 6.678970
52 -86.768305 46.387767
53 -65.577266 -86.768305
54 -11.301925 -65.577266
55 27.532195 -11.301925
> 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/70jr51258726945.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/88hac1258726945.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/97fdv1258726945.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/10w5f71258726945.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/11g0xd1258726945.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/12gw7o1258726945.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/13qjpm1258726945.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/14z0wa1258726945.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/157tcy1258726945.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/16wgzm1258726945.tab")
+ }
>
> system("convert tmp/1ex0k1258726945.ps tmp/1ex0k1258726945.png")
> system("convert tmp/29ifh1258726945.ps tmp/29ifh1258726945.png")
> system("convert tmp/3axl71258726945.ps tmp/3axl71258726945.png")
> system("convert tmp/4n1zv1258726945.ps tmp/4n1zv1258726945.png")
> system("convert tmp/51wri1258726945.ps tmp/51wri1258726945.png")
> system("convert tmp/6ihx21258726945.ps tmp/6ihx21258726945.png")
> system("convert tmp/70jr51258726945.ps tmp/70jr51258726945.png")
> system("convert tmp/88hac1258726945.ps tmp/88hac1258726945.png")
> system("convert tmp/97fdv1258726945.ps tmp/97fdv1258726945.png")
> system("convert tmp/10w5f71258726945.ps tmp/10w5f71258726945.png")
>
>
> proc.time()
user system elapsed
2.340 1.561 2.705