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,8.6,8.9,8.5,8.9,8.3,8.9,7.8,9,7.8,9,8,9,8.6,9,8.9,9,8.9,9,8.6,9,8.3,9.1,8.3,9,8.3,9.1,8.4,9.1,8.5,9,8.4,9,8.6,9,8.5,9,8.5,8.9,8.4,8.9,8.5,8.9,8.5,8.9,8.5,8.8,8.5,8.8,8.5,8.7,8.5,8.7,8.5,8.5,8.5,8.5,8.6,8.4,8.4,8.2,8.1,8.2,8,8.1,8,8.1,8,8,8,7.9,7.9,7.8,7.8,7.7,7.8,7.6,7.9,7.5,8.1,7.5,8,7.5,7.6,7.5,7.3,7.5,7,7.4,6.8,7.4,7,7.3,7.1,7.3,7.2,7.3,7.1,7.2,6.9,7.2,6.7,7.3,6.7,7.4,6.6,7.4,6.9,7.5,7.3,7.6,7.5,7.7,7.3,7.9,7.1,8,6.9,8.2,7.1),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 8.9 8.6 1 0 0 0 0 0 0 0 0 0 0 1
2 8.9 8.5 0 1 0 0 0 0 0 0 0 0 0 2
3 8.9 8.3 0 0 1 0 0 0 0 0 0 0 0 3
4 8.9 7.8 0 0 0 1 0 0 0 0 0 0 0 4
5 9.0 7.8 0 0 0 0 1 0 0 0 0 0 0 5
6 9.0 8.0 0 0 0 0 0 1 0 0 0 0 0 6
7 9.0 8.6 0 0 0 0 0 0 1 0 0 0 0 7
8 9.0 8.9 0 0 0 0 0 0 0 1 0 0 0 8
9 9.0 8.9 0 0 0 0 0 0 0 0 1 0 0 9
10 9.0 8.6 0 0 0 0 0 0 0 0 0 1 0 10
11 9.0 8.3 0 0 0 0 0 0 0 0 0 0 1 11
12 9.1 8.3 0 0 0 0 0 0 0 0 0 0 0 12
13 9.0 8.3 1 0 0 0 0 0 0 0 0 0 0 13
14 9.1 8.4 0 1 0 0 0 0 0 0 0 0 0 14
15 9.1 8.5 0 0 1 0 0 0 0 0 0 0 0 15
16 9.0 8.4 0 0 0 1 0 0 0 0 0 0 0 16
17 9.0 8.6 0 0 0 0 1 0 0 0 0 0 0 17
18 9.0 8.5 0 0 0 0 0 1 0 0 0 0 0 18
19 9.0 8.5 0 0 0 0 0 0 1 0 0 0 0 19
20 8.9 8.4 0 0 0 0 0 0 0 1 0 0 0 20
21 8.9 8.5 0 0 0 0 0 0 0 0 1 0 0 21
22 8.9 8.5 0 0 0 0 0 0 0 0 0 1 0 22
23 8.9 8.5 0 0 0 0 0 0 0 0 0 0 1 23
24 8.8 8.5 0 0 0 0 0 0 0 0 0 0 0 24
25 8.8 8.5 1 0 0 0 0 0 0 0 0 0 0 25
26 8.7 8.5 0 1 0 0 0 0 0 0 0 0 0 26
27 8.7 8.5 0 0 1 0 0 0 0 0 0 0 0 27
28 8.5 8.5 0 0 0 1 0 0 0 0 0 0 0 28
29 8.5 8.6 0 0 0 0 1 0 0 0 0 0 0 29
30 8.4 8.4 0 0 0 0 0 1 0 0 0 0 0 30
31 8.2 8.1 0 0 0 0 0 0 1 0 0 0 0 31
32 8.2 8.0 0 0 0 0 0 0 0 1 0 0 0 32
33 8.1 8.0 0 0 0 0 0 0 0 0 1 0 0 33
34 8.1 8.0 0 0 0 0 0 0 0 0 0 1 0 34
35 8.0 8.0 0 0 0 0 0 0 0 0 0 0 1 35
36 7.9 7.9 0 0 0 0 0 0 0 0 0 0 0 36
37 7.8 7.8 1 0 0 0 0 0 0 0 0 0 0 37
38 7.7 7.8 0 1 0 0 0 0 0 0 0 0 0 38
39 7.6 7.9 0 0 1 0 0 0 0 0 0 0 0 39
40 7.5 8.1 0 0 0 1 0 0 0 0 0 0 0 40
41 7.5 8.0 0 0 0 0 1 0 0 0 0 0 0 41
42 7.5 7.6 0 0 0 0 0 1 0 0 0 0 0 42
43 7.5 7.3 0 0 0 0 0 0 1 0 0 0 0 43
44 7.5 7.0 0 0 0 0 0 0 0 1 0 0 0 44
45 7.4 6.8 0 0 0 0 0 0 0 0 1 0 0 45
46 7.4 7.0 0 0 0 0 0 0 0 0 0 1 0 46
47 7.3 7.1 0 0 0 0 0 0 0 0 0 0 1 47
48 7.3 7.2 0 0 0 0 0 0 0 0 0 0 0 48
49 7.3 7.1 1 0 0 0 0 0 0 0 0 0 0 49
50 7.2 6.9 0 1 0 0 0 0 0 0 0 0 0 50
51 7.2 6.7 0 0 1 0 0 0 0 0 0 0 0 51
52 7.3 6.7 0 0 0 1 0 0 0 0 0 0 0 52
53 7.4 6.6 0 0 0 0 1 0 0 0 0 0 0 53
54 7.4 6.9 0 0 0 0 0 1 0 0 0 0 0 54
55 7.5 7.3 0 0 0 0 0 0 1 0 0 0 0 55
56 7.6 7.5 0 0 0 0 0 0 0 1 0 0 0 56
57 7.7 7.3 0 0 0 0 0 0 0 0 1 0 0 57
58 7.9 7.1 0 0 0 0 0 0 0 0 0 1 0 58
59 8.0 6.9 0 0 0 0 0 0 0 0 0 0 1 59
60 8.2 7.1 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
6.31542 0.36242 -0.26381 -0.26481 -0.24581 -0.25230
M5 M6 M7 M8 M9 M10
-0.19504 -0.17604 -0.20053 -0.17602 -0.14977 -0.06351
M11 t
-0.03001 -0.02451
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.51839 -0.17087 -0.02517 0.18124 0.78188
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.315421 0.969361 6.515 4.88e-08 ***
X 0.362418 0.107132 3.383 0.00147 **
M1 -0.263814 0.194209 -1.358 0.18096
M2 -0.264810 0.193965 -1.365 0.17882
M3 -0.245805 0.193753 -1.269 0.21095
M4 -0.252304 0.193881 -1.301 0.19962
M5 -0.195045 0.193296 -1.009 0.31823
M6 -0.176040 0.193169 -0.911 0.36687
M7 -0.200526 0.192673 -1.041 0.30343
M8 -0.176018 0.192584 -0.914 0.36549
M9 -0.149765 0.192452 -0.778 0.44044
M10 -0.063512 0.192396 -0.330 0.74281
M11 -0.030011 0.192491 -0.156 0.87679
t -0.024508 0.004069 -6.023 2.68e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3041 on 46 degrees of freedom
Multiple R-squared: 0.8458, Adjusted R-squared: 0.8023
F-statistic: 19.41 on 13 and 46 DF, p-value: 1.926e-14
> 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,] 5.654604e-03 0.0113092090 0.9943454
[2,] 1.455921e-03 0.0029118412 0.9985441
[3,] 8.334932e-04 0.0016669865 0.9991665
[4,] 9.767254e-04 0.0019534507 0.9990233
[5,] 4.400092e-04 0.0008800185 0.9995600
[6,] 1.857475e-04 0.0003714949 0.9998143
[7,] 8.734737e-05 0.0001746947 0.9999127
[8,] 2.822441e-04 0.0005644882 0.9997178
[9,] 1.898273e-04 0.0003796546 0.9998102
[10,] 3.577569e-04 0.0007155138 0.9996422
[11,] 5.857363e-04 0.0011714727 0.9994143
[12,] 2.254933e-03 0.0045098654 0.9977451
[13,] 5.233853e-03 0.0104677063 0.9947661
[14,] 1.706396e-02 0.0341279181 0.9829360
[15,] 6.933055e-02 0.1386611051 0.9306694
[16,] 1.099684e-01 0.2199367590 0.8900316
[17,] 1.316919e-01 0.2633838265 0.8683081
[18,] 1.405343e-01 0.2810685720 0.8594657
[19,] 1.602773e-01 0.3205545416 0.8397227
[20,] 1.671072e-01 0.3342144594 0.8328928
[21,] 1.955486e-01 0.3910971481 0.8044514
[22,] 2.399498e-01 0.4798996990 0.7600502
[23,] 2.661128e-01 0.5322256870 0.7338872
[24,] 2.787056e-01 0.5574112321 0.7212944
[25,] 2.457225e-01 0.4914450480 0.7542775
[26,] 5.017617e-01 0.9964765460 0.4982383
[27,] 8.972930e-01 0.2054140815 0.1027070
> postscript(file="/var/www/html/rcomp/tmp/1makf1258709658.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/2cp3d1258709658.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/3x6fr1258709658.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/4x09g1258709658.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/5ln021258709658.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 = 60
Frequency = 1
1 2 3 4 5 6
-0.243892268 -0.182147201 -0.104160357 0.108055107 0.175303463 0.108323197
7 8 9 10 11 12
-0.060134045 -0.168859377 -0.170604443 -0.123624178 -0.023892268 0.070604443
13 14 15 16 17 18
0.258926399 0.348187911 0.317449423 0.284697778 0.179462579 0.221207646
19 20 21 22 23 24
0.270201067 0.206442845 0.168456001 0.106710935 0.097717513 -0.007785776
25 26 27 28 29 30
0.280536180 0.206039469 0.211542758 0.042549336 -0.026444086 -0.048457242
31 32 33 34 35 36
-0.090738488 -0.054496711 -0.156241777 -0.217986844 -0.326980266 -0.396241777
37 38 39 40 41 42
-0.171678044 -0.246174755 -0.376913243 -0.518390220 -0.514900087 -0.364429688
43 44 45 46 47 48
-0.206710935 -0.097985603 -0.127247114 -0.261475735 -0.406710935 -0.448456001
49 50 51 52 53 54
-0.123892268 -0.125905424 -0.047918580 0.083087998 0.186578131 0.083356088
55 56 57 58 59 60
0.087382400 0.114898846 0.285637334 0.496375822 0.659865955 0.781879111
> postscript(file="/var/www/html/rcomp/tmp/6xhw91258709658.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.243892268 NA
1 -0.182147201 -0.243892268
2 -0.104160357 -0.182147201
3 0.108055107 -0.104160357
4 0.175303463 0.108055107
5 0.108323197 0.175303463
6 -0.060134045 0.108323197
7 -0.168859377 -0.060134045
8 -0.170604443 -0.168859377
9 -0.123624178 -0.170604443
10 -0.023892268 -0.123624178
11 0.070604443 -0.023892268
12 0.258926399 0.070604443
13 0.348187911 0.258926399
14 0.317449423 0.348187911
15 0.284697778 0.317449423
16 0.179462579 0.284697778
17 0.221207646 0.179462579
18 0.270201067 0.221207646
19 0.206442845 0.270201067
20 0.168456001 0.206442845
21 0.106710935 0.168456001
22 0.097717513 0.106710935
23 -0.007785776 0.097717513
24 0.280536180 -0.007785776
25 0.206039469 0.280536180
26 0.211542758 0.206039469
27 0.042549336 0.211542758
28 -0.026444086 0.042549336
29 -0.048457242 -0.026444086
30 -0.090738488 -0.048457242
31 -0.054496711 -0.090738488
32 -0.156241777 -0.054496711
33 -0.217986844 -0.156241777
34 -0.326980266 -0.217986844
35 -0.396241777 -0.326980266
36 -0.171678044 -0.396241777
37 -0.246174755 -0.171678044
38 -0.376913243 -0.246174755
39 -0.518390220 -0.376913243
40 -0.514900087 -0.518390220
41 -0.364429688 -0.514900087
42 -0.206710935 -0.364429688
43 -0.097985603 -0.206710935
44 -0.127247114 -0.097985603
45 -0.261475735 -0.127247114
46 -0.406710935 -0.261475735
47 -0.448456001 -0.406710935
48 -0.123892268 -0.448456001
49 -0.125905424 -0.123892268
50 -0.047918580 -0.125905424
51 0.083087998 -0.047918580
52 0.186578131 0.083087998
53 0.083356088 0.186578131
54 0.087382400 0.083356088
55 0.114898846 0.087382400
56 0.285637334 0.114898846
57 0.496375822 0.285637334
58 0.659865955 0.496375822
59 0.781879111 0.659865955
60 NA 0.781879111
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.182147201 -0.243892268
[2,] -0.104160357 -0.182147201
[3,] 0.108055107 -0.104160357
[4,] 0.175303463 0.108055107
[5,] 0.108323197 0.175303463
[6,] -0.060134045 0.108323197
[7,] -0.168859377 -0.060134045
[8,] -0.170604443 -0.168859377
[9,] -0.123624178 -0.170604443
[10,] -0.023892268 -0.123624178
[11,] 0.070604443 -0.023892268
[12,] 0.258926399 0.070604443
[13,] 0.348187911 0.258926399
[14,] 0.317449423 0.348187911
[15,] 0.284697778 0.317449423
[16,] 0.179462579 0.284697778
[17,] 0.221207646 0.179462579
[18,] 0.270201067 0.221207646
[19,] 0.206442845 0.270201067
[20,] 0.168456001 0.206442845
[21,] 0.106710935 0.168456001
[22,] 0.097717513 0.106710935
[23,] -0.007785776 0.097717513
[24,] 0.280536180 -0.007785776
[25,] 0.206039469 0.280536180
[26,] 0.211542758 0.206039469
[27,] 0.042549336 0.211542758
[28,] -0.026444086 0.042549336
[29,] -0.048457242 -0.026444086
[30,] -0.090738488 -0.048457242
[31,] -0.054496711 -0.090738488
[32,] -0.156241777 -0.054496711
[33,] -0.217986844 -0.156241777
[34,] -0.326980266 -0.217986844
[35,] -0.396241777 -0.326980266
[36,] -0.171678044 -0.396241777
[37,] -0.246174755 -0.171678044
[38,] -0.376913243 -0.246174755
[39,] -0.518390220 -0.376913243
[40,] -0.514900087 -0.518390220
[41,] -0.364429688 -0.514900087
[42,] -0.206710935 -0.364429688
[43,] -0.097985603 -0.206710935
[44,] -0.127247114 -0.097985603
[45,] -0.261475735 -0.127247114
[46,] -0.406710935 -0.261475735
[47,] -0.448456001 -0.406710935
[48,] -0.123892268 -0.448456001
[49,] -0.125905424 -0.123892268
[50,] -0.047918580 -0.125905424
[51,] 0.083087998 -0.047918580
[52,] 0.186578131 0.083087998
[53,] 0.083356088 0.186578131
[54,] 0.087382400 0.083356088
[55,] 0.114898846 0.087382400
[56,] 0.285637334 0.114898846
[57,] 0.496375822 0.285637334
[58,] 0.659865955 0.496375822
[59,] 0.781879111 0.659865955
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.182147201 -0.243892268
2 -0.104160357 -0.182147201
3 0.108055107 -0.104160357
4 0.175303463 0.108055107
5 0.108323197 0.175303463
6 -0.060134045 0.108323197
7 -0.168859377 -0.060134045
8 -0.170604443 -0.168859377
9 -0.123624178 -0.170604443
10 -0.023892268 -0.123624178
11 0.070604443 -0.023892268
12 0.258926399 0.070604443
13 0.348187911 0.258926399
14 0.317449423 0.348187911
15 0.284697778 0.317449423
16 0.179462579 0.284697778
17 0.221207646 0.179462579
18 0.270201067 0.221207646
19 0.206442845 0.270201067
20 0.168456001 0.206442845
21 0.106710935 0.168456001
22 0.097717513 0.106710935
23 -0.007785776 0.097717513
24 0.280536180 -0.007785776
25 0.206039469 0.280536180
26 0.211542758 0.206039469
27 0.042549336 0.211542758
28 -0.026444086 0.042549336
29 -0.048457242 -0.026444086
30 -0.090738488 -0.048457242
31 -0.054496711 -0.090738488
32 -0.156241777 -0.054496711
33 -0.217986844 -0.156241777
34 -0.326980266 -0.217986844
35 -0.396241777 -0.326980266
36 -0.171678044 -0.396241777
37 -0.246174755 -0.171678044
38 -0.376913243 -0.246174755
39 -0.518390220 -0.376913243
40 -0.514900087 -0.518390220
41 -0.364429688 -0.514900087
42 -0.206710935 -0.364429688
43 -0.097985603 -0.206710935
44 -0.127247114 -0.097985603
45 -0.261475735 -0.127247114
46 -0.406710935 -0.261475735
47 -0.448456001 -0.406710935
48 -0.123892268 -0.448456001
49 -0.125905424 -0.123892268
50 -0.047918580 -0.125905424
51 0.083087998 -0.047918580
52 0.186578131 0.083087998
53 0.083356088 0.186578131
54 0.087382400 0.083356088
55 0.114898846 0.087382400
56 0.285637334 0.114898846
57 0.496375822 0.285637334
58 0.659865955 0.496375822
59 0.781879111 0.659865955
> 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/7o7zl1258709658.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/86vm81258709658.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/9y5a11258709658.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/10wtkl1258709658.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/116ged1258709658.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/120vnq1258709658.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/13zyx11258709658.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/1487rb1258709658.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/158ren1258709658.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/16a2hl1258709658.tab")
+ }
>
> system("convert tmp/1makf1258709658.ps tmp/1makf1258709658.png")
> system("convert tmp/2cp3d1258709658.ps tmp/2cp3d1258709658.png")
> system("convert tmp/3x6fr1258709658.ps tmp/3x6fr1258709658.png")
> system("convert tmp/4x09g1258709658.ps tmp/4x09g1258709658.png")
> system("convert tmp/5ln021258709658.ps tmp/5ln021258709658.png")
> system("convert tmp/6xhw91258709658.ps tmp/6xhw91258709658.png")
> system("convert tmp/7o7zl1258709658.ps tmp/7o7zl1258709658.png")
> system("convert tmp/86vm81258709658.ps tmp/86vm81258709658.png")
> system("convert tmp/9y5a11258709658.ps tmp/9y5a11258709658.png")
> system("convert tmp/10wtkl1258709658.ps tmp/10wtkl1258709658.png")
>
>
> proc.time()
user system elapsed
2.343 1.528 3.149