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(1.4
+ ,1.9
+ ,3
+ ,1.5
+ ,-0.7
+ ,-0.7
+ ,-2.9
+ ,-0.8
+ ,1
+ ,1
+ ,1.6
+ ,3.2
+ ,3
+ ,1.5
+ ,-0.7
+ ,-0.7
+ ,-2.9
+ ,-0.8
+ ,-0.8
+ ,0
+ ,3.1
+ ,3.2
+ ,3
+ ,1.5
+ ,-0.7
+ ,-0.7
+ ,-2.9
+ ,-2.9
+ ,-1.3
+ ,3.9
+ ,3.1
+ ,3.2
+ ,3
+ ,1.5
+ ,-0.7
+ ,-0.7
+ ,-0.7
+ ,-0.4
+ ,1
+ ,3.9
+ ,3.1
+ ,3.2
+ ,3
+ ,1.5
+ ,-0.7
+ ,-0.7
+ ,-0.3
+ ,1.3
+ ,1
+ ,3.9
+ ,3.1
+ ,3.2
+ ,3
+ ,1.5
+ ,1.5
+ ,1.4
+ ,0.8
+ ,1.3
+ ,1
+ ,3.9
+ ,3.1
+ ,3.2
+ ,3
+ ,3
+ ,2.6
+ ,1.2
+ ,0.8
+ ,1.3
+ ,1
+ ,3.9
+ ,3.1
+ ,3.2
+ ,3.2
+ ,2.8
+ ,2.9
+ ,1.2
+ ,0.8
+ ,1.3
+ ,1
+ ,3.9
+ ,3.1
+ ,3.1
+ ,2.6
+ ,3.9
+ ,2.9
+ ,1.2
+ ,0.8
+ ,1.3
+ ,1
+ ,3.9
+ ,3.9
+ ,3.4
+ ,4.5
+ ,3.9
+ ,2.9
+ ,1.2
+ ,0.8
+ ,1.3
+ ,1
+ ,1
+ ,1.7
+ ,4.5
+ ,4.5
+ ,3.9
+ ,2.9
+ ,1.2
+ ,0.8
+ ,1.3
+ ,1.3
+ ,1.2
+ ,3.3
+ ,4.5
+ ,4.5
+ ,3.9
+ ,2.9
+ ,1.2
+ ,0.8
+ ,0.8
+ ,0
+ ,2
+ ,3.3
+ ,4.5
+ ,4.5
+ ,3.9
+ ,2.9
+ ,1.2
+ ,1.2
+ ,0
+ ,1.5
+ ,2
+ ,3.3
+ ,4.5
+ ,4.5
+ ,3.9
+ ,2.9
+ ,2.9
+ ,1.6
+ ,1
+ ,1.5
+ ,2
+ ,3.3
+ ,4.5
+ ,4.5
+ ,3.9
+ ,3.9
+ ,2.5
+ ,2.1
+ ,1
+ ,1.5
+ ,2
+ ,3.3
+ ,4.5
+ ,4.5
+ ,4.5
+ ,3.2
+ ,3
+ ,2.1
+ ,1
+ ,1.5
+ ,2
+ ,3.3
+ ,4.5
+ ,4.5
+ ,3.4
+ ,4
+ ,3
+ ,2.1
+ ,1
+ ,1.5
+ ,2
+ ,3.3
+ ,3.3
+ ,2.3
+ ,5.1
+ ,4
+ ,3
+ ,2.1
+ ,1
+ ,1.5
+ ,2
+ ,2
+ ,1.9
+ ,4.5
+ ,5.1
+ ,4
+ ,3
+ ,2.1
+ ,1
+ ,1.5
+ ,1.5
+ ,1.7
+ ,4.2
+ ,4.5
+ ,5.1
+ ,4
+ ,3
+ ,2.1
+ ,1
+ ,1
+ ,1.9
+ ,3.3
+ ,4.2
+ ,4.5
+ ,5.1
+ ,4
+ ,3
+ ,2.1
+ ,2.1
+ ,3.3
+ ,2.7
+ ,3.3
+ ,4.2
+ ,4.5
+ ,5.1
+ ,4
+ ,3
+ ,3
+ ,3.8
+ ,1.8
+ ,2.7
+ ,3.3
+ ,4.2
+ ,4.5
+ ,5.1
+ ,4
+ ,4
+ ,4.4
+ ,1.4
+ ,1.8
+ ,2.7
+ ,3.3
+ ,4.2
+ ,4.5
+ ,5.1
+ ,5.1
+ ,4.5
+ ,0.5
+ ,1.4
+ ,1.8
+ ,2.7
+ ,3.3
+ ,4.2
+ ,4.5
+ ,4.5
+ ,3.5
+ ,-0.4
+ ,0.5
+ ,1.4
+ ,1.8
+ ,2.7
+ ,3.3
+ ,4.2
+ ,4.2
+ ,3
+ ,0.8
+ ,-0.4
+ ,0.5
+ ,1.4
+ ,1.8
+ ,2.7
+ ,3.3
+ ,3.3
+ ,2.8
+ ,0.7
+ ,0.8
+ ,-0.4
+ ,0.5
+ ,1.4
+ ,1.8
+ ,2.7
+ ,2.7
+ ,2.9
+ ,1.9
+ ,0.7
+ ,0.8
+ ,-0.4
+ ,0.5
+ ,1.4
+ ,1.8
+ ,1.8
+ ,2.6
+ ,2
+ ,1.9
+ ,0.7
+ ,0.8
+ ,-0.4
+ ,0.5
+ ,1.4
+ ,1.4
+ ,2.1
+ ,1.1
+ ,2
+ ,1.9
+ ,0.7
+ ,0.8
+ ,-0.4
+ ,0.5
+ ,0.5
+ ,1.5
+ ,0.9
+ ,1.1
+ ,2
+ ,1.9
+ ,0.7
+ ,0.8
+ ,-0.4
+ ,-0.4
+ ,1.1
+ ,0.4
+ ,0.9
+ ,1.1
+ ,2
+ ,1.9
+ ,0.7
+ ,0.8
+ ,0.8
+ ,1.5
+ ,0.7
+ ,0.4
+ ,0.9
+ ,1.1
+ ,2
+ ,1.9
+ ,0.7
+ ,0.7
+ ,1.7
+ ,2.1
+ ,0.7
+ ,0.4
+ ,0.9
+ ,1.1
+ ,2
+ ,1.9
+ ,1.9
+ ,2.3
+ ,2.8
+ ,2.1
+ ,0.7
+ ,0.4
+ ,0.9
+ ,1.1
+ ,2
+ ,2
+ ,2.3
+ ,3.9
+ ,2.8
+ ,2.1
+ ,0.7
+ ,0.4
+ ,0.9
+ ,1.1
+ ,1.1
+ ,1.9
+ ,3.5
+ ,3.9
+ ,2.8
+ ,2.1
+ ,0.7
+ ,0.4
+ ,0.9
+ ,0.9
+ ,2
+ ,2
+ ,3.5
+ ,3.9
+ ,2.8
+ ,2.1
+ ,0.7
+ ,0.4
+ ,0.4
+ ,1.6
+ ,2
+ ,2
+ ,3.5
+ ,3.9
+ ,2.8
+ ,2.1
+ ,0.7
+ ,0.7
+ ,1.2
+ ,1.5
+ ,2
+ ,2
+ ,3.5
+ ,3.9
+ ,2.8
+ ,2.1
+ ,2.1
+ ,1.9
+ ,2.5
+ ,1.5
+ ,2
+ ,2
+ ,3.5
+ ,3.9
+ ,2.8
+ ,2.8
+ ,2.1
+ ,3.1
+ ,2.5
+ ,1.5
+ ,2
+ ,2
+ ,3.5
+ ,3.9
+ ,3.9
+ ,2.4
+ ,2.7
+ ,3.1
+ ,2.5
+ ,1.5
+ ,2
+ ,2
+ ,3.5
+ ,3.5
+ ,2.9
+ ,2.8
+ ,2.7
+ ,3.1
+ ,2.5
+ ,1.5
+ ,2
+ ,2
+ ,2
+ ,2.5
+ ,2.5
+ ,2.8
+ ,2.7
+ ,3.1
+ ,2.5
+ ,1.5
+ ,2
+ ,2
+ ,2.3
+ ,3
+ ,2.5
+ ,2.8
+ ,2.7
+ ,3.1
+ ,2.5
+ ,1.5
+ ,1.5
+ ,2.5
+ ,3.2
+ ,3
+ ,2.5
+ ,2.8
+ ,2.7
+ ,3.1
+ ,2.5
+ ,2.5
+ ,2.6
+ ,2.8
+ ,3.2
+ ,3
+ ,2.5
+ ,2.8
+ ,2.7
+ ,3.1
+ ,3.1
+ ,2.4
+ ,2.4
+ ,2.8
+ ,3.2
+ ,3
+ ,2.5
+ ,2.8
+ ,2.7
+ ,2.7
+ ,2.5
+ ,2
+ ,2.4
+ ,2.8
+ ,3.2
+ ,3
+ ,2.5
+ ,2.8
+ ,2.8
+ ,2.1
+ ,1.8
+ ,2
+ ,2.4
+ ,2.8
+ ,3.2
+ ,3
+ ,2.5
+ ,2.5
+ ,2.2
+ ,1.1
+ ,1.8
+ ,2
+ ,2.4
+ ,2.8
+ ,3.2
+ ,3
+ ,3
+ ,2.7
+ ,-1.5
+ ,1.1
+ ,1.8
+ ,2
+ ,2.4
+ ,2.8
+ ,3.2
+ ,3.2
+ ,3
+ ,-3.7
+ ,-1.5
+ ,1.1
+ ,1.8
+ ,2
+ ,2.4
+ ,2.8)
+ ,dim=c(9
+ ,57)
+ ,dimnames=list(c('bbp'
+ ,'dnst'
+ ,'y1'
+ ,'y2'
+ ,'y3'
+ ,'y4'
+ ,'y5'
+ ,'y6'
+ ,'y7')
+ ,1:57))
> y <- array(NA,dim=c(9,57),dimnames=list(c('bbp','dnst','y1','y2','y3','y4','y5','y6','y7'),1:57))
> 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
bbp dnst y1 y2 y3 y4 y5 y6 y7 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10
1 1.4 1.9 3.0 1.5 -0.7 -0.7 -2.9 -0.8 1.0 1 0 0 0 0 0 0 0 0 0
2 1.0 1.6 3.2 3.0 1.5 -0.7 -0.7 -2.9 -0.8 0 1 0 0 0 0 0 0 0 0
3 -0.8 0.0 3.1 3.2 3.0 1.5 -0.7 -0.7 -2.9 0 0 1 0 0 0 0 0 0 0
4 -2.9 -1.3 3.9 3.1 3.2 3.0 1.5 -0.7 -0.7 0 0 0 1 0 0 0 0 0 0
5 -0.7 -0.4 1.0 3.9 3.1 3.2 3.0 1.5 -0.7 0 0 0 0 1 0 0 0 0 0
6 -0.7 -0.3 1.3 1.0 3.9 3.1 3.2 3.0 1.5 0 0 0 0 0 1 0 0 0 0
7 1.5 1.4 0.8 1.3 1.0 3.9 3.1 3.2 3.0 0 0 0 0 0 0 1 0 0 0
8 3.0 2.6 1.2 0.8 1.3 1.0 3.9 3.1 3.2 0 0 0 0 0 0 0 1 0 0
9 3.2 2.8 2.9 1.2 0.8 1.3 1.0 3.9 3.1 0 0 0 0 0 0 0 0 1 0
10 3.1 2.6 3.9 2.9 1.2 0.8 1.3 1.0 3.9 0 0 0 0 0 0 0 0 0 1
11 3.9 3.4 4.5 3.9 2.9 1.2 0.8 1.3 1.0 0 0 0 0 0 0 0 0 0 0
12 1.0 1.7 4.5 4.5 3.9 2.9 1.2 0.8 1.3 0 0 0 0 0 0 0 0 0 0
13 1.3 1.2 3.3 4.5 4.5 3.9 2.9 1.2 0.8 1 0 0 0 0 0 0 0 0 0
14 0.8 0.0 2.0 3.3 4.5 4.5 3.9 2.9 1.2 0 1 0 0 0 0 0 0 0 0
15 1.2 0.0 1.5 2.0 3.3 4.5 4.5 3.9 2.9 0 0 1 0 0 0 0 0 0 0
16 2.9 1.6 1.0 1.5 2.0 3.3 4.5 4.5 3.9 0 0 0 1 0 0 0 0 0 0
17 3.9 2.5 2.1 1.0 1.5 2.0 3.3 4.5 4.5 0 0 0 0 1 0 0 0 0 0
18 4.5 3.2 3.0 2.1 1.0 1.5 2.0 3.3 4.5 0 0 0 0 0 1 0 0 0 0
19 4.5 3.4 4.0 3.0 2.1 1.0 1.5 2.0 3.3 0 0 0 0 0 0 1 0 0 0
20 3.3 2.3 5.1 4.0 3.0 2.1 1.0 1.5 2.0 0 0 0 0 0 0 0 1 0 0
21 2.0 1.9 4.5 5.1 4.0 3.0 2.1 1.0 1.5 0 0 0 0 0 0 0 0 1 0
22 1.5 1.7 4.2 4.5 5.1 4.0 3.0 2.1 1.0 0 0 0 0 0 0 0 0 0 1
23 1.0 1.9 3.3 4.2 4.5 5.1 4.0 3.0 2.1 0 0 0 0 0 0 0 0 0 0
24 2.1 3.3 2.7 3.3 4.2 4.5 5.1 4.0 3.0 0 0 0 0 0 0 0 0 0 0
25 3.0 3.8 1.8 2.7 3.3 4.2 4.5 5.1 4.0 1 0 0 0 0 0 0 0 0 0
26 4.0 4.4 1.4 1.8 2.7 3.3 4.2 4.5 5.1 0 1 0 0 0 0 0 0 0 0
27 5.1 4.5 0.5 1.4 1.8 2.7 3.3 4.2 4.5 0 0 1 0 0 0 0 0 0 0
28 4.5 3.5 -0.4 0.5 1.4 1.8 2.7 3.3 4.2 0 0 0 1 0 0 0 0 0 0
29 4.2 3.0 0.8 -0.4 0.5 1.4 1.8 2.7 3.3 0 0 0 0 1 0 0 0 0 0
30 3.3 2.8 0.7 0.8 -0.4 0.5 1.4 1.8 2.7 0 0 0 0 0 1 0 0 0 0
31 2.7 2.9 1.9 0.7 0.8 -0.4 0.5 1.4 1.8 0 0 0 0 0 0 1 0 0 0
32 1.8 2.6 2.0 1.9 0.7 0.8 -0.4 0.5 1.4 0 0 0 0 0 0 0 1 0 0
33 1.4 2.1 1.1 2.0 1.9 0.7 0.8 -0.4 0.5 0 0 0 0 0 0 0 0 1 0
34 0.5 1.5 0.9 1.1 2.0 1.9 0.7 0.8 -0.4 0 0 0 0 0 0 0 0 0 1
35 -0.4 1.1 0.4 0.9 1.1 2.0 1.9 0.7 0.8 0 0 0 0 0 0 0 0 0 0
36 0.8 1.5 0.7 0.4 0.9 1.1 2.0 1.9 0.7 0 0 0 0 0 0 0 0 0 0
37 0.7 1.7 2.1 0.7 0.4 0.9 1.1 2.0 1.9 1 0 0 0 0 0 0 0 0 0
38 1.9 2.3 2.8 2.1 0.7 0.4 0.9 1.1 2.0 0 1 0 0 0 0 0 0 0 0
39 2.0 2.3 3.9 2.8 2.1 0.7 0.4 0.9 1.1 0 0 1 0 0 0 0 0 0 0
40 1.1 1.9 3.5 3.9 2.8 2.1 0.7 0.4 0.9 0 0 0 1 0 0 0 0 0 0
41 0.9 2.0 2.0 3.5 3.9 2.8 2.1 0.7 0.4 0 0 0 0 1 0 0 0 0 0
42 0.4 1.6 2.0 2.0 3.5 3.9 2.8 2.1 0.7 0 0 0 0 0 1 0 0 0 0
43 0.7 1.2 1.5 2.0 2.0 3.5 3.9 2.8 2.1 0 0 0 0 0 0 1 0 0 0
44 2.1 1.9 2.5 1.5 2.0 2.0 3.5 3.9 2.8 0 0 0 0 0 0 0 1 0 0
45 2.8 2.1 3.1 2.5 1.5 2.0 2.0 3.5 3.9 0 0 0 0 0 0 0 0 1 0
46 3.9 2.4 2.7 3.1 2.5 1.5 2.0 2.0 3.5 0 0 0 0 0 0 0 0 0 1
47 3.5 2.9 2.8 2.7 3.1 2.5 1.5 2.0 2.0 0 0 0 0 0 0 0 0 0 0
48 2.0 2.5 2.5 2.8 2.7 3.1 2.5 1.5 2.0 0 0 0 0 0 0 0 0 0 0
49 2.0 2.3 3.0 2.5 2.8 2.7 3.1 2.5 1.5 1 0 0 0 0 0 0 0 0 0
50 1.5 2.5 3.2 3.0 2.5 2.8 2.7 3.1 2.5 0 1 0 0 0 0 0 0 0 0
51 2.5 2.6 2.8 3.2 3.0 2.5 2.8 2.7 3.1 0 0 1 0 0 0 0 0 0 0
52 3.1 2.4 2.4 2.8 3.2 3.0 2.5 2.8 2.7 0 0 0 1 0 0 0 0 0 0
53 2.7 2.5 2.0 2.4 2.8 3.2 3.0 2.5 2.8 0 0 0 0 1 0 0 0 0 0
54 2.8 2.1 1.8 2.0 2.4 2.8 3.2 3.0 2.5 0 0 0 0 0 1 0 0 0 0
55 2.5 2.2 1.1 1.8 2.0 2.4 2.8 3.2 3.0 0 0 0 0 0 0 1 0 0 0
56 3.0 2.7 -1.5 1.1 1.8 2.0 2.4 2.8 3.2 0 0 0 0 0 0 0 1 0 0
57 3.2 3.0 -3.7 -1.5 1.1 1.8 2.0 2.4 2.8 0 0 0 0 0 0 0 0 1 0
M11 t
1 0 1
2 0 2
3 0 3
4 0 4
5 0 5
6 0 6
7 0 7
8 0 8
9 0 9
10 0 10
11 1 11
12 0 12
13 0 13
14 0 14
15 0 15
16 0 16
17 0 17
18 0 18
19 0 19
20 0 20
21 0 21
22 0 22
23 1 23
24 0 24
25 0 25
26 0 26
27 0 27
28 0 28
29 0 29
30 0 30
31 0 31
32 0 32
33 0 33
34 0 34
35 1 35
36 0 36
37 0 37
38 0 38
39 0 39
40 0 40
41 0 41
42 0 42
43 0 43
44 0 44
45 0 45
46 0 46
47 1 47
48 0 48
49 0 49
50 0 50
51 0 51
52 0 52
53 0 53
54 0 54
55 0 55
56 0 56
57 0 57
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) dnst y1 y2 y3 y4
-0.568424 0.703035 -0.162054 0.244116 0.266239 -0.441380
y5 y6 y7 M1 M2 M3
-0.171511 0.249940 0.470408 0.051251 0.069341 0.490219
M4 M5 M6 M7 M8 M9
0.432208 0.694046 0.506547 0.454879 0.334984 0.315369
M10 M11 t
0.468923 0.459781 -0.005321
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.00538 -0.42371 -0.05287 0.34148 1.05982
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.568424 0.503963 -1.128 0.26682
dnst 0.703035 0.140534 5.003 1.49e-05 ***
y1 -0.162054 0.122481 -1.323 0.19414
y2 0.244116 0.185397 1.317 0.19625
y3 0.266239 0.182138 1.462 0.15249
y4 -0.441380 0.181964 -2.426 0.02042 *
y5 -0.171511 0.160095 -1.071 0.29116
y6 0.249940 0.164760 1.517 0.13800
y7 0.470408 0.146134 3.219 0.00272 **
M1 0.051251 0.436683 0.117 0.90722
M2 0.069341 0.436292 0.159 0.87461
M3 0.490219 0.447715 1.095 0.28082
M4 0.432208 0.450169 0.960 0.34341
M5 0.694046 0.440497 1.576 0.12387
M6 0.506547 0.446008 1.136 0.26357
M7 0.454879 0.444369 1.024 0.31283
M8 0.334984 0.455761 0.735 0.46710
M9 0.315369 0.463585 0.680 0.50068
M10 0.468923 0.471853 0.994 0.32696
M11 0.459781 0.450034 1.022 0.31376
t -0.005321 0.005760 -0.924 0.36171
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6298 on 36 degrees of freedom
Multiple R-squared: 0.8951, Adjusted R-squared: 0.8369
F-statistic: 15.37 on 20 and 36 DF, p-value: 4.299e-12
> 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.4580515 0.91610300 0.54194850
[2,] 0.9666377 0.06672459 0.03336230
[3,] 0.9608574 0.07828514 0.03914257
[4,] 0.9353073 0.12938541 0.06469271
[5,] 0.9065577 0.18688467 0.09344233
[6,] 0.9568855 0.08622903 0.04311452
[7,] 0.9178635 0.16427303 0.08213651
[8,] 0.9596658 0.08066838 0.04033419
[9,] 0.9485105 0.10297900 0.05148950
[10,] 0.8883368 0.22332647 0.11166324
> postscript(file="/var/www/html/rcomp/tmp/1qjky1258646517.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/29kn61258646517.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/3ud881258646517.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/462i61258646517.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/5f2ll1258646517.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 = 57
Frequency = 1
1 2 3 4 5 6
-0.183720851 0.443863349 0.297800121 -0.719616825 -0.251816022 -1.005378955
7 8 9 10 11 12
0.254636708 -0.068662321 -0.191246701 -0.479049106 0.553493378 -0.296156537
13 14 15 16 17 18
0.823399397 1.059820288 0.653348762 0.528939956 0.011276887 0.178686755
19 20 21 22 23 24
0.327520604 0.856865077 0.177094164 0.030469617 -0.593908840 -0.560184011
25 26 27 28 29 30
-0.798057485 -0.734200582 0.009390412 0.222027794 0.913231880 0.318473397
31 32 33 34 35 36
-0.423710711 -0.449289944 -0.152595473 0.017607260 -0.668847633 0.305984728
37 38 39 40 41 42
-0.525910825 -0.145783966 -0.306683762 -0.493417499 -0.749232197 -0.187989157
43 44 45 46 47 48
-0.052872686 -0.070562671 -0.174736588 0.430972230 0.709263094 0.550355819
49 50 51 52 53 54
0.684289764 -0.623699089 -0.653855532 0.462066574 0.076539451 0.696207960
55 56 57
-0.105573915 -0.268350140 0.341484598
> postscript(file="/var/www/html/rcomp/tmp/6d6bs1258646517.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 = 57
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.183720851 NA
1 0.443863349 -0.183720851
2 0.297800121 0.443863349
3 -0.719616825 0.297800121
4 -0.251816022 -0.719616825
5 -1.005378955 -0.251816022
6 0.254636708 -1.005378955
7 -0.068662321 0.254636708
8 -0.191246701 -0.068662321
9 -0.479049106 -0.191246701
10 0.553493378 -0.479049106
11 -0.296156537 0.553493378
12 0.823399397 -0.296156537
13 1.059820288 0.823399397
14 0.653348762 1.059820288
15 0.528939956 0.653348762
16 0.011276887 0.528939956
17 0.178686755 0.011276887
18 0.327520604 0.178686755
19 0.856865077 0.327520604
20 0.177094164 0.856865077
21 0.030469617 0.177094164
22 -0.593908840 0.030469617
23 -0.560184011 -0.593908840
24 -0.798057485 -0.560184011
25 -0.734200582 -0.798057485
26 0.009390412 -0.734200582
27 0.222027794 0.009390412
28 0.913231880 0.222027794
29 0.318473397 0.913231880
30 -0.423710711 0.318473397
31 -0.449289944 -0.423710711
32 -0.152595473 -0.449289944
33 0.017607260 -0.152595473
34 -0.668847633 0.017607260
35 0.305984728 -0.668847633
36 -0.525910825 0.305984728
37 -0.145783966 -0.525910825
38 -0.306683762 -0.145783966
39 -0.493417499 -0.306683762
40 -0.749232197 -0.493417499
41 -0.187989157 -0.749232197
42 -0.052872686 -0.187989157
43 -0.070562671 -0.052872686
44 -0.174736588 -0.070562671
45 0.430972230 -0.174736588
46 0.709263094 0.430972230
47 0.550355819 0.709263094
48 0.684289764 0.550355819
49 -0.623699089 0.684289764
50 -0.653855532 -0.623699089
51 0.462066574 -0.653855532
52 0.076539451 0.462066574
53 0.696207960 0.076539451
54 -0.105573915 0.696207960
55 -0.268350140 -0.105573915
56 0.341484598 -0.268350140
57 NA 0.341484598
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.443863349 -0.183720851
[2,] 0.297800121 0.443863349
[3,] -0.719616825 0.297800121
[4,] -0.251816022 -0.719616825
[5,] -1.005378955 -0.251816022
[6,] 0.254636708 -1.005378955
[7,] -0.068662321 0.254636708
[8,] -0.191246701 -0.068662321
[9,] -0.479049106 -0.191246701
[10,] 0.553493378 -0.479049106
[11,] -0.296156537 0.553493378
[12,] 0.823399397 -0.296156537
[13,] 1.059820288 0.823399397
[14,] 0.653348762 1.059820288
[15,] 0.528939956 0.653348762
[16,] 0.011276887 0.528939956
[17,] 0.178686755 0.011276887
[18,] 0.327520604 0.178686755
[19,] 0.856865077 0.327520604
[20,] 0.177094164 0.856865077
[21,] 0.030469617 0.177094164
[22,] -0.593908840 0.030469617
[23,] -0.560184011 -0.593908840
[24,] -0.798057485 -0.560184011
[25,] -0.734200582 -0.798057485
[26,] 0.009390412 -0.734200582
[27,] 0.222027794 0.009390412
[28,] 0.913231880 0.222027794
[29,] 0.318473397 0.913231880
[30,] -0.423710711 0.318473397
[31,] -0.449289944 -0.423710711
[32,] -0.152595473 -0.449289944
[33,] 0.017607260 -0.152595473
[34,] -0.668847633 0.017607260
[35,] 0.305984728 -0.668847633
[36,] -0.525910825 0.305984728
[37,] -0.145783966 -0.525910825
[38,] -0.306683762 -0.145783966
[39,] -0.493417499 -0.306683762
[40,] -0.749232197 -0.493417499
[41,] -0.187989157 -0.749232197
[42,] -0.052872686 -0.187989157
[43,] -0.070562671 -0.052872686
[44,] -0.174736588 -0.070562671
[45,] 0.430972230 -0.174736588
[46,] 0.709263094 0.430972230
[47,] 0.550355819 0.709263094
[48,] 0.684289764 0.550355819
[49,] -0.623699089 0.684289764
[50,] -0.653855532 -0.623699089
[51,] 0.462066574 -0.653855532
[52,] 0.076539451 0.462066574
[53,] 0.696207960 0.076539451
[54,] -0.105573915 0.696207960
[55,] -0.268350140 -0.105573915
[56,] 0.341484598 -0.268350140
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.443863349 -0.183720851
2 0.297800121 0.443863349
3 -0.719616825 0.297800121
4 -0.251816022 -0.719616825
5 -1.005378955 -0.251816022
6 0.254636708 -1.005378955
7 -0.068662321 0.254636708
8 -0.191246701 -0.068662321
9 -0.479049106 -0.191246701
10 0.553493378 -0.479049106
11 -0.296156537 0.553493378
12 0.823399397 -0.296156537
13 1.059820288 0.823399397
14 0.653348762 1.059820288
15 0.528939956 0.653348762
16 0.011276887 0.528939956
17 0.178686755 0.011276887
18 0.327520604 0.178686755
19 0.856865077 0.327520604
20 0.177094164 0.856865077
21 0.030469617 0.177094164
22 -0.593908840 0.030469617
23 -0.560184011 -0.593908840
24 -0.798057485 -0.560184011
25 -0.734200582 -0.798057485
26 0.009390412 -0.734200582
27 0.222027794 0.009390412
28 0.913231880 0.222027794
29 0.318473397 0.913231880
30 -0.423710711 0.318473397
31 -0.449289944 -0.423710711
32 -0.152595473 -0.449289944
33 0.017607260 -0.152595473
34 -0.668847633 0.017607260
35 0.305984728 -0.668847633
36 -0.525910825 0.305984728
37 -0.145783966 -0.525910825
38 -0.306683762 -0.145783966
39 -0.493417499 -0.306683762
40 -0.749232197 -0.493417499
41 -0.187989157 -0.749232197
42 -0.052872686 -0.187989157
43 -0.070562671 -0.052872686
44 -0.174736588 -0.070562671
45 0.430972230 -0.174736588
46 0.709263094 0.430972230
47 0.550355819 0.709263094
48 0.684289764 0.550355819
49 -0.623699089 0.684289764
50 -0.653855532 -0.623699089
51 0.462066574 -0.653855532
52 0.076539451 0.462066574
53 0.696207960 0.076539451
54 -0.105573915 0.696207960
55 -0.268350140 -0.105573915
56 0.341484598 -0.268350140
> 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/7umqg1258646517.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/864h21258646517.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/9te6c1258646517.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/10oa8n1258646517.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/1187lh1258646517.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/12mraq1258646517.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/13t1si1258646517.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/14vfy81258646517.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/15nyg11258646517.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/16xvxp1258646517.tab")
+ }
>
> system("convert tmp/1qjky1258646517.ps tmp/1qjky1258646517.png")
> system("convert tmp/29kn61258646517.ps tmp/29kn61258646517.png")
> system("convert tmp/3ud881258646517.ps tmp/3ud881258646517.png")
> system("convert tmp/462i61258646517.ps tmp/462i61258646517.png")
> system("convert tmp/5f2ll1258646517.ps tmp/5f2ll1258646517.png")
> system("convert tmp/6d6bs1258646517.ps tmp/6d6bs1258646517.png")
> system("convert tmp/7umqg1258646517.ps tmp/7umqg1258646517.png")
> system("convert tmp/864h21258646517.ps tmp/864h21258646517.png")
> system("convert tmp/9te6c1258646517.ps tmp/9te6c1258646517.png")
> system("convert tmp/10oa8n1258646517.ps tmp/10oa8n1258646517.png")
>
>
> proc.time()
user system elapsed
2.296 1.525 2.906