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(2.97
+ ,101.1
+ ,2.98
+ ,3.01
+ ,3.06
+ ,3.12
+ ,3.58
+ ,3.02
+ ,100.93
+ ,2.97
+ ,2.98
+ ,3.01
+ ,3.06
+ ,3.52
+ ,3.07
+ ,100.85
+ ,3.02
+ ,2.97
+ ,2.98
+ ,3.01
+ ,3.45
+ ,3.18
+ ,100.93
+ ,3.07
+ ,3.02
+ ,2.97
+ ,2.98
+ ,3.36
+ ,3.29
+ ,99.6
+ ,3.18
+ ,3.07
+ ,3.02
+ ,2.97
+ ,3.27
+ ,3.43
+ ,101.88
+ ,3.29
+ ,3.18
+ ,3.07
+ ,3.02
+ ,3.21
+ ,3.61
+ ,101.81
+ ,3.43
+ ,3.29
+ ,3.18
+ ,3.07
+ ,3.19
+ ,3.74
+ ,102.38
+ ,3.61
+ ,3.43
+ ,3.29
+ ,3.18
+ ,3.16
+ ,3.87
+ ,102.74
+ ,3.74
+ ,3.61
+ ,3.43
+ ,3.29
+ ,3.12
+ ,3.88
+ ,102.82
+ ,3.87
+ ,3.74
+ ,3.61
+ ,3.43
+ ,3.06
+ ,4.09
+ ,101.72
+ ,3.88
+ ,3.87
+ ,3.74
+ ,3.61
+ ,3.01
+ ,4.19
+ ,103.47
+ ,4.09
+ ,3.88
+ ,3.87
+ ,3.74
+ ,2.98
+ ,4.2
+ ,102.98
+ ,4.19
+ ,4.09
+ ,3.88
+ ,3.87
+ ,2.97
+ ,4.29
+ ,102.68
+ ,4.2
+ ,4.19
+ ,4.09
+ ,3.88
+ ,3.02
+ ,4.37
+ ,102.9
+ ,4.29
+ ,4.2
+ ,4.19
+ ,4.09
+ ,3.07
+ ,4.47
+ ,103.03
+ ,4.37
+ ,4.29
+ ,4.2
+ ,4.19
+ ,3.18
+ ,4.61
+ ,101.29
+ ,4.47
+ ,4.37
+ ,4.29
+ ,4.2
+ ,3.29
+ ,4.65
+ ,103.69
+ ,4.61
+ ,4.47
+ ,4.37
+ ,4.29
+ ,3.43
+ ,4.69
+ ,103.68
+ ,4.65
+ ,4.61
+ ,4.47
+ ,4.37
+ ,3.61
+ ,4.82
+ ,104.2
+ ,4.69
+ ,4.65
+ ,4.61
+ ,4.47
+ ,3.74
+ ,4.86
+ ,104.08
+ ,4.82
+ ,4.69
+ ,4.65
+ ,4.61
+ ,3.87
+ ,4.87
+ ,104.16
+ ,4.86
+ ,4.82
+ ,4.69
+ ,4.65
+ ,3.88
+ ,5.01
+ ,103.05
+ ,4.87
+ ,4.86
+ ,4.82
+ ,4.69
+ ,4.09
+ ,5.03
+ ,104.66
+ ,5.01
+ ,4.87
+ ,4.86
+ ,4.82
+ ,4.19
+ ,5.13
+ ,104.46
+ ,5.03
+ ,5.01
+ ,4.87
+ ,4.86
+ ,4.2
+ ,5.18
+ ,104.95
+ ,5.13
+ ,5.03
+ ,5.01
+ ,4.87
+ ,4.29
+ ,5.21
+ ,105.85
+ ,5.18
+ ,5.13
+ ,5.03
+ ,5.01
+ ,4.37
+ ,5.26
+ ,106.23
+ ,5.21
+ ,5.18
+ ,5.13
+ ,5.03
+ ,4.47
+ ,5.25
+ ,104.86
+ ,5.26
+ ,5.21
+ ,5.18
+ ,5.13
+ ,4.61
+ ,5.2
+ ,107.44
+ ,5.25
+ ,5.26
+ ,5.21
+ ,5.18
+ ,4.65
+ ,5.16
+ ,108.23
+ ,5.2
+ ,5.25
+ ,5.26
+ ,5.21
+ ,4.69
+ ,5.19
+ ,108.45
+ ,5.16
+ ,5.2
+ ,5.25
+ ,5.26
+ ,4.82
+ ,5.39
+ ,109.39
+ ,5.19
+ ,5.16
+ ,5.2
+ ,5.25
+ ,4.86
+ ,5.58
+ ,110.15
+ ,5.39
+ ,5.19
+ ,5.16
+ ,5.2
+ ,4.87
+ ,5.76
+ ,109.13
+ ,5.58
+ ,5.39
+ ,5.19
+ ,5.16
+ ,5.01
+ ,5.89
+ ,110.28
+ ,5.76
+ ,5.58
+ ,5.39
+ ,5.19
+ ,5.03
+ ,5.98
+ ,110.17
+ ,5.89
+ ,5.76
+ ,5.58
+ ,5.39
+ ,5.13
+ ,6.02
+ ,109.99
+ ,5.98
+ ,5.89
+ ,5.76
+ ,5.58
+ ,5.18
+ ,5.62
+ ,109.26
+ ,6.02
+ ,5.98
+ ,5.89
+ ,5.76
+ ,5.21
+ ,4.87
+ ,109.11
+ ,5.62
+ ,6.02
+ ,5.98
+ ,5.89
+ ,5.26
+ ,4.24
+ ,107.06
+ ,4.87
+ ,5.62
+ ,6.02
+ ,5.98
+ ,5.25
+ ,4.02
+ ,109.53
+ ,4.24
+ ,4.87
+ ,5.62
+ ,6.02
+ ,5.2
+ ,3.74
+ ,108.92
+ ,4.02
+ ,4.24
+ ,4.87
+ ,5.62
+ ,5.16
+ ,3.45
+ ,109.24
+ ,3.74
+ ,4.02
+ ,4.24
+ ,4.87
+ ,5.19
+ ,3.34
+ ,109.12
+ ,3.45
+ ,3.74
+ ,4.02
+ ,4.24
+ ,5.39
+ ,3.21
+ ,109
+ ,3.34
+ ,3.45
+ ,3.74
+ ,4.02
+ ,5.58
+ ,3.12
+ ,107.23
+ ,3.21
+ ,3.34
+ ,3.45
+ ,3.74
+ ,5.76
+ ,3.04
+ ,109.49
+ ,3.12
+ ,3.21
+ ,3.34
+ ,3.45
+ ,5.89)
+ ,dim=c(7
+ ,48)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'y1'
+ ,'y2'
+ ,'y3'
+ ,'y4'
+ ,'y12')
+ ,1:48))
> y <- array(NA,dim=c(7,48),dimnames=list(c('Y','X','y1','y2','y3','y4','y12'),1:48))
> 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 y1 y2 y3 y4 y12 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 2.97 101.10 2.98 3.01 3.06 3.12 3.58 1 0 0 0 0 0 0 0 0 0 0 1
2 3.02 100.93 2.97 2.98 3.01 3.06 3.52 0 1 0 0 0 0 0 0 0 0 0 2
3 3.07 100.85 3.02 2.97 2.98 3.01 3.45 0 0 1 0 0 0 0 0 0 0 0 3
4 3.18 100.93 3.07 3.02 2.97 2.98 3.36 0 0 0 1 0 0 0 0 0 0 0 4
5 3.29 99.60 3.18 3.07 3.02 2.97 3.27 0 0 0 0 1 0 0 0 0 0 0 5
6 3.43 101.88 3.29 3.18 3.07 3.02 3.21 0 0 0 0 0 1 0 0 0 0 0 6
7 3.61 101.81 3.43 3.29 3.18 3.07 3.19 0 0 0 0 0 0 1 0 0 0 0 7
8 3.74 102.38 3.61 3.43 3.29 3.18 3.16 0 0 0 0 0 0 0 1 0 0 0 8
9 3.87 102.74 3.74 3.61 3.43 3.29 3.12 0 0 0 0 0 0 0 0 1 0 0 9
10 3.88 102.82 3.87 3.74 3.61 3.43 3.06 0 0 0 0 0 0 0 0 0 1 0 10
11 4.09 101.72 3.88 3.87 3.74 3.61 3.01 0 0 0 0 0 0 0 0 0 0 1 11
12 4.19 103.47 4.09 3.88 3.87 3.74 2.98 0 0 0 0 0 0 0 0 0 0 0 12
13 4.20 102.98 4.19 4.09 3.88 3.87 2.97 1 0 0 0 0 0 0 0 0 0 0 13
14 4.29 102.68 4.20 4.19 4.09 3.88 3.02 0 1 0 0 0 0 0 0 0 0 0 14
15 4.37 102.90 4.29 4.20 4.19 4.09 3.07 0 0 1 0 0 0 0 0 0 0 0 15
16 4.47 103.03 4.37 4.29 4.20 4.19 3.18 0 0 0 1 0 0 0 0 0 0 0 16
17 4.61 101.29 4.47 4.37 4.29 4.20 3.29 0 0 0 0 1 0 0 0 0 0 0 17
18 4.65 103.69 4.61 4.47 4.37 4.29 3.43 0 0 0 0 0 1 0 0 0 0 0 18
19 4.69 103.68 4.65 4.61 4.47 4.37 3.61 0 0 0 0 0 0 1 0 0 0 0 19
20 4.82 104.20 4.69 4.65 4.61 4.47 3.74 0 0 0 0 0 0 0 1 0 0 0 20
21 4.86 104.08 4.82 4.69 4.65 4.61 3.87 0 0 0 0 0 0 0 0 1 0 0 21
22 4.87 104.16 4.86 4.82 4.69 4.65 3.88 0 0 0 0 0 0 0 0 0 1 0 22
23 5.01 103.05 4.87 4.86 4.82 4.69 4.09 0 0 0 0 0 0 0 0 0 0 1 23
24 5.03 104.66 5.01 4.87 4.86 4.82 4.19 0 0 0 0 0 0 0 0 0 0 0 24
25 5.13 104.46 5.03 5.01 4.87 4.86 4.20 1 0 0 0 0 0 0 0 0 0 0 25
26 5.18 104.95 5.13 5.03 5.01 4.87 4.29 0 1 0 0 0 0 0 0 0 0 0 26
27 5.21 105.85 5.18 5.13 5.03 5.01 4.37 0 0 1 0 0 0 0 0 0 0 0 27
28 5.26 106.23 5.21 5.18 5.13 5.03 4.47 0 0 0 1 0 0 0 0 0 0 0 28
29 5.25 104.86 5.26 5.21 5.18 5.13 4.61 0 0 0 0 1 0 0 0 0 0 0 29
30 5.20 107.44 5.25 5.26 5.21 5.18 4.65 0 0 0 0 0 1 0 0 0 0 0 30
31 5.16 108.23 5.20 5.25 5.26 5.21 4.69 0 0 0 0 0 0 1 0 0 0 0 31
32 5.19 108.45 5.16 5.20 5.25 5.26 4.82 0 0 0 0 0 0 0 1 0 0 0 32
33 5.39 109.39 5.19 5.16 5.20 5.25 4.86 0 0 0 0 0 0 0 0 1 0 0 33
34 5.58 110.15 5.39 5.19 5.16 5.20 4.87 0 0 0 0 0 0 0 0 0 1 0 34
35 5.76 109.13 5.58 5.39 5.19 5.16 5.01 0 0 0 0 0 0 0 0 0 0 1 35
36 5.89 110.28 5.76 5.58 5.39 5.19 5.03 0 0 0 0 0 0 0 0 0 0 0 36
37 5.98 110.17 5.89 5.76 5.58 5.39 5.13 1 0 0 0 0 0 0 0 0 0 0 37
38 6.02 109.99 5.98 5.89 5.76 5.58 5.18 0 1 0 0 0 0 0 0 0 0 0 38
39 5.62 109.26 6.02 5.98 5.89 5.76 5.21 0 0 1 0 0 0 0 0 0 0 0 39
40 4.87 109.11 5.62 6.02 5.98 5.89 5.26 0 0 0 1 0 0 0 0 0 0 0 40
41 4.24 107.06 4.87 5.62 6.02 5.98 5.25 0 0 0 0 1 0 0 0 0 0 0 41
42 4.02 109.53 4.24 4.87 5.62 6.02 5.20 0 0 0 0 0 1 0 0 0 0 0 42
43 3.74 108.92 4.02 4.24 4.87 5.62 5.16 0 0 0 0 0 0 1 0 0 0 0 43
44 3.45 109.24 3.74 4.02 4.24 4.87 5.19 0 0 0 0 0 0 0 1 0 0 0 44
45 3.34 109.12 3.45 3.74 4.02 4.24 5.39 0 0 0 0 0 0 0 0 1 0 0 45
46 3.21 109.00 3.34 3.45 3.74 4.02 5.58 0 0 0 0 0 0 0 0 0 1 0 46
47 3.12 107.23 3.21 3.34 3.45 3.74 5.76 0 0 0 0 0 0 0 0 0 0 1 47
48 3.04 109.49 3.12 3.21 3.34 3.45 5.89 0 0 0 0 0 0 0 0 0 0 0 48
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X y1 y2 y3 y4
-1.397320 0.016172 2.001573 -1.564241 0.790011 -0.262257
y12 M1 M2 M3 M4 M5
-0.042426 0.113985 0.081565 -0.029441 0.037083 0.087844
M6 M7 M8 M9 M10 M11
0.094825 0.046743 0.104716 0.124092 0.028258 0.220714
t
-0.002062
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.25659 -0.04269 0.01443 0.05732 0.16260
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.397320 2.096764 -0.666 0.510411
X 0.016172 0.022280 0.726 0.473752
y1 2.001573 0.184845 10.828 1.05e-11 ***
y2 -1.564241 0.396628 -3.944 0.000466 ***
y3 0.790011 0.402015 1.965 0.059043 .
y4 -0.262257 0.201902 -1.299 0.204201
y12 -0.042426 0.067502 -0.629 0.534581
M1 0.113985 0.097171 1.173 0.250325
M2 0.081565 0.088295 0.924 0.363229
M3 -0.029441 0.090153 -0.327 0.746342
M4 0.037083 0.098302 0.377 0.708748
M5 0.087844 0.105279 0.834 0.410880
M6 0.094825 0.092153 1.029 0.311986
M7 0.046743 0.091547 0.511 0.613502
M8 0.104716 0.093043 1.125 0.269622
M9 0.124092 0.085637 1.449 0.158050
M10 0.028258 0.083982 0.336 0.738931
M11 0.220714 0.093076 2.371 0.024582 *
t -0.002062 0.005175 -0.399 0.693167
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.115 on 29 degrees of freedom
Multiple R-squared: 0.9901, Adjusted R-squared: 0.9839
F-statistic: 160.6 on 18 and 29 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.09726460 0.19452920 0.9027354
[2,] 0.07038161 0.14076321 0.9296184
[3,] 0.02468726 0.04937453 0.9753127
[4,] 0.03348004 0.06696008 0.9665200
[5,] 0.18882535 0.37765071 0.8111746
> postscript(file="/var/www/html/rcomp/tmp/14hxu1258561689.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/23i8k1258561689.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/3acn41258561689.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/4q7kv1258561689.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/591xv1258561689.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 = 48
Frequency = 1
1 2 3 4 5 6
-0.083210928 -0.001671155 0.054586742 0.073179392 -0.031914152 -0.010744522
7 8 9 10 11 12
0.037741225 -0.098003184 -0.053230392 -0.111512846 0.051601458 -0.128491746
13 14 15 16 17 18
-0.068388071 0.036195843 0.039400500 0.076485412 0.057095250 -0.104090165
19 20 21 22 23 24
0.074762888 0.044088504 -0.118288638 0.090916243 0.017725463 -0.023377118
25 26 27 28 29 30
0.149911102 -0.046563804 0.162604309 0.090649439 -0.006381374 -0.013686146
31 32 33 34 35 36
0.038182858 0.037093060 0.120535879 0.061665864 -0.027934224 0.093881144
37 38 39 40 41 42
0.001687897 0.012039117 -0.256591550 -0.240314242 -0.018799724 0.128520833
43 44 45 46 47 48
-0.150686971 0.016821620 0.050983150 -0.041069261 -0.041392696 0.057987720
> postscript(file="/var/www/html/rcomp/tmp/64uep1258561689.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 = 48
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.083210928 NA
1 -0.001671155 -0.083210928
2 0.054586742 -0.001671155
3 0.073179392 0.054586742
4 -0.031914152 0.073179392
5 -0.010744522 -0.031914152
6 0.037741225 -0.010744522
7 -0.098003184 0.037741225
8 -0.053230392 -0.098003184
9 -0.111512846 -0.053230392
10 0.051601458 -0.111512846
11 -0.128491746 0.051601458
12 -0.068388071 -0.128491746
13 0.036195843 -0.068388071
14 0.039400500 0.036195843
15 0.076485412 0.039400500
16 0.057095250 0.076485412
17 -0.104090165 0.057095250
18 0.074762888 -0.104090165
19 0.044088504 0.074762888
20 -0.118288638 0.044088504
21 0.090916243 -0.118288638
22 0.017725463 0.090916243
23 -0.023377118 0.017725463
24 0.149911102 -0.023377118
25 -0.046563804 0.149911102
26 0.162604309 -0.046563804
27 0.090649439 0.162604309
28 -0.006381374 0.090649439
29 -0.013686146 -0.006381374
30 0.038182858 -0.013686146
31 0.037093060 0.038182858
32 0.120535879 0.037093060
33 0.061665864 0.120535879
34 -0.027934224 0.061665864
35 0.093881144 -0.027934224
36 0.001687897 0.093881144
37 0.012039117 0.001687897
38 -0.256591550 0.012039117
39 -0.240314242 -0.256591550
40 -0.018799724 -0.240314242
41 0.128520833 -0.018799724
42 -0.150686971 0.128520833
43 0.016821620 -0.150686971
44 0.050983150 0.016821620
45 -0.041069261 0.050983150
46 -0.041392696 -0.041069261
47 0.057987720 -0.041392696
48 NA 0.057987720
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.001671155 -0.083210928
[2,] 0.054586742 -0.001671155
[3,] 0.073179392 0.054586742
[4,] -0.031914152 0.073179392
[5,] -0.010744522 -0.031914152
[6,] 0.037741225 -0.010744522
[7,] -0.098003184 0.037741225
[8,] -0.053230392 -0.098003184
[9,] -0.111512846 -0.053230392
[10,] 0.051601458 -0.111512846
[11,] -0.128491746 0.051601458
[12,] -0.068388071 -0.128491746
[13,] 0.036195843 -0.068388071
[14,] 0.039400500 0.036195843
[15,] 0.076485412 0.039400500
[16,] 0.057095250 0.076485412
[17,] -0.104090165 0.057095250
[18,] 0.074762888 -0.104090165
[19,] 0.044088504 0.074762888
[20,] -0.118288638 0.044088504
[21,] 0.090916243 -0.118288638
[22,] 0.017725463 0.090916243
[23,] -0.023377118 0.017725463
[24,] 0.149911102 -0.023377118
[25,] -0.046563804 0.149911102
[26,] 0.162604309 -0.046563804
[27,] 0.090649439 0.162604309
[28,] -0.006381374 0.090649439
[29,] -0.013686146 -0.006381374
[30,] 0.038182858 -0.013686146
[31,] 0.037093060 0.038182858
[32,] 0.120535879 0.037093060
[33,] 0.061665864 0.120535879
[34,] -0.027934224 0.061665864
[35,] 0.093881144 -0.027934224
[36,] 0.001687897 0.093881144
[37,] 0.012039117 0.001687897
[38,] -0.256591550 0.012039117
[39,] -0.240314242 -0.256591550
[40,] -0.018799724 -0.240314242
[41,] 0.128520833 -0.018799724
[42,] -0.150686971 0.128520833
[43,] 0.016821620 -0.150686971
[44,] 0.050983150 0.016821620
[45,] -0.041069261 0.050983150
[46,] -0.041392696 -0.041069261
[47,] 0.057987720 -0.041392696
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.001671155 -0.083210928
2 0.054586742 -0.001671155
3 0.073179392 0.054586742
4 -0.031914152 0.073179392
5 -0.010744522 -0.031914152
6 0.037741225 -0.010744522
7 -0.098003184 0.037741225
8 -0.053230392 -0.098003184
9 -0.111512846 -0.053230392
10 0.051601458 -0.111512846
11 -0.128491746 0.051601458
12 -0.068388071 -0.128491746
13 0.036195843 -0.068388071
14 0.039400500 0.036195843
15 0.076485412 0.039400500
16 0.057095250 0.076485412
17 -0.104090165 0.057095250
18 0.074762888 -0.104090165
19 0.044088504 0.074762888
20 -0.118288638 0.044088504
21 0.090916243 -0.118288638
22 0.017725463 0.090916243
23 -0.023377118 0.017725463
24 0.149911102 -0.023377118
25 -0.046563804 0.149911102
26 0.162604309 -0.046563804
27 0.090649439 0.162604309
28 -0.006381374 0.090649439
29 -0.013686146 -0.006381374
30 0.038182858 -0.013686146
31 0.037093060 0.038182858
32 0.120535879 0.037093060
33 0.061665864 0.120535879
34 -0.027934224 0.061665864
35 0.093881144 -0.027934224
36 0.001687897 0.093881144
37 0.012039117 0.001687897
38 -0.256591550 0.012039117
39 -0.240314242 -0.256591550
40 -0.018799724 -0.240314242
41 0.128520833 -0.018799724
42 -0.150686971 0.128520833
43 0.016821620 -0.150686971
44 0.050983150 0.016821620
45 -0.041069261 0.050983150
46 -0.041392696 -0.041069261
47 0.057987720 -0.041392696
> 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/79qa81258561689.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/8swyc1258561689.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/9mb1l1258561689.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/1018n01258561689.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/11bjlv1258561689.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/12s5md1258561689.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/13t4mg1258561689.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/140ffb1258561689.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/15p0py1258561689.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/16tl5r1258561689.tab")
+ }
>
> system("convert tmp/14hxu1258561689.ps tmp/14hxu1258561689.png")
> system("convert tmp/23i8k1258561689.ps tmp/23i8k1258561689.png")
> system("convert tmp/3acn41258561689.ps tmp/3acn41258561689.png")
> system("convert tmp/4q7kv1258561689.ps tmp/4q7kv1258561689.png")
> system("convert tmp/591xv1258561689.ps tmp/591xv1258561689.png")
> system("convert tmp/64uep1258561689.ps tmp/64uep1258561689.png")
> system("convert tmp/79qa81258561689.ps tmp/79qa81258561689.png")
> system("convert tmp/8swyc1258561689.ps tmp/8swyc1258561689.png")
> system("convert tmp/9mb1l1258561689.ps tmp/9mb1l1258561689.png")
> system("convert tmp/1018n01258561689.ps tmp/1018n01258561689.png")
>
>
> proc.time()
user system elapsed
2.239 1.564 2.658