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(6.5
+ ,0
+ ,6.3
+ ,6.1
+ ,6.2
+ ,6.3
+ ,6.6
+ ,0
+ ,6.5
+ ,6.3
+ ,6.1
+ ,6.2
+ ,6.5
+ ,0
+ ,6.6
+ ,6.5
+ ,6.3
+ ,6.1
+ ,6.2
+ ,0
+ ,6.5
+ ,6.6
+ ,6.5
+ ,6.3
+ ,6.2
+ ,0
+ ,6.2
+ ,6.5
+ ,6.6
+ ,6.5
+ ,5.9
+ ,0
+ ,6.2
+ ,6.2
+ ,6.5
+ ,6.6
+ ,6.1
+ ,0
+ ,5.9
+ ,6.2
+ ,6.2
+ ,6.5
+ ,6.1
+ ,0
+ ,6.1
+ ,5.9
+ ,6.2
+ ,6.2
+ ,6.1
+ ,0
+ ,6.1
+ ,6.1
+ ,5.9
+ ,6.2
+ ,6.1
+ ,0
+ ,6.1
+ ,6.1
+ ,6.1
+ ,5.9
+ ,6.1
+ ,0
+ ,6.1
+ ,6.1
+ ,6.1
+ ,6.1
+ ,6.4
+ ,0
+ ,6.1
+ ,6.1
+ ,6.1
+ ,6.1
+ ,6.7
+ ,0
+ ,6.4
+ ,6.1
+ ,6.1
+ ,6.1
+ ,6.9
+ ,0
+ ,6.7
+ ,6.4
+ ,6.1
+ ,6.1
+ ,7
+ ,0
+ ,6.9
+ ,6.7
+ ,6.4
+ ,6.1
+ ,7
+ ,0
+ ,7
+ ,6.9
+ ,6.7
+ ,6.4
+ ,6.8
+ ,0
+ ,7
+ ,7
+ ,6.9
+ ,6.7
+ ,6.4
+ ,0
+ ,6.8
+ ,7
+ ,7
+ ,6.9
+ ,5.9
+ ,0
+ ,6.4
+ ,6.8
+ ,7
+ ,7
+ ,5.5
+ ,0
+ ,5.9
+ ,6.4
+ ,6.8
+ ,7
+ ,5.5
+ ,0
+ ,5.5
+ ,5.9
+ ,6.4
+ ,6.8
+ ,5.6
+ ,0
+ ,5.5
+ ,5.5
+ ,5.9
+ ,6.4
+ ,5.8
+ ,0
+ ,5.6
+ ,5.5
+ ,5.5
+ ,5.9
+ ,5.9
+ ,0
+ ,5.8
+ ,5.6
+ ,5.5
+ ,5.5
+ ,6.1
+ ,0
+ ,5.9
+ ,5.8
+ ,5.6
+ ,5.5
+ ,6.1
+ ,0
+ ,6.1
+ ,5.9
+ ,5.8
+ ,5.6
+ ,6
+ ,0
+ ,6.1
+ ,6.1
+ ,5.9
+ ,5.8
+ ,6
+ ,0
+ ,6
+ ,6.1
+ ,6.1
+ ,5.9
+ ,5.9
+ ,0
+ ,6
+ ,6
+ ,6.1
+ ,6.1
+ ,5.5
+ ,0
+ ,5.9
+ ,6
+ ,6
+ ,6.1
+ ,5.6
+ ,0
+ ,5.5
+ ,5.9
+ ,6
+ ,6
+ ,5.4
+ ,0
+ ,5.6
+ ,5.5
+ ,5.9
+ ,6
+ ,5.2
+ ,0
+ ,5.4
+ ,5.6
+ ,5.5
+ ,5.9
+ ,5.2
+ ,0
+ ,5.2
+ ,5.4
+ ,5.6
+ ,5.5
+ ,5.2
+ ,0
+ ,5.2
+ ,5.2
+ ,5.4
+ ,5.6
+ ,5.5
+ ,0
+ ,5.2
+ ,5.2
+ ,5.2
+ ,5.4
+ ,5.8
+ ,1
+ ,5.5
+ ,5.2
+ ,5.2
+ ,5.2
+ ,5.8
+ ,1
+ ,5.8
+ ,5.5
+ ,5.2
+ ,5.2
+ ,5.5
+ ,1
+ ,5.8
+ ,5.8
+ ,5.5
+ ,5.2
+ ,5.3
+ ,1
+ ,5.5
+ ,5.8
+ ,5.8
+ ,5.5
+ ,5.1
+ ,1
+ ,5.3
+ ,5.5
+ ,5.8
+ ,5.8
+ ,5.2
+ ,1
+ ,5.1
+ ,5.3
+ ,5.5
+ ,5.8
+ ,5.8
+ ,1
+ ,5.2
+ ,5.1
+ ,5.3
+ ,5.5
+ ,5.8
+ ,1
+ ,5.8
+ ,5.2
+ ,5.1
+ ,5.3
+ ,5.5
+ ,1
+ ,5.8
+ ,5.8
+ ,5.2
+ ,5.1
+ ,5
+ ,1
+ ,5.5
+ ,5.8
+ ,5.8
+ ,5.2
+ ,4.9
+ ,1
+ ,5
+ ,5.5
+ ,5.8
+ ,5.8
+ ,5.3
+ ,1
+ ,4.9
+ ,5
+ ,5.5
+ ,5.8
+ ,6.1
+ ,1
+ ,5.3
+ ,4.9
+ ,5
+ ,5.5
+ ,6.5
+ ,1
+ ,6.1
+ ,5.3
+ ,4.9
+ ,5
+ ,6.8
+ ,1
+ ,6.5
+ ,6.1
+ ,5.3
+ ,4.9
+ ,6.6
+ ,1
+ ,6.8
+ ,6.5
+ ,6.1
+ ,5.3
+ ,6.4
+ ,1
+ ,6.6
+ ,6.8
+ ,6.5
+ ,6.1
+ ,6.4
+ ,1
+ ,6.4
+ ,6.6
+ ,6.8
+ ,6.5)
+ ,dim=c(6
+ ,54)
+ ,dimnames=list(c('y'
+ ,'x'
+ ,'y-1'
+ ,'y-2'
+ ,'y-3'
+ ,'y-4')
+ ,1:54))
> y <- array(NA,dim=c(6,54),dimnames=list(c('y','x','y-1','y-2','y-3','y-4'),1:54))
> 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 = '2'
> #'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 y-1 y-2 y-3 y-4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 0 6.5 6.3 6.1 6.2 6.3 1 0 0 0 0 0 0 0 0 0 0 1
2 0 6.6 6.5 6.3 6.1 6.2 0 1 0 0 0 0 0 0 0 0 0 2
3 0 6.5 6.6 6.5 6.3 6.1 0 0 1 0 0 0 0 0 0 0 0 3
4 0 6.2 6.5 6.6 6.5 6.3 0 0 0 1 0 0 0 0 0 0 0 4
5 0 6.2 6.2 6.5 6.6 6.5 0 0 0 0 1 0 0 0 0 0 0 5
6 0 5.9 6.2 6.2 6.5 6.6 0 0 0 0 0 1 0 0 0 0 0 6
7 0 6.1 5.9 6.2 6.2 6.5 0 0 0 0 0 0 1 0 0 0 0 7
8 0 6.1 6.1 5.9 6.2 6.2 0 0 0 0 0 0 0 1 0 0 0 8
9 0 6.1 6.1 6.1 5.9 6.2 0 0 0 0 0 0 0 0 1 0 0 9
10 0 6.1 6.1 6.1 6.1 5.9 0 0 0 0 0 0 0 0 0 1 0 10
11 0 6.1 6.1 6.1 6.1 6.1 0 0 0 0 0 0 0 0 0 0 1 11
12 0 6.4 6.1 6.1 6.1 6.1 0 0 0 0 0 0 0 0 0 0 0 12
13 0 6.7 6.4 6.1 6.1 6.1 1 0 0 0 0 0 0 0 0 0 0 13
14 0 6.9 6.7 6.4 6.1 6.1 0 1 0 0 0 0 0 0 0 0 0 14
15 0 7.0 6.9 6.7 6.4 6.1 0 0 1 0 0 0 0 0 0 0 0 15
16 0 7.0 7.0 6.9 6.7 6.4 0 0 0 1 0 0 0 0 0 0 0 16
17 0 6.8 7.0 7.0 6.9 6.7 0 0 0 0 1 0 0 0 0 0 0 17
18 0 6.4 6.8 7.0 7.0 6.9 0 0 0 0 0 1 0 0 0 0 0 18
19 0 5.9 6.4 6.8 7.0 7.0 0 0 0 0 0 0 1 0 0 0 0 19
20 0 5.5 5.9 6.4 6.8 7.0 0 0 0 0 0 0 0 1 0 0 0 20
21 0 5.5 5.5 5.9 6.4 6.8 0 0 0 0 0 0 0 0 1 0 0 21
22 0 5.6 5.5 5.5 5.9 6.4 0 0 0 0 0 0 0 0 0 1 0 22
23 0 5.8 5.6 5.5 5.5 5.9 0 0 0 0 0 0 0 0 0 0 1 23
24 0 5.9 5.8 5.6 5.5 5.5 0 0 0 0 0 0 0 0 0 0 0 24
25 0 6.1 5.9 5.8 5.6 5.5 1 0 0 0 0 0 0 0 0 0 0 25
26 0 6.1 6.1 5.9 5.8 5.6 0 1 0 0 0 0 0 0 0 0 0 26
27 0 6.0 6.1 6.1 5.9 5.8 0 0 1 0 0 0 0 0 0 0 0 27
28 0 6.0 6.0 6.1 6.1 5.9 0 0 0 1 0 0 0 0 0 0 0 28
29 0 5.9 6.0 6.0 6.1 6.1 0 0 0 0 1 0 0 0 0 0 0 29
30 0 5.5 5.9 6.0 6.0 6.1 0 0 0 0 0 1 0 0 0 0 0 30
31 0 5.6 5.5 5.9 6.0 6.0 0 0 0 0 0 0 1 0 0 0 0 31
32 0 5.4 5.6 5.5 5.9 6.0 0 0 0 0 0 0 0 1 0 0 0 32
33 0 5.2 5.4 5.6 5.5 5.9 0 0 0 0 0 0 0 0 1 0 0 33
34 0 5.2 5.2 5.4 5.6 5.5 0 0 0 0 0 0 0 0 0 1 0 34
35 0 5.2 5.2 5.2 5.4 5.6 0 0 0 0 0 0 0 0 0 0 1 35
36 0 5.5 5.2 5.2 5.2 5.4 0 0 0 0 0 0 0 0 0 0 0 36
37 1 5.8 5.5 5.2 5.2 5.2 1 0 0 0 0 0 0 0 0 0 0 37
38 1 5.8 5.8 5.5 5.2 5.2 0 1 0 0 0 0 0 0 0 0 0 38
39 1 5.5 5.8 5.8 5.5 5.2 0 0 1 0 0 0 0 0 0 0 0 39
40 1 5.3 5.5 5.8 5.8 5.5 0 0 0 1 0 0 0 0 0 0 0 40
41 1 5.1 5.3 5.5 5.8 5.8 0 0 0 0 1 0 0 0 0 0 0 41
42 1 5.2 5.1 5.3 5.5 5.8 0 0 0 0 0 1 0 0 0 0 0 42
43 1 5.8 5.2 5.1 5.3 5.5 0 0 0 0 0 0 1 0 0 0 0 43
44 1 5.8 5.8 5.2 5.1 5.3 0 0 0 0 0 0 0 1 0 0 0 44
45 1 5.5 5.8 5.8 5.2 5.1 0 0 0 0 0 0 0 0 1 0 0 45
46 1 5.0 5.5 5.8 5.8 5.2 0 0 0 0 0 0 0 0 0 1 0 46
47 1 4.9 5.0 5.5 5.8 5.8 0 0 0 0 0 0 0 0 0 0 1 47
48 1 5.3 4.9 5.0 5.5 5.8 0 0 0 0 0 0 0 0 0 0 0 48
49 1 6.1 5.3 4.9 5.0 5.5 1 0 0 0 0 0 0 0 0 0 0 49
50 1 6.5 6.1 5.3 4.9 5.0 0 1 0 0 0 0 0 0 0 0 0 50
51 1 6.8 6.5 6.1 5.3 4.9 0 0 1 0 0 0 0 0 0 0 0 51
52 1 6.6 6.8 6.5 6.1 5.3 0 0 0 1 0 0 0 0 0 0 0 52
53 1 6.4 6.6 6.8 6.5 6.1 0 0 0 0 1 0 0 0 0 0 0 53
54 1 6.4 6.4 6.6 6.8 6.5 0 0 0 0 0 1 0 0 0 0 0 54
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) y `y-1` `y-2` `y-3` `y-4`
0.82827 0.18169 -0.23108 -0.11398 0.40312 -0.45191
M1 M2 M3 M4 M5 M6
0.26008 0.28241 0.23476 0.22404 0.30123 0.34072
M7 M8 M9 M10 M11 t
0.23016 0.22606 0.24993 0.04896 0.09285 0.01987
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.46080 -0.09651 0.07203 0.20782 0.36134
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.828266 1.081117 0.766 0.449
y 0.181691 0.290964 0.624 0.536
`y-1` -0.231080 0.523453 -0.441 0.662
`y-2` -0.113985 0.549486 -0.207 0.837
`y-3` 0.403116 0.528222 0.763 0.450
`y-4` -0.451911 0.342055 -1.321 0.195
M1 0.260083 0.213724 1.217 0.232
M2 0.282414 0.234449 1.205 0.236
M3 0.234756 0.238767 0.983 0.332
M4 0.224038 0.241102 0.929 0.359
M5 0.301229 0.240233 1.254 0.218
M6 0.340721 0.255050 1.336 0.190
M7 0.230156 0.233398 0.986 0.331
M8 0.226055 0.269844 0.838 0.408
M9 0.249926 0.268449 0.931 0.358
M10 0.048957 0.234486 0.209 0.836
M11 0.092850 0.226469 0.410 0.684
t 0.019867 0.004439 4.476 7.37e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3059 on 36 degrees of freedom
Multiple R-squared: 0.7193, Adjusted R-squared: 0.5868
F-statistic: 5.427 on 17 and 36 DF, p-value: 1.017e-05
> 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 0 1
[2,] 0 0 1
[3,] 0 0 1
[4,] 0 0 1
[5,] 0 0 1
[6,] 0 0 1
[7,] 0 0 1
[8,] 0 0 1
[9,] 0 0 1
[10,] 0 0 1
[11,] 0 0 1
[12,] 0 0 1
[13,] 0 0 1
> postscript(file="/var/www/html/rcomp/tmp/1ow4q1259009398.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/2abg41259009398.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/3t1vb1259009398.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/4xfm51259009398.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/5v5ae1259009398.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 = 54
Frequency = 1
1 2 3 4 5 6
0.20962737 0.21339473 0.17944504 0.22285365 0.09514421 0.14159972
7 8 9 10 11 12
0.20237996 0.06306110 0.16305496 0.12796142 0.15458351 0.17305993
13 14 15 16 17 18
-0.09207314 -0.06708884 -0.09799065 -0.04659525 -0.04096611 -0.02379415
19 20 21 22 23 24
0.08771208 0.06411173 -0.05818637 0.07994721 -0.06175244 -0.13008801
25 26 27 28 29 30
-0.44078265 -0.46079723 -0.34196986 -0.40965807 -0.40956288 -0.37904201
31 32 33 34 35 36
-0.45553398 -0.41713613 -0.34329814 -0.45228415 -0.41302695 -0.40430960
37 38 39 40 41 42
0.24017504 0.30149747 0.29705594 0.26956056 0.26400304 0.23839686
43 44 45 46 47 48
0.16544195 0.28996330 0.23842955 0.24437552 0.32019588 0.36133769
49 50 51 52 53 54
0.08305339 0.01299387 -0.03654047 -0.03616089 0.09138174 0.02283957
> postscript(file="/var/www/html/rcomp/tmp/6bvdd1259009398.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 = 54
Frequency = 1
lag(myerror, k = 1) myerror
0 0.20962737 NA
1 0.21339473 0.20962737
2 0.17944504 0.21339473
3 0.22285365 0.17944504
4 0.09514421 0.22285365
5 0.14159972 0.09514421
6 0.20237996 0.14159972
7 0.06306110 0.20237996
8 0.16305496 0.06306110
9 0.12796142 0.16305496
10 0.15458351 0.12796142
11 0.17305993 0.15458351
12 -0.09207314 0.17305993
13 -0.06708884 -0.09207314
14 -0.09799065 -0.06708884
15 -0.04659525 -0.09799065
16 -0.04096611 -0.04659525
17 -0.02379415 -0.04096611
18 0.08771208 -0.02379415
19 0.06411173 0.08771208
20 -0.05818637 0.06411173
21 0.07994721 -0.05818637
22 -0.06175244 0.07994721
23 -0.13008801 -0.06175244
24 -0.44078265 -0.13008801
25 -0.46079723 -0.44078265
26 -0.34196986 -0.46079723
27 -0.40965807 -0.34196986
28 -0.40956288 -0.40965807
29 -0.37904201 -0.40956288
30 -0.45553398 -0.37904201
31 -0.41713613 -0.45553398
32 -0.34329814 -0.41713613
33 -0.45228415 -0.34329814
34 -0.41302695 -0.45228415
35 -0.40430960 -0.41302695
36 0.24017504 -0.40430960
37 0.30149747 0.24017504
38 0.29705594 0.30149747
39 0.26956056 0.29705594
40 0.26400304 0.26956056
41 0.23839686 0.26400304
42 0.16544195 0.23839686
43 0.28996330 0.16544195
44 0.23842955 0.28996330
45 0.24437552 0.23842955
46 0.32019588 0.24437552
47 0.36133769 0.32019588
48 0.08305339 0.36133769
49 0.01299387 0.08305339
50 -0.03654047 0.01299387
51 -0.03616089 -0.03654047
52 0.09138174 -0.03616089
53 0.02283957 0.09138174
54 NA 0.02283957
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.21339473 0.20962737
[2,] 0.17944504 0.21339473
[3,] 0.22285365 0.17944504
[4,] 0.09514421 0.22285365
[5,] 0.14159972 0.09514421
[6,] 0.20237996 0.14159972
[7,] 0.06306110 0.20237996
[8,] 0.16305496 0.06306110
[9,] 0.12796142 0.16305496
[10,] 0.15458351 0.12796142
[11,] 0.17305993 0.15458351
[12,] -0.09207314 0.17305993
[13,] -0.06708884 -0.09207314
[14,] -0.09799065 -0.06708884
[15,] -0.04659525 -0.09799065
[16,] -0.04096611 -0.04659525
[17,] -0.02379415 -0.04096611
[18,] 0.08771208 -0.02379415
[19,] 0.06411173 0.08771208
[20,] -0.05818637 0.06411173
[21,] 0.07994721 -0.05818637
[22,] -0.06175244 0.07994721
[23,] -0.13008801 -0.06175244
[24,] -0.44078265 -0.13008801
[25,] -0.46079723 -0.44078265
[26,] -0.34196986 -0.46079723
[27,] -0.40965807 -0.34196986
[28,] -0.40956288 -0.40965807
[29,] -0.37904201 -0.40956288
[30,] -0.45553398 -0.37904201
[31,] -0.41713613 -0.45553398
[32,] -0.34329814 -0.41713613
[33,] -0.45228415 -0.34329814
[34,] -0.41302695 -0.45228415
[35,] -0.40430960 -0.41302695
[36,] 0.24017504 -0.40430960
[37,] 0.30149747 0.24017504
[38,] 0.29705594 0.30149747
[39,] 0.26956056 0.29705594
[40,] 0.26400304 0.26956056
[41,] 0.23839686 0.26400304
[42,] 0.16544195 0.23839686
[43,] 0.28996330 0.16544195
[44,] 0.23842955 0.28996330
[45,] 0.24437552 0.23842955
[46,] 0.32019588 0.24437552
[47,] 0.36133769 0.32019588
[48,] 0.08305339 0.36133769
[49,] 0.01299387 0.08305339
[50,] -0.03654047 0.01299387
[51,] -0.03616089 -0.03654047
[52,] 0.09138174 -0.03616089
[53,] 0.02283957 0.09138174
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.21339473 0.20962737
2 0.17944504 0.21339473
3 0.22285365 0.17944504
4 0.09514421 0.22285365
5 0.14159972 0.09514421
6 0.20237996 0.14159972
7 0.06306110 0.20237996
8 0.16305496 0.06306110
9 0.12796142 0.16305496
10 0.15458351 0.12796142
11 0.17305993 0.15458351
12 -0.09207314 0.17305993
13 -0.06708884 -0.09207314
14 -0.09799065 -0.06708884
15 -0.04659525 -0.09799065
16 -0.04096611 -0.04659525
17 -0.02379415 -0.04096611
18 0.08771208 -0.02379415
19 0.06411173 0.08771208
20 -0.05818637 0.06411173
21 0.07994721 -0.05818637
22 -0.06175244 0.07994721
23 -0.13008801 -0.06175244
24 -0.44078265 -0.13008801
25 -0.46079723 -0.44078265
26 -0.34196986 -0.46079723
27 -0.40965807 -0.34196986
28 -0.40956288 -0.40965807
29 -0.37904201 -0.40956288
30 -0.45553398 -0.37904201
31 -0.41713613 -0.45553398
32 -0.34329814 -0.41713613
33 -0.45228415 -0.34329814
34 -0.41302695 -0.45228415
35 -0.40430960 -0.41302695
36 0.24017504 -0.40430960
37 0.30149747 0.24017504
38 0.29705594 0.30149747
39 0.26956056 0.29705594
40 0.26400304 0.26956056
41 0.23839686 0.26400304
42 0.16544195 0.23839686
43 0.28996330 0.16544195
44 0.23842955 0.28996330
45 0.24437552 0.23842955
46 0.32019588 0.24437552
47 0.36133769 0.32019588
48 0.08305339 0.36133769
49 0.01299387 0.08305339
50 -0.03654047 0.01299387
51 -0.03616089 -0.03654047
52 0.09138174 -0.03616089
53 0.02283957 0.09138174
> 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/7y1ji1259009398.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/8a0vw1259009398.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/9yw4g1259009398.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/10v7k71259009398.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/11772q1259009398.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/12gig01259009398.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/13ooki1259009398.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/14q3dc1259009398.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/157zkn1259009398.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/16ggfm1259009398.tab")
+ }
>
> system("convert tmp/1ow4q1259009398.ps tmp/1ow4q1259009398.png")
> system("convert tmp/2abg41259009398.ps tmp/2abg41259009398.png")
> system("convert tmp/3t1vb1259009398.ps tmp/3t1vb1259009398.png")
> system("convert tmp/4xfm51259009398.ps tmp/4xfm51259009398.png")
> system("convert tmp/5v5ae1259009398.ps tmp/5v5ae1259009398.png")
> system("convert tmp/6bvdd1259009398.ps tmp/6bvdd1259009398.png")
> system("convert tmp/7y1ji1259009398.ps tmp/7y1ji1259009398.png")
> system("convert tmp/8a0vw1259009398.ps tmp/8a0vw1259009398.png")
> system("convert tmp/9yw4g1259009398.ps tmp/9yw4g1259009398.png")
> system("convert tmp/10v7k71259009398.ps tmp/10v7k71259009398.png")
>
>
> proc.time()
user system elapsed
2.351 1.591 2.930