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(8.9,6.3,8.2,6.2,7.6,6.1,7.7,6.3,8.1,6.5,8.3,6.6,8.3,6.5,7.9,6.2,7.8,6.2,8,5.9,8.5,6.1,8.6,6.1,8.5,6.1,8,6.1,7.8,6.1,8,6.4,8.2,6.7,8.3,6.9,8.2,7,8.1,7,8,6.8,7.8,6.4,7.8,5.9,7.7,5.5,7.6,5.5,7.6,5.6,7.6,5.8,7.8,5.9,8,6.1,8,6.1,7.9,6,7.7,6,7.4,5.9,6.9,5.5,6.7,5.6,6.5,5.4,6.4,5.2,6.7,5.2,6.8,5.2,6.9,5.5,6.9,5.8,6.7,5.8,6.4,5.5,6.2,5.3,5.9,5.1,6.1,5.2,6.7,5.8,6.8,5.8,6.6,5.5,6.4,5,6.4,4.9,6.7,5.3,7.1,6.1,7.1,6.5,6.9,6.8,6.4,6.6,6,6.4,6,6.4),dim=c(2,58),dimnames=list(c('X','Y'),1:58))
> y <- array(NA,dim=c(2,58),dimnames=list(c('X','Y'),1:58))
> 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
X Y M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 8.9 6.3 1 0 0 0 0 0 0 0 0 0 0
2 8.2 6.2 0 1 0 0 0 0 0 0 0 0 0
3 7.6 6.1 0 0 1 0 0 0 0 0 0 0 0
4 7.7 6.3 0 0 0 1 0 0 0 0 0 0 0
5 8.1 6.5 0 0 0 0 1 0 0 0 0 0 0
6 8.3 6.6 0 0 0 0 0 1 0 0 0 0 0
7 8.3 6.5 0 0 0 0 0 0 1 0 0 0 0
8 7.9 6.2 0 0 0 0 0 0 0 1 0 0 0
9 7.8 6.2 0 0 0 0 0 0 0 0 1 0 0
10 8.0 5.9 0 0 0 0 0 0 0 0 0 1 0
11 8.5 6.1 0 0 0 0 0 0 0 0 0 0 1
12 8.6 6.1 0 0 0 0 0 0 0 0 0 0 0
13 8.5 6.1 1 0 0 0 0 0 0 0 0 0 0
14 8.0 6.1 0 1 0 0 0 0 0 0 0 0 0
15 7.8 6.1 0 0 1 0 0 0 0 0 0 0 0
16 8.0 6.4 0 0 0 1 0 0 0 0 0 0 0
17 8.2 6.7 0 0 0 0 1 0 0 0 0 0 0
18 8.3 6.9 0 0 0 0 0 1 0 0 0 0 0
19 8.2 7.0 0 0 0 0 0 0 1 0 0 0 0
20 8.1 7.0 0 0 0 0 0 0 0 1 0 0 0
21 8.0 6.8 0 0 0 0 0 0 0 0 1 0 0
22 7.8 6.4 0 0 0 0 0 0 0 0 0 1 0
23 7.8 5.9 0 0 0 0 0 0 0 0 0 0 1
24 7.7 5.5 0 0 0 0 0 0 0 0 0 0 0
25 7.6 5.5 1 0 0 0 0 0 0 0 0 0 0
26 7.6 5.6 0 1 0 0 0 0 0 0 0 0 0
27 7.6 5.8 0 0 1 0 0 0 0 0 0 0 0
28 7.8 5.9 0 0 0 1 0 0 0 0 0 0 0
29 8.0 6.1 0 0 0 0 1 0 0 0 0 0 0
30 8.0 6.1 0 0 0 0 0 1 0 0 0 0 0
31 7.9 6.0 0 0 0 0 0 0 1 0 0 0 0
32 7.7 6.0 0 0 0 0 0 0 0 1 0 0 0
33 7.4 5.9 0 0 0 0 0 0 0 0 1 0 0
34 6.9 5.5 0 0 0 0 0 0 0 0 0 1 0
35 6.7 5.6 0 0 0 0 0 0 0 0 0 0 1
36 6.5 5.4 0 0 0 0 0 0 0 0 0 0 0
37 6.4 5.2 1 0 0 0 0 0 0 0 0 0 0
38 6.7 5.2 0 1 0 0 0 0 0 0 0 0 0
39 6.8 5.2 0 0 1 0 0 0 0 0 0 0 0
40 6.9 5.5 0 0 0 1 0 0 0 0 0 0 0
41 6.9 5.8 0 0 0 0 1 0 0 0 0 0 0
42 6.7 5.8 0 0 0 0 0 1 0 0 0 0 0
43 6.4 5.5 0 0 0 0 0 0 1 0 0 0 0
44 6.2 5.3 0 0 0 0 0 0 0 1 0 0 0
45 5.9 5.1 0 0 0 0 0 0 0 0 1 0 0
46 6.1 5.2 0 0 0 0 0 0 0 0 0 1 0
47 6.7 5.8 0 0 0 0 0 0 0 0 0 0 1
48 6.8 5.8 0 0 0 0 0 0 0 0 0 0 0
49 6.6 5.5 1 0 0 0 0 0 0 0 0 0 0
50 6.4 5.0 0 1 0 0 0 0 0 0 0 0 0
51 6.4 4.9 0 0 1 0 0 0 0 0 0 0 0
52 6.7 5.3 0 0 0 1 0 0 0 0 0 0 0
53 7.1 6.1 0 0 0 0 1 0 0 0 0 0 0
54 7.1 6.5 0 0 0 0 0 1 0 0 0 0 0
55 6.9 6.8 0 0 0 0 0 0 1 0 0 0 0
56 6.4 6.6 0 0 0 0 0 0 0 1 0 0 0
57 6.0 6.4 0 0 0 0 0 0 0 0 1 0 0
58 6.0 6.4 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) Y M1 M2 M3 M4
1.05048 1.11395 0.17772 0.06912 -0.07088 -0.18051
M5 M6 M7 M8 M9 M10
-0.34153 -0.47749 -0.59521 -0.71925 -0.80330 -0.64051
M11
-0.14209
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.53925 -0.25543 0.02644 0.37202 1.01772
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.05048 1.09932 0.956 0.344
Y 1.11395 0.18473 6.030 2.82e-07 ***
M1 0.17772 0.42385 0.419 0.677
M2 0.06912 0.42409 0.163 0.871
M3 -0.07088 0.42409 -0.167 0.868
M4 -0.18051 0.42514 -0.425 0.673
M5 -0.34153 0.43542 -0.784 0.437
M6 -0.47749 0.44206 -1.080 0.286
M7 -0.59521 0.44102 -1.350 0.184
M8 -0.71925 0.43458 -1.655 0.105
M9 -0.80330 0.42961 -1.870 0.068 .
M10 -0.64051 0.42514 -1.507 0.139
M11 -0.14209 0.44762 -0.317 0.752
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6318 on 45 degrees of freedom
Multiple R-squared: 0.4927, Adjusted R-squared: 0.3574
F-statistic: 3.642 on 12 and 45 DF, p-value: 0.0007628
> 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.950589e-03 0.0079011788 0.996049411
[2,] 2.452565e-03 0.0049051306 0.997547435
[3,] 1.888930e-03 0.0037778593 0.998111070
[4,] 1.149644e-03 0.0022992872 0.998850356
[5,] 2.657852e-04 0.0005315704 0.999734215
[6,] 6.835388e-05 0.0001367078 0.999931646
[7,] 5.604556e-05 0.0001120911 0.999943954
[8,] 4.184652e-04 0.0008369304 0.999581535
[9,] 1.477604e-03 0.0029552081 0.998522396
[10,] 4.853358e-03 0.0097067158 0.995146642
[11,] 2.885292e-03 0.0057705833 0.997114708
[12,] 1.661119e-03 0.0033222381 0.998338881
[13,] 1.348453e-03 0.0026969058 0.998651547
[14,] 1.148499e-03 0.0022969984 0.998851501
[15,] 1.291079e-03 0.0025821584 0.998708921
[16,] 2.824112e-03 0.0056482241 0.997175888
[17,] 1.588093e-02 0.0317618672 0.984119066
[18,] 3.061961e-01 0.6123922763 0.693803862
[19,] 9.348666e-01 0.1302667894 0.065133395
[20,] 9.788143e-01 0.0423714164 0.021185708
[21,] 9.921400e-01 0.0157199969 0.007859998
[22,] 9.930768e-01 0.0138464276 0.006923214
[23,] 9.891335e-01 0.0217329994 0.010866500
[24,] 9.874735e-01 0.0250529354 0.012526468
[25,] 9.715381e-01 0.0569238781 0.028461939
[26,] 9.382934e-01 0.1234131366 0.061706568
[27,] 9.131852e-01 0.1736296287 0.086814814
> postscript(file="/var/www/html/rcomp/tmp/1y7ix1258661063.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/2l6nh1258661063.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/33c521258661063.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/4hiiq1258661063.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/58kkl1258661063.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 = 58
Frequency = 1
1 2 3 4 5
0.6539083604 0.1739083604 -0.1746965293 -0.1878594632 0.1503727133
6 7 8 9 10
0.3749307574 0.6040468456 0.6622790221 0.6463258677 1.0177209779
11 12 13 14 15
0.7965122243 0.7544195589 0.4766985810 0.0853034707 0.0253034707
16 17 18 19 20
0.0007454266 0.0275824927 0.0407454266 -0.0529287058 -0.0288818601
21 22 23 24 25
0.1779552060 0.2607454266 0.3193024449 0.5227902206 0.2450692426
26 27 28 29 30
0.2422790221 0.1594888015 0.3577209779 0.4959531544 0.6319063088
31 32 33 34 35
0.7610223970 0.6850692426 0.5805111985 0.3633014190 -0.4465122243
36 37 38 39 40
-0.5658146692 -0.6207454266 -0.2121405368 0.0278594632 -0.0966985810
41 42 43 44 45
-0.2698615148 -0.3339083604 -0.1820020516 -0.0351649855 -0.0283279193
46 47 48 49 50
-0.1025132501 -0.6693024449 -0.7113951103 -0.7549307574 -0.2893503163
51 52 53 54 55
-0.0379552060 -0.0739083604 -0.4040468456 -0.7136741323 -1.1301384852
56 57 58
-1.2833014190 -1.3764643529 -1.5392545734
> postscript(file="/var/www/html/rcomp/tmp/6j8181258661063.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 = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 0.6539083604 NA
1 0.1739083604 0.6539083604
2 -0.1746965293 0.1739083604
3 -0.1878594632 -0.1746965293
4 0.1503727133 -0.1878594632
5 0.3749307574 0.1503727133
6 0.6040468456 0.3749307574
7 0.6622790221 0.6040468456
8 0.6463258677 0.6622790221
9 1.0177209779 0.6463258677
10 0.7965122243 1.0177209779
11 0.7544195589 0.7965122243
12 0.4766985810 0.7544195589
13 0.0853034707 0.4766985810
14 0.0253034707 0.0853034707
15 0.0007454266 0.0253034707
16 0.0275824927 0.0007454266
17 0.0407454266 0.0275824927
18 -0.0529287058 0.0407454266
19 -0.0288818601 -0.0529287058
20 0.1779552060 -0.0288818601
21 0.2607454266 0.1779552060
22 0.3193024449 0.2607454266
23 0.5227902206 0.3193024449
24 0.2450692426 0.5227902206
25 0.2422790221 0.2450692426
26 0.1594888015 0.2422790221
27 0.3577209779 0.1594888015
28 0.4959531544 0.3577209779
29 0.6319063088 0.4959531544
30 0.7610223970 0.6319063088
31 0.6850692426 0.7610223970
32 0.5805111985 0.6850692426
33 0.3633014190 0.5805111985
34 -0.4465122243 0.3633014190
35 -0.5658146692 -0.4465122243
36 -0.6207454266 -0.5658146692
37 -0.2121405368 -0.6207454266
38 0.0278594632 -0.2121405368
39 -0.0966985810 0.0278594632
40 -0.2698615148 -0.0966985810
41 -0.3339083604 -0.2698615148
42 -0.1820020516 -0.3339083604
43 -0.0351649855 -0.1820020516
44 -0.0283279193 -0.0351649855
45 -0.1025132501 -0.0283279193
46 -0.6693024449 -0.1025132501
47 -0.7113951103 -0.6693024449
48 -0.7549307574 -0.7113951103
49 -0.2893503163 -0.7549307574
50 -0.0379552060 -0.2893503163
51 -0.0739083604 -0.0379552060
52 -0.4040468456 -0.0739083604
53 -0.7136741323 -0.4040468456
54 -1.1301384852 -0.7136741323
55 -1.2833014190 -1.1301384852
56 -1.3764643529 -1.2833014190
57 -1.5392545734 -1.3764643529
58 NA -1.5392545734
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.1739083604 0.6539083604
[2,] -0.1746965293 0.1739083604
[3,] -0.1878594632 -0.1746965293
[4,] 0.1503727133 -0.1878594632
[5,] 0.3749307574 0.1503727133
[6,] 0.6040468456 0.3749307574
[7,] 0.6622790221 0.6040468456
[8,] 0.6463258677 0.6622790221
[9,] 1.0177209779 0.6463258677
[10,] 0.7965122243 1.0177209779
[11,] 0.7544195589 0.7965122243
[12,] 0.4766985810 0.7544195589
[13,] 0.0853034707 0.4766985810
[14,] 0.0253034707 0.0853034707
[15,] 0.0007454266 0.0253034707
[16,] 0.0275824927 0.0007454266
[17,] 0.0407454266 0.0275824927
[18,] -0.0529287058 0.0407454266
[19,] -0.0288818601 -0.0529287058
[20,] 0.1779552060 -0.0288818601
[21,] 0.2607454266 0.1779552060
[22,] 0.3193024449 0.2607454266
[23,] 0.5227902206 0.3193024449
[24,] 0.2450692426 0.5227902206
[25,] 0.2422790221 0.2450692426
[26,] 0.1594888015 0.2422790221
[27,] 0.3577209779 0.1594888015
[28,] 0.4959531544 0.3577209779
[29,] 0.6319063088 0.4959531544
[30,] 0.7610223970 0.6319063088
[31,] 0.6850692426 0.7610223970
[32,] 0.5805111985 0.6850692426
[33,] 0.3633014190 0.5805111985
[34,] -0.4465122243 0.3633014190
[35,] -0.5658146692 -0.4465122243
[36,] -0.6207454266 -0.5658146692
[37,] -0.2121405368 -0.6207454266
[38,] 0.0278594632 -0.2121405368
[39,] -0.0966985810 0.0278594632
[40,] -0.2698615148 -0.0966985810
[41,] -0.3339083604 -0.2698615148
[42,] -0.1820020516 -0.3339083604
[43,] -0.0351649855 -0.1820020516
[44,] -0.0283279193 -0.0351649855
[45,] -0.1025132501 -0.0283279193
[46,] -0.6693024449 -0.1025132501
[47,] -0.7113951103 -0.6693024449
[48,] -0.7549307574 -0.7113951103
[49,] -0.2893503163 -0.7549307574
[50,] -0.0379552060 -0.2893503163
[51,] -0.0739083604 -0.0379552060
[52,] -0.4040468456 -0.0739083604
[53,] -0.7136741323 -0.4040468456
[54,] -1.1301384852 -0.7136741323
[55,] -1.2833014190 -1.1301384852
[56,] -1.3764643529 -1.2833014190
[57,] -1.5392545734 -1.3764643529
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.1739083604 0.6539083604
2 -0.1746965293 0.1739083604
3 -0.1878594632 -0.1746965293
4 0.1503727133 -0.1878594632
5 0.3749307574 0.1503727133
6 0.6040468456 0.3749307574
7 0.6622790221 0.6040468456
8 0.6463258677 0.6622790221
9 1.0177209779 0.6463258677
10 0.7965122243 1.0177209779
11 0.7544195589 0.7965122243
12 0.4766985810 0.7544195589
13 0.0853034707 0.4766985810
14 0.0253034707 0.0853034707
15 0.0007454266 0.0253034707
16 0.0275824927 0.0007454266
17 0.0407454266 0.0275824927
18 -0.0529287058 0.0407454266
19 -0.0288818601 -0.0529287058
20 0.1779552060 -0.0288818601
21 0.2607454266 0.1779552060
22 0.3193024449 0.2607454266
23 0.5227902206 0.3193024449
24 0.2450692426 0.5227902206
25 0.2422790221 0.2450692426
26 0.1594888015 0.2422790221
27 0.3577209779 0.1594888015
28 0.4959531544 0.3577209779
29 0.6319063088 0.4959531544
30 0.7610223970 0.6319063088
31 0.6850692426 0.7610223970
32 0.5805111985 0.6850692426
33 0.3633014190 0.5805111985
34 -0.4465122243 0.3633014190
35 -0.5658146692 -0.4465122243
36 -0.6207454266 -0.5658146692
37 -0.2121405368 -0.6207454266
38 0.0278594632 -0.2121405368
39 -0.0966985810 0.0278594632
40 -0.2698615148 -0.0966985810
41 -0.3339083604 -0.2698615148
42 -0.1820020516 -0.3339083604
43 -0.0351649855 -0.1820020516
44 -0.0283279193 -0.0351649855
45 -0.1025132501 -0.0283279193
46 -0.6693024449 -0.1025132501
47 -0.7113951103 -0.6693024449
48 -0.7549307574 -0.7113951103
49 -0.2893503163 -0.7549307574
50 -0.0379552060 -0.2893503163
51 -0.0739083604 -0.0379552060
52 -0.4040468456 -0.0739083604
53 -0.7136741323 -0.4040468456
54 -1.1301384852 -0.7136741323
55 -1.2833014190 -1.1301384852
56 -1.3764643529 -1.2833014190
57 -1.5392545734 -1.3764643529
> 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/787tt1258661063.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/8ppjy1258661063.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/9e4o51258661063.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/10l9e21258661063.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/11b74k1258661063.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/12746k1258661063.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/13ig4g1258661063.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/14xus31258661063.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/15lvoo1258661063.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/16edk31258661063.tab")
+ }
>
> system("convert tmp/1y7ix1258661063.ps tmp/1y7ix1258661063.png")
> system("convert tmp/2l6nh1258661063.ps tmp/2l6nh1258661063.png")
> system("convert tmp/33c521258661063.ps tmp/33c521258661063.png")
> system("convert tmp/4hiiq1258661063.ps tmp/4hiiq1258661063.png")
> system("convert tmp/58kkl1258661063.ps tmp/58kkl1258661063.png")
> system("convert tmp/6j8181258661063.ps tmp/6j8181258661063.png")
> system("convert tmp/787tt1258661063.ps tmp/787tt1258661063.png")
> system("convert tmp/8ppjy1258661063.ps tmp/8ppjy1258661063.png")
> system("convert tmp/9e4o51258661063.ps tmp/9e4o51258661063.png")
> system("convert tmp/10l9e21258661063.ps tmp/10l9e21258661063.png")
>
>
> proc.time()
user system elapsed
2.338 1.541 2.979