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(627,0,696,0,825,0,677,0,656,0,785,0,412,0,352,0,839,0,729,0,696,0,641,0,695,0,638,0,762,0,635,0,721,0,854,0,418,0,367,0,824,0,687,0,601,0,676,0,740,0,691,0,683,0,594,0,729,0,731,0,386,0,331,0,707,0,715,0,657,0,653,0,642,0,643,0,718,0,654,0,632,0,731,0,392,1,344,1,792,1,852,1,649,1,629,1,685,1,617,1,715,1,715,1,629,1,916,1,531,1,357,1,917,1,828,1,708,1,858,1,775,1,785,1,1006,1,789,1,734,1,906,1,532,1,387,1,991,1,841,1),dim=c(2,70),dimnames=list(c('Y','X'),1:70))
> y <- array(NA,dim=c(2,70),dimnames=list(c('Y','X'),1:70))
> 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'
> #'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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 627 0 1 0 0 0 0 0 0 0 0 0 0
2 696 0 0 1 0 0 0 0 0 0 0 0 0
3 825 0 0 0 1 0 0 0 0 0 0 0 0
4 677 0 0 0 0 1 0 0 0 0 0 0 0
5 656 0 0 0 0 0 1 0 0 0 0 0 0
6 785 0 0 0 0 0 0 1 0 0 0 0 0
7 412 0 0 0 0 0 0 0 1 0 0 0 0
8 352 0 0 0 0 0 0 0 0 1 0 0 0
9 839 0 0 0 0 0 0 0 0 0 1 0 0
10 729 0 0 0 0 0 0 0 0 0 0 1 0
11 696 0 0 0 0 0 0 0 0 0 0 0 1
12 641 0 0 0 0 0 0 0 0 0 0 0 0
13 695 0 1 0 0 0 0 0 0 0 0 0 0
14 638 0 0 1 0 0 0 0 0 0 0 0 0
15 762 0 0 0 1 0 0 0 0 0 0 0 0
16 635 0 0 0 0 1 0 0 0 0 0 0 0
17 721 0 0 0 0 0 1 0 0 0 0 0 0
18 854 0 0 0 0 0 0 1 0 0 0 0 0
19 418 0 0 0 0 0 0 0 1 0 0 0 0
20 367 0 0 0 0 0 0 0 0 1 0 0 0
21 824 0 0 0 0 0 0 0 0 0 1 0 0
22 687 0 0 0 0 0 0 0 0 0 0 1 0
23 601 0 0 0 0 0 0 0 0 0 0 0 1
24 676 0 0 0 0 0 0 0 0 0 0 0 0
25 740 0 1 0 0 0 0 0 0 0 0 0 0
26 691 0 0 1 0 0 0 0 0 0 0 0 0
27 683 0 0 0 1 0 0 0 0 0 0 0 0
28 594 0 0 0 0 1 0 0 0 0 0 0 0
29 729 0 0 0 0 0 1 0 0 0 0 0 0
30 731 0 0 0 0 0 0 1 0 0 0 0 0
31 386 0 0 0 0 0 0 0 1 0 0 0 0
32 331 0 0 0 0 0 0 0 0 1 0 0 0
33 707 0 0 0 0 0 0 0 0 0 1 0 0
34 715 0 0 0 0 0 0 0 0 0 0 1 0
35 657 0 0 0 0 0 0 0 0 0 0 0 1
36 653 0 0 0 0 0 0 0 0 0 0 0 0
37 642 0 1 0 0 0 0 0 0 0 0 0 0
38 643 0 0 1 0 0 0 0 0 0 0 0 0
39 718 0 0 0 1 0 0 0 0 0 0 0 0
40 654 0 0 0 0 1 0 0 0 0 0 0 0
41 632 0 0 0 0 0 1 0 0 0 0 0 0
42 731 0 0 0 0 0 0 1 0 0 0 0 0
43 392 1 0 0 0 0 0 0 1 0 0 0 0
44 344 1 0 0 0 0 0 0 0 1 0 0 0
45 792 1 0 0 0 0 0 0 0 0 1 0 0
46 852 1 0 0 0 0 0 0 0 0 0 1 0
47 649 1 0 0 0 0 0 0 0 0 0 0 1
48 629 1 0 0 0 0 0 0 0 0 0 0 0
49 685 1 1 0 0 0 0 0 0 0 0 0 0
50 617 1 0 1 0 0 0 0 0 0 0 0 0
51 715 1 0 0 1 0 0 0 0 0 0 0 0
52 715 1 0 0 0 1 0 0 0 0 0 0 0
53 629 1 0 0 0 0 1 0 0 0 0 0 0
54 916 1 0 0 0 0 0 1 0 0 0 0 0
55 531 1 0 0 0 0 0 0 1 0 0 0 0
56 357 1 0 0 0 0 0 0 0 1 0 0 0
57 917 1 0 0 0 0 0 0 0 0 1 0 0
58 828 1 0 0 0 0 0 0 0 0 0 1 0
59 708 1 0 0 0 0 0 0 0 0 0 0 1
60 858 1 0 0 0 0 0 0 0 0 0 0 0
61 775 1 1 0 0 0 0 0 0 0 0 0 0
62 785 1 0 1 0 0 0 0 0 0 0 0 0
63 1006 1 0 0 1 0 0 0 0 0 0 0 0
64 789 1 0 0 0 1 0 0 0 0 0 0 0
65 734 1 0 0 0 0 1 0 0 0 0 0 0
66 906 1 0 0 0 0 0 1 0 0 0 0 0
67 532 1 0 0 0 0 0 0 1 0 0 0 0
68 387 1 0 0 0 0 0 0 0 1 0 0 0
69 991 1 0 0 0 0 0 0 0 0 1 0 0
70 841 1 0 0 0 0 0 0 0 0 0 1 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
661.393 75.018 7.601 -8.065 98.435 -9.065
M5 M6 M7 M8 M9 M10
-2.899 134.101 -253.735 -342.568 146.098 76.432
M11
-29.200
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-119.846 -35.430 1.231 37.127 171.154
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 661.393 28.885 22.897 < 2e-16 ***
X 75.018 15.574 4.817 1.12e-05 ***
M1 7.601 38.204 0.199 0.842999
M2 -8.065 38.204 -0.211 0.833551
M3 98.435 38.204 2.577 0.012594 *
M4 -9.065 38.204 -0.237 0.813283
M5 -2.899 38.204 -0.076 0.939783
M6 134.101 38.204 3.510 0.000883 ***
M7 -253.735 38.222 -6.638 1.27e-08 ***
M8 -342.568 38.222 -8.963 1.78e-12 ***
M9 146.098 38.222 3.822 0.000329 ***
M10 76.432 38.222 2.000 0.050308 .
M11 -29.200 39.888 -0.732 0.467140
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 63.07 on 57 degrees of freedom
Multiple R-squared: 0.8653, Adjusted R-squared: 0.8369
F-statistic: 30.51 on 12 and 57 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.2839949488 0.5679898976 0.7160051
[2,] 0.2252729706 0.4505459412 0.7747270
[3,] 0.1895958306 0.3791916611 0.8104042
[4,] 0.1048017631 0.2096035262 0.8951982
[5,] 0.0587301109 0.1174602217 0.9412699
[6,] 0.0299888352 0.0599776703 0.9700112
[7,] 0.0184325136 0.0368650272 0.9815675
[8,] 0.0261332836 0.0522665673 0.9738667
[9,] 0.0149217915 0.0298435829 0.9850782
[10,] 0.0208588885 0.0417177769 0.9791411
[11,] 0.0127749202 0.0255498403 0.9872251
[12,] 0.0286287560 0.0572575121 0.9713712
[13,] 0.0246289238 0.0492578475 0.9753711
[14,] 0.0229282514 0.0458565028 0.9770717
[15,] 0.0273216717 0.0546433434 0.9726783
[16,] 0.0167567238 0.0335134476 0.9832433
[17,] 0.0113592599 0.0227185198 0.9886407
[18,] 0.0267879254 0.0535758508 0.9732121
[19,] 0.0161453401 0.0322906801 0.9838547
[20,] 0.0105986108 0.0211972215 0.9894014
[21,] 0.0058082256 0.0116164512 0.9941918
[22,] 0.0037296764 0.0074593529 0.9962703
[23,] 0.0023039576 0.0046079151 0.9976960
[24,] 0.0013607458 0.0027214916 0.9986393
[25,] 0.0006757843 0.0013515685 0.9993242
[26,] 0.0007281279 0.0014562558 0.9992719
[27,] 0.0004750305 0.0009500609 0.9995250
[28,] 0.0004995146 0.0009990292 0.9995005
[29,] 0.0002292704 0.0004585408 0.9997707
[30,] 0.0003887308 0.0007774616 0.9996113
[31,] 0.0007696859 0.0015393718 0.9992303
[32,] 0.0004182597 0.0008365193 0.9995817
[33,] 0.0023295295 0.0046590591 0.9976705
[34,] 0.0015943559 0.0031887118 0.9984056
[35,] 0.0043913111 0.0087826222 0.9956087
[36,] 0.4888262925 0.9776525849 0.5111737
[37,] 0.5372287024 0.9255425953 0.4627713
[38,] 0.8177980652 0.3644038697 0.1822019
[39,] 0.7206681732 0.5586636535 0.2793318
> postscript(file="/var/www/html/rcomp/tmp/1uyvl1260387826.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/2wiip1260387826.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/356kf1260387826.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/4razq1260387826.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/5su391260387826.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 = 70
Frequency = 1
1 2 3 4 5 6
-41.9939024 42.6727642 65.1727642 24.6727642 -2.4939024 -10.4939024
7 8 9 10 11 12
4.3424797 33.1758130 31.5091463 -8.8241870 63.8073171 -20.3926829
13 14 15 16 17 18
26.0060976 -15.3272358 2.1727642 -17.3272358 62.5060976 58.5060976
19 20 21 22 23 24
10.3424797 48.1758130 16.5091463 -50.8241870 -31.1926829 14.6073171
25 26 27 28 29 30
71.0060976 37.6727642 -76.8272358 -58.3272358 70.5060976 -64.4939024
31 32 33 34 35 36
-21.6575203 12.1758130 -100.4908537 -22.8241870 24.8073171 -8.3926829
37 38 39 40 41 42
-26.9939024 -10.3272358 -41.8272358 1.6727642 -26.4939024 -64.4939024
43 44 45 46 47 48
-90.6758130 -49.8424797 -90.5091463 39.1575203 -58.2109756 -107.4109756
49 50 51 52 53 54
-59.0121951 -111.3455285 -119.8455285 -12.3455285 -104.5121951 45.4878049
55 56 57 58 59 60
48.3241870 -36.8424797 34.4908537 15.1575203 0.7890244 121.5890244
61 62 63 64 65 66
30.9878049 56.6544715 171.1544715 61.6544715 0.4878049 35.4878049
67 68 69 70
49.3241870 -6.8424797 108.4908537 28.1575203
> postscript(file="/var/www/html/rcomp/tmp/6w82a1260387826.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 = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 -41.9939024 NA
1 42.6727642 -41.9939024
2 65.1727642 42.6727642
3 24.6727642 65.1727642
4 -2.4939024 24.6727642
5 -10.4939024 -2.4939024
6 4.3424797 -10.4939024
7 33.1758130 4.3424797
8 31.5091463 33.1758130
9 -8.8241870 31.5091463
10 63.8073171 -8.8241870
11 -20.3926829 63.8073171
12 26.0060976 -20.3926829
13 -15.3272358 26.0060976
14 2.1727642 -15.3272358
15 -17.3272358 2.1727642
16 62.5060976 -17.3272358
17 58.5060976 62.5060976
18 10.3424797 58.5060976
19 48.1758130 10.3424797
20 16.5091463 48.1758130
21 -50.8241870 16.5091463
22 -31.1926829 -50.8241870
23 14.6073171 -31.1926829
24 71.0060976 14.6073171
25 37.6727642 71.0060976
26 -76.8272358 37.6727642
27 -58.3272358 -76.8272358
28 70.5060976 -58.3272358
29 -64.4939024 70.5060976
30 -21.6575203 -64.4939024
31 12.1758130 -21.6575203
32 -100.4908537 12.1758130
33 -22.8241870 -100.4908537
34 24.8073171 -22.8241870
35 -8.3926829 24.8073171
36 -26.9939024 -8.3926829
37 -10.3272358 -26.9939024
38 -41.8272358 -10.3272358
39 1.6727642 -41.8272358
40 -26.4939024 1.6727642
41 -64.4939024 -26.4939024
42 -90.6758130 -64.4939024
43 -49.8424797 -90.6758130
44 -90.5091463 -49.8424797
45 39.1575203 -90.5091463
46 -58.2109756 39.1575203
47 -107.4109756 -58.2109756
48 -59.0121951 -107.4109756
49 -111.3455285 -59.0121951
50 -119.8455285 -111.3455285
51 -12.3455285 -119.8455285
52 -104.5121951 -12.3455285
53 45.4878049 -104.5121951
54 48.3241870 45.4878049
55 -36.8424797 48.3241870
56 34.4908537 -36.8424797
57 15.1575203 34.4908537
58 0.7890244 15.1575203
59 121.5890244 0.7890244
60 30.9878049 121.5890244
61 56.6544715 30.9878049
62 171.1544715 56.6544715
63 61.6544715 171.1544715
64 0.4878049 61.6544715
65 35.4878049 0.4878049
66 49.3241870 35.4878049
67 -6.8424797 49.3241870
68 108.4908537 -6.8424797
69 28.1575203 108.4908537
70 NA 28.1575203
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 42.6727642 -41.9939024
[2,] 65.1727642 42.6727642
[3,] 24.6727642 65.1727642
[4,] -2.4939024 24.6727642
[5,] -10.4939024 -2.4939024
[6,] 4.3424797 -10.4939024
[7,] 33.1758130 4.3424797
[8,] 31.5091463 33.1758130
[9,] -8.8241870 31.5091463
[10,] 63.8073171 -8.8241870
[11,] -20.3926829 63.8073171
[12,] 26.0060976 -20.3926829
[13,] -15.3272358 26.0060976
[14,] 2.1727642 -15.3272358
[15,] -17.3272358 2.1727642
[16,] 62.5060976 -17.3272358
[17,] 58.5060976 62.5060976
[18,] 10.3424797 58.5060976
[19,] 48.1758130 10.3424797
[20,] 16.5091463 48.1758130
[21,] -50.8241870 16.5091463
[22,] -31.1926829 -50.8241870
[23,] 14.6073171 -31.1926829
[24,] 71.0060976 14.6073171
[25,] 37.6727642 71.0060976
[26,] -76.8272358 37.6727642
[27,] -58.3272358 -76.8272358
[28,] 70.5060976 -58.3272358
[29,] -64.4939024 70.5060976
[30,] -21.6575203 -64.4939024
[31,] 12.1758130 -21.6575203
[32,] -100.4908537 12.1758130
[33,] -22.8241870 -100.4908537
[34,] 24.8073171 -22.8241870
[35,] -8.3926829 24.8073171
[36,] -26.9939024 -8.3926829
[37,] -10.3272358 -26.9939024
[38,] -41.8272358 -10.3272358
[39,] 1.6727642 -41.8272358
[40,] -26.4939024 1.6727642
[41,] -64.4939024 -26.4939024
[42,] -90.6758130 -64.4939024
[43,] -49.8424797 -90.6758130
[44,] -90.5091463 -49.8424797
[45,] 39.1575203 -90.5091463
[46,] -58.2109756 39.1575203
[47,] -107.4109756 -58.2109756
[48,] -59.0121951 -107.4109756
[49,] -111.3455285 -59.0121951
[50,] -119.8455285 -111.3455285
[51,] -12.3455285 -119.8455285
[52,] -104.5121951 -12.3455285
[53,] 45.4878049 -104.5121951
[54,] 48.3241870 45.4878049
[55,] -36.8424797 48.3241870
[56,] 34.4908537 -36.8424797
[57,] 15.1575203 34.4908537
[58,] 0.7890244 15.1575203
[59,] 121.5890244 0.7890244
[60,] 30.9878049 121.5890244
[61,] 56.6544715 30.9878049
[62,] 171.1544715 56.6544715
[63,] 61.6544715 171.1544715
[64,] 0.4878049 61.6544715
[65,] 35.4878049 0.4878049
[66,] 49.3241870 35.4878049
[67,] -6.8424797 49.3241870
[68,] 108.4908537 -6.8424797
[69,] 28.1575203 108.4908537
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 42.6727642 -41.9939024
2 65.1727642 42.6727642
3 24.6727642 65.1727642
4 -2.4939024 24.6727642
5 -10.4939024 -2.4939024
6 4.3424797 -10.4939024
7 33.1758130 4.3424797
8 31.5091463 33.1758130
9 -8.8241870 31.5091463
10 63.8073171 -8.8241870
11 -20.3926829 63.8073171
12 26.0060976 -20.3926829
13 -15.3272358 26.0060976
14 2.1727642 -15.3272358
15 -17.3272358 2.1727642
16 62.5060976 -17.3272358
17 58.5060976 62.5060976
18 10.3424797 58.5060976
19 48.1758130 10.3424797
20 16.5091463 48.1758130
21 -50.8241870 16.5091463
22 -31.1926829 -50.8241870
23 14.6073171 -31.1926829
24 71.0060976 14.6073171
25 37.6727642 71.0060976
26 -76.8272358 37.6727642
27 -58.3272358 -76.8272358
28 70.5060976 -58.3272358
29 -64.4939024 70.5060976
30 -21.6575203 -64.4939024
31 12.1758130 -21.6575203
32 -100.4908537 12.1758130
33 -22.8241870 -100.4908537
34 24.8073171 -22.8241870
35 -8.3926829 24.8073171
36 -26.9939024 -8.3926829
37 -10.3272358 -26.9939024
38 -41.8272358 -10.3272358
39 1.6727642 -41.8272358
40 -26.4939024 1.6727642
41 -64.4939024 -26.4939024
42 -90.6758130 -64.4939024
43 -49.8424797 -90.6758130
44 -90.5091463 -49.8424797
45 39.1575203 -90.5091463
46 -58.2109756 39.1575203
47 -107.4109756 -58.2109756
48 -59.0121951 -107.4109756
49 -111.3455285 -59.0121951
50 -119.8455285 -111.3455285
51 -12.3455285 -119.8455285
52 -104.5121951 -12.3455285
53 45.4878049 -104.5121951
54 48.3241870 45.4878049
55 -36.8424797 48.3241870
56 34.4908537 -36.8424797
57 15.1575203 34.4908537
58 0.7890244 15.1575203
59 121.5890244 0.7890244
60 30.9878049 121.5890244
61 56.6544715 30.9878049
62 171.1544715 56.6544715
63 61.6544715 171.1544715
64 0.4878049 61.6544715
65 35.4878049 0.4878049
66 49.3241870 35.4878049
67 -6.8424797 49.3241870
68 108.4908537 -6.8424797
69 28.1575203 108.4908537
> 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/7owcn1260387826.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/8u57y1260387826.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/9umjg1260387826.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/10zx571260387826.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/11djtg1260387826.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/12n6bn1260387826.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/13yeeh1260387826.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/1499331260387826.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/15eh021260387826.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/16v67o1260387826.tab")
+ }
>
> system("convert tmp/1uyvl1260387826.ps tmp/1uyvl1260387826.png")
> system("convert tmp/2wiip1260387826.ps tmp/2wiip1260387826.png")
> system("convert tmp/356kf1260387826.ps tmp/356kf1260387826.png")
> system("convert tmp/4razq1260387826.ps tmp/4razq1260387826.png")
> system("convert tmp/5su391260387826.ps tmp/5su391260387826.png")
> system("convert tmp/6w82a1260387826.ps tmp/6w82a1260387826.png")
> system("convert tmp/7owcn1260387826.ps tmp/7owcn1260387826.png")
> system("convert tmp/8u57y1260387826.ps tmp/8u57y1260387826.png")
> system("convert tmp/9umjg1260387826.ps tmp/9umjg1260387826.png")
> system("convert tmp/10zx571260387826.ps tmp/10zx571260387826.png")
>
>
> proc.time()
user system elapsed
2.601 1.615 3.313