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(105.7
+ ,102.9
+ ,105.7
+ ,105.6
+ ,105.4
+ ,105.4
+ ,105.8
+ ,103.1
+ ,105.8
+ ,105.7
+ ,105.6
+ ,105.4
+ ,105.8
+ ,103
+ ,105.8
+ ,105.8
+ ,105.7
+ ,105.6
+ ,105.8
+ ,102.8
+ ,105.8
+ ,105.8
+ ,105.8
+ ,105.7
+ ,105.9
+ ,102.5
+ ,105.9
+ ,105.8
+ ,105.8
+ ,105.8
+ ,106.1
+ ,101.9
+ ,106.1
+ ,105.9
+ ,105.8
+ ,105.8
+ ,106.4
+ ,101.9
+ ,106.4
+ ,106.1
+ ,105.9
+ ,105.8
+ ,106.4
+ ,101.8
+ ,106.4
+ ,106.4
+ ,106.1
+ ,105.9
+ ,106.3
+ ,102
+ ,106.3
+ ,106.4
+ ,106.4
+ ,106.1
+ ,106.2
+ ,102.6
+ ,106.2
+ ,106.3
+ ,106.4
+ ,106.4
+ ,106.2
+ ,102.5
+ ,106.2
+ ,106.2
+ ,106.3
+ ,106.4
+ ,106.3
+ ,102.5
+ ,106.3
+ ,106.2
+ ,106.2
+ ,106.3
+ ,106.4
+ ,101.6
+ ,106.4
+ ,106.3
+ ,106.2
+ ,106.2
+ ,106.5
+ ,101.4
+ ,106.5
+ ,106.4
+ ,106.3
+ ,106.2
+ ,106.6
+ ,100.8
+ ,106.6
+ ,106.5
+ ,106.4
+ ,106.3
+ ,106.6
+ ,101.1
+ ,106.6
+ ,106.6
+ ,106.5
+ ,106.4
+ ,106.6
+ ,101.3
+ ,106.6
+ ,106.6
+ ,106.6
+ ,106.5
+ ,106.8
+ ,101.2
+ ,106.8
+ ,106.6
+ ,106.6
+ ,106.6
+ ,107
+ ,101.3
+ ,107
+ ,106.8
+ ,106.6
+ ,106.6
+ ,107.2
+ ,101.1
+ ,107.2
+ ,107
+ ,106.8
+ ,106.6
+ ,107.3
+ ,101.3
+ ,107.3
+ ,107.2
+ ,107
+ ,106.8
+ ,107.5
+ ,101.2
+ ,107.5
+ ,107.3
+ ,107.2
+ ,107
+ ,107.6
+ ,101.6
+ ,107.6
+ ,107.5
+ ,107.3
+ ,107.2
+ ,107.6
+ ,101.7
+ ,107.6
+ ,107.6
+ ,107.5
+ ,107.3
+ ,107.7
+ ,101.5
+ ,107.7
+ ,107.6
+ ,107.6
+ ,107.5
+ ,107.7
+ ,100.9
+ ,107.7
+ ,107.7
+ ,107.6
+ ,107.6
+ ,107.7
+ ,101.5
+ ,107.7
+ ,107.7
+ ,107.7
+ ,107.6
+ ,107.7
+ ,101.4
+ ,107.7
+ ,107.7
+ ,107.7
+ ,107.7
+ ,107.6
+ ,101.6
+ ,107.6
+ ,107.7
+ ,107.7
+ ,107.7
+ ,107.7
+ ,101.7
+ ,107.7
+ ,107.6
+ ,107.7
+ ,107.7
+ ,107.9
+ ,101.4
+ ,107.9
+ ,107.7
+ ,107.6
+ ,107.7
+ ,107.9
+ ,101.8
+ ,107.9
+ ,107.9
+ ,107.7
+ ,107.6
+ ,107.9
+ ,101.7
+ ,107.9
+ ,107.9
+ ,107.9
+ ,107.7
+ ,107.8
+ ,101.4
+ ,107.8
+ ,107.9
+ ,107.9
+ ,107.9
+ ,107.6
+ ,101.2
+ ,107.6
+ ,107.8
+ ,107.9
+ ,107.9
+ ,107.4
+ ,101
+ ,107.4
+ ,107.6
+ ,107.8
+ ,107.9
+ ,107
+ ,101.7
+ ,107
+ ,107.4
+ ,107.6
+ ,107.8
+ ,107
+ ,102.4
+ ,107
+ ,107
+ ,107.4
+ ,107.6
+ ,107.2
+ ,102
+ ,107.2
+ ,107
+ ,107
+ ,107.4
+ ,107.5
+ ,102.1
+ ,107.5
+ ,107.2
+ ,107
+ ,107
+ ,107.8
+ ,102
+ ,107.8
+ ,107.5
+ ,107.2
+ ,107
+ ,107.8
+ ,101.8
+ ,107.8
+ ,107.8
+ ,107.5
+ ,107.2
+ ,107.7
+ ,102.7
+ ,107.7
+ ,107.8
+ ,107.8
+ ,107.5
+ ,107.6
+ ,102.3
+ ,107.6
+ ,107.7
+ ,107.8
+ ,107.8
+ ,107.6
+ ,101.9
+ ,107.6
+ ,107.6
+ ,107.7
+ ,107.8
+ ,107.5
+ ,102
+ ,107.5
+ ,107.6
+ ,107.6
+ ,107.7
+ ,107.5
+ ,102.3
+ ,107.5
+ ,107.5
+ ,107.6
+ ,107.6
+ ,107.6
+ ,102.8
+ ,107.6
+ ,107.5
+ ,107.5
+ ,107.6
+ ,107.6
+ ,102.4
+ ,107.6
+ ,107.6
+ ,107.5
+ ,107.5
+ ,107.9
+ ,102.3
+ ,107.9
+ ,107.6
+ ,107.6
+ ,107.5
+ ,107.6
+ ,102.7
+ ,107.6
+ ,107.9
+ ,107.6
+ ,107.6
+ ,107.5
+ ,102.7
+ ,107.5
+ ,107.6
+ ,107.9
+ ,107.6
+ ,107.5
+ ,102.9
+ ,107.5
+ ,107.5
+ ,107.6
+ ,107.9
+ ,107.6
+ ,103
+ ,107.6
+ ,107.5
+ ,107.5
+ ,107.6
+ ,107.7
+ ,102.2
+ ,107.7
+ ,107.6
+ ,107.5
+ ,107.5
+ ,107.8
+ ,102.3
+ ,107.8
+ ,107.7
+ ,107.6
+ ,107.5
+ ,107.9
+ ,102.8
+ ,107.9
+ ,107.8
+ ,107.7
+ ,107.6
+ ,107.9
+ ,102.8
+ ,107.9
+ ,107.9
+ ,107.8
+ ,107.7)
+ ,dim=c(6
+ ,58)
+ ,dimnames=list(c('Werkl'
+ ,'Infl'
+ ,'Yt-1'
+ ,'Yt-2'
+ ,'Yt-3'
+ ,'Yt-4')
+ ,1:58))
> y <- array(NA,dim=c(6,58),dimnames=list(c('Werkl','Infl','Yt-1','Yt-2','Yt-3','Yt-4'),1:58))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = '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
Werkl Infl Yt-1 Yt-2 Yt-3 Yt-4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 105.7 102.9 105.7 105.6 105.4 105.4 1 0 0 0 0 0 0 0 0 0 0 1
2 105.8 103.1 105.8 105.7 105.6 105.4 0 1 0 0 0 0 0 0 0 0 0 2
3 105.8 103.0 105.8 105.8 105.7 105.6 0 0 1 0 0 0 0 0 0 0 0 3
4 105.8 102.8 105.8 105.8 105.8 105.7 0 0 0 1 0 0 0 0 0 0 0 4
5 105.9 102.5 105.9 105.8 105.8 105.8 0 0 0 0 1 0 0 0 0 0 0 5
6 106.1 101.9 106.1 105.9 105.8 105.8 0 0 0 0 0 1 0 0 0 0 0 6
7 106.4 101.9 106.4 106.1 105.9 105.8 0 0 0 0 0 0 1 0 0 0 0 7
8 106.4 101.8 106.4 106.4 106.1 105.9 0 0 0 0 0 0 0 1 0 0 0 8
9 106.3 102.0 106.3 106.4 106.4 106.1 0 0 0 0 0 0 0 0 1 0 0 9
10 106.2 102.6 106.2 106.3 106.4 106.4 0 0 0 0 0 0 0 0 0 1 0 10
11 106.2 102.5 106.2 106.2 106.3 106.4 0 0 0 0 0 0 0 0 0 0 1 11
12 106.3 102.5 106.3 106.2 106.2 106.3 0 0 0 0 0 0 0 0 0 0 0 12
13 106.4 101.6 106.4 106.3 106.2 106.2 1 0 0 0 0 0 0 0 0 0 0 13
14 106.5 101.4 106.5 106.4 106.3 106.2 0 1 0 0 0 0 0 0 0 0 0 14
15 106.6 100.8 106.6 106.5 106.4 106.3 0 0 1 0 0 0 0 0 0 0 0 15
16 106.6 101.1 106.6 106.6 106.5 106.4 0 0 0 1 0 0 0 0 0 0 0 16
17 106.6 101.3 106.6 106.6 106.6 106.5 0 0 0 0 1 0 0 0 0 0 0 17
18 106.8 101.2 106.8 106.6 106.6 106.6 0 0 0 0 0 1 0 0 0 0 0 18
19 107.0 101.3 107.0 106.8 106.6 106.6 0 0 0 0 0 0 1 0 0 0 0 19
20 107.2 101.1 107.2 107.0 106.8 106.6 0 0 0 0 0 0 0 1 0 0 0 20
21 107.3 101.3 107.3 107.2 107.0 106.8 0 0 0 0 0 0 0 0 1 0 0 21
22 107.5 101.2 107.5 107.3 107.2 107.0 0 0 0 0 0 0 0 0 0 1 0 22
23 107.6 101.6 107.6 107.5 107.3 107.2 0 0 0 0 0 0 0 0 0 0 1 23
24 107.6 101.7 107.6 107.6 107.5 107.3 0 0 0 0 0 0 0 0 0 0 0 24
25 107.7 101.5 107.7 107.6 107.6 107.5 1 0 0 0 0 0 0 0 0 0 0 25
26 107.7 100.9 107.7 107.7 107.6 107.6 0 1 0 0 0 0 0 0 0 0 0 26
27 107.7 101.5 107.7 107.7 107.7 107.6 0 0 1 0 0 0 0 0 0 0 0 27
28 107.7 101.4 107.7 107.7 107.7 107.7 0 0 0 1 0 0 0 0 0 0 0 28
29 107.6 101.6 107.6 107.7 107.7 107.7 0 0 0 0 1 0 0 0 0 0 0 29
30 107.7 101.7 107.7 107.6 107.7 107.7 0 0 0 0 0 1 0 0 0 0 0 30
31 107.9 101.4 107.9 107.7 107.6 107.7 0 0 0 0 0 0 1 0 0 0 0 31
32 107.9 101.8 107.9 107.9 107.7 107.6 0 0 0 0 0 0 0 1 0 0 0 32
33 107.9 101.7 107.9 107.9 107.9 107.7 0 0 0 0 0 0 0 0 1 0 0 33
34 107.8 101.4 107.8 107.9 107.9 107.9 0 0 0 0 0 0 0 0 0 1 0 34
35 107.6 101.2 107.6 107.8 107.9 107.9 0 0 0 0 0 0 0 0 0 0 1 35
36 107.4 101.0 107.4 107.6 107.8 107.9 0 0 0 0 0 0 0 0 0 0 0 36
37 107.0 101.7 107.0 107.4 107.6 107.8 1 0 0 0 0 0 0 0 0 0 0 37
38 107.0 102.4 107.0 107.0 107.4 107.6 0 1 0 0 0 0 0 0 0 0 0 38
39 107.2 102.0 107.2 107.0 107.0 107.4 0 0 1 0 0 0 0 0 0 0 0 39
40 107.5 102.1 107.5 107.2 107.0 107.0 0 0 0 1 0 0 0 0 0 0 0 40
41 107.8 102.0 107.8 107.5 107.2 107.0 0 0 0 0 1 0 0 0 0 0 0 41
42 107.8 101.8 107.8 107.8 107.5 107.2 0 0 0 0 0 1 0 0 0 0 0 42
43 107.7 102.7 107.7 107.8 107.8 107.5 0 0 0 0 0 0 1 0 0 0 0 43
44 107.6 102.3 107.6 107.7 107.8 107.8 0 0 0 0 0 0 0 1 0 0 0 44
45 107.6 101.9 107.6 107.6 107.7 107.8 0 0 0 0 0 0 0 0 1 0 0 45
46 107.5 102.0 107.5 107.6 107.6 107.7 0 0 0 0 0 0 0 0 0 1 0 46
47 107.5 102.3 107.5 107.5 107.6 107.6 0 0 0 0 0 0 0 0 0 0 1 47
48 107.6 102.8 107.6 107.5 107.5 107.6 0 0 0 0 0 0 0 0 0 0 0 48
49 107.6 102.4 107.6 107.6 107.5 107.5 1 0 0 0 0 0 0 0 0 0 0 49
50 107.9 102.3 107.9 107.6 107.6 107.5 0 1 0 0 0 0 0 0 0 0 0 50
51 107.6 102.7 107.6 107.9 107.6 107.6 0 0 1 0 0 0 0 0 0 0 0 51
52 107.5 102.7 107.5 107.6 107.9 107.6 0 0 0 1 0 0 0 0 0 0 0 52
53 107.5 102.9 107.5 107.5 107.6 107.9 0 0 0 0 1 0 0 0 0 0 0 53
54 107.6 103.0 107.6 107.5 107.5 107.6 0 0 0 0 0 1 0 0 0 0 0 54
55 107.7 102.2 107.7 107.6 107.5 107.5 0 0 0 0 0 0 1 0 0 0 0 55
56 107.8 102.3 107.8 107.7 107.6 107.5 0 0 0 0 0 0 0 1 0 0 0 56
57 107.9 102.8 107.9 107.8 107.7 107.6 0 0 0 0 0 0 0 0 1 0 0 57
58 107.9 102.8 107.9 107.9 107.8 107.7 0 0 0 0 0 0 0 0 0 1 0 58
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Infl `Yt-1` `Yt-2` `Yt-3` `Yt-4`
-2.986e-14 0.000e+00 1.000e+00 1.815e-16 -2.455e-16 7.414e-17
M1 M2 M3 M4 M5 M6
6.598e-17 9.081e-17 -2.666e-16 7.092e-17 1.637e-17 3.817e-17
M7 M8 M9 M10 M11 t
2.985e-17 -2.710e-18 6.796e-17 3.474e-17 5.463e-17 6.356e-18
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.087e-15 -5.305e-17 -6.422e-18 5.491e-17 3.547e-16
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.986e-14 1.762e-14 -1.695e+00 0.0979 .
Infl 0.000e+00 6.693e-17 0.000e+00 1.0000
`Yt-1` 1.000e+00 2.673e-16 3.742e+15 <2e-16 ***
`Yt-2` 1.815e-16 4.499e-16 4.030e-01 0.6887
`Yt-3` -2.455e-16 4.496e-16 -5.460e-01 0.5881
`Yt-4` 7.414e-17 2.639e-16 2.810e-01 0.7802
M1 6.598e-17 1.484e-16 4.450e-01 0.6591
M2 9.081e-17 1.500e-16 6.050e-01 0.5484
M3 -2.666e-16 1.507e-16 -1.769e+00 0.0845 .
M4 7.092e-17 1.511e-16 4.690e-01 0.6413
M5 1.637e-17 1.474e-16 1.110e-01 0.9121
M6 3.817e-17 1.500e-16 2.550e-01 0.8004
M7 2.985e-17 1.516e-16 1.970e-01 0.8449
M8 -2.710e-18 1.527e-16 -1.800e-02 0.9859
M9 6.796e-17 1.514e-16 4.490e-01 0.6559
M10 3.474e-17 1.471e-16 2.360e-01 0.8145
M11 5.463e-17 1.538e-16 3.550e-01 0.7243
t 6.356e-18 4.775e-18 1.331e+00 0.1907
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.171e-16 on 40 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 3.332e+31 on 17 and 40 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.5910697855 8.178604e-01 4.089302e-01
[2,] 0.9999999970 6.072225e-09 3.036113e-09
[3,] 0.0010222101 2.044420e-03 9.989778e-01
[4,] 0.9999975504 4.899193e-06 2.449597e-06
[5,] 0.6235077849 7.529844e-01 3.764922e-01
[6,] 0.0006504264 1.300853e-03 9.993496e-01
[7,] 0.8599394637 2.801211e-01 1.400605e-01
[8,] 0.0009543922 1.908784e-03 9.990456e-01
[9,] 1.0000000000 2.437190e-11 1.218595e-11
[10,] 0.8245564440 3.508871e-01 1.754436e-01
[11,] 0.0054653944 1.093079e-02 9.945346e-01
[12,] 0.9304293699 1.391413e-01 6.957063e-02
[13,] 0.9997316965 5.366070e-04 2.683035e-04
[14,] 0.9999995283 9.433821e-07 4.716910e-07
[15,] 0.9541822988 9.163540e-02 4.581770e-02
[16,] 0.8695621604 2.608757e-01 1.304378e-01
[17,] 0.9999238612 1.522775e-04 7.613875e-05
> postscript(file="/var/www/html/rcomp/tmp/1bvnr1258724081.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/2b2wz1258724081.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/3tl9o1258724081.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/4mr4z1258724081.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/575lh1258724081.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 58
Frequency = 1
1 2 3 4 5
1.495525e-16 3.339375e-16 -1.087200e-15 1.394266e-16 1.552651e-16
6 7 8 9 10
-5.195247e-17 6.004128e-17 5.163499e-17 3.652647e-17 3.876108e-17
11 12 13 14 15
9.528622e-17 5.599780e-17 3.896025e-17 -9.994550e-17 2.706027e-16
16 17 18 19 20
-8.716663e-18 -1.071866e-17 1.305002e-17 -2.927880e-17 -8.286291e-18
21 22 23 24 25
-5.900478e-17 1.600961e-17 -4.762680e-17 4.603022e-17 3.583665e-17
26 27 28 29 30
-9.647552e-17 3.546522e-16 -7.395779e-17 -2.266913e-17 -2.575028e-17
31 32 33 34 35
-5.083625e-17 -5.261027e-17 1.225801e-18 -9.286188e-17 -5.930945e-17
36 37 38 39 40
-6.517936e-17 -1.280720e-16 -2.351986e-17 2.101752e-16 -5.219528e-17
41 42 43 44 45
-1.314622e-16 2.309913e-17 1.214621e-16 1.266266e-16 -4.418155e-17
46 47 48 49 50
-5.319796e-17 1.165003e-17 -3.684866e-17 -9.627741e-17 -1.139967e-16
51 52 53 54 55
2.517698e-16 -4.556888e-18 9.584887e-18 4.155360e-17 -1.013884e-16
56 57 58
-1.173651e-16 6.543406e-17 9.128916e-17
> postscript(file="/var/www/html/rcomp/tmp/6pwkf1258724081.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 1.495525e-16 NA
1 3.339375e-16 1.495525e-16
2 -1.087200e-15 3.339375e-16
3 1.394266e-16 -1.087200e-15
4 1.552651e-16 1.394266e-16
5 -5.195247e-17 1.552651e-16
6 6.004128e-17 -5.195247e-17
7 5.163499e-17 6.004128e-17
8 3.652647e-17 5.163499e-17
9 3.876108e-17 3.652647e-17
10 9.528622e-17 3.876108e-17
11 5.599780e-17 9.528622e-17
12 3.896025e-17 5.599780e-17
13 -9.994550e-17 3.896025e-17
14 2.706027e-16 -9.994550e-17
15 -8.716663e-18 2.706027e-16
16 -1.071866e-17 -8.716663e-18
17 1.305002e-17 -1.071866e-17
18 -2.927880e-17 1.305002e-17
19 -8.286291e-18 -2.927880e-17
20 -5.900478e-17 -8.286291e-18
21 1.600961e-17 -5.900478e-17
22 -4.762680e-17 1.600961e-17
23 4.603022e-17 -4.762680e-17
24 3.583665e-17 4.603022e-17
25 -9.647552e-17 3.583665e-17
26 3.546522e-16 -9.647552e-17
27 -7.395779e-17 3.546522e-16
28 -2.266913e-17 -7.395779e-17
29 -2.575028e-17 -2.266913e-17
30 -5.083625e-17 -2.575028e-17
31 -5.261027e-17 -5.083625e-17
32 1.225801e-18 -5.261027e-17
33 -9.286188e-17 1.225801e-18
34 -5.930945e-17 -9.286188e-17
35 -6.517936e-17 -5.930945e-17
36 -1.280720e-16 -6.517936e-17
37 -2.351986e-17 -1.280720e-16
38 2.101752e-16 -2.351986e-17
39 -5.219528e-17 2.101752e-16
40 -1.314622e-16 -5.219528e-17
41 2.309913e-17 -1.314622e-16
42 1.214621e-16 2.309913e-17
43 1.266266e-16 1.214621e-16
44 -4.418155e-17 1.266266e-16
45 -5.319796e-17 -4.418155e-17
46 1.165003e-17 -5.319796e-17
47 -3.684866e-17 1.165003e-17
48 -9.627741e-17 -3.684866e-17
49 -1.139967e-16 -9.627741e-17
50 2.517698e-16 -1.139967e-16
51 -4.556888e-18 2.517698e-16
52 9.584887e-18 -4.556888e-18
53 4.155360e-17 9.584887e-18
54 -1.013884e-16 4.155360e-17
55 -1.173651e-16 -1.013884e-16
56 6.543406e-17 -1.173651e-16
57 9.128916e-17 6.543406e-17
58 NA 9.128916e-17
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.339375e-16 1.495525e-16
[2,] -1.087200e-15 3.339375e-16
[3,] 1.394266e-16 -1.087200e-15
[4,] 1.552651e-16 1.394266e-16
[5,] -5.195247e-17 1.552651e-16
[6,] 6.004128e-17 -5.195247e-17
[7,] 5.163499e-17 6.004128e-17
[8,] 3.652647e-17 5.163499e-17
[9,] 3.876108e-17 3.652647e-17
[10,] 9.528622e-17 3.876108e-17
[11,] 5.599780e-17 9.528622e-17
[12,] 3.896025e-17 5.599780e-17
[13,] -9.994550e-17 3.896025e-17
[14,] 2.706027e-16 -9.994550e-17
[15,] -8.716663e-18 2.706027e-16
[16,] -1.071866e-17 -8.716663e-18
[17,] 1.305002e-17 -1.071866e-17
[18,] -2.927880e-17 1.305002e-17
[19,] -8.286291e-18 -2.927880e-17
[20,] -5.900478e-17 -8.286291e-18
[21,] 1.600961e-17 -5.900478e-17
[22,] -4.762680e-17 1.600961e-17
[23,] 4.603022e-17 -4.762680e-17
[24,] 3.583665e-17 4.603022e-17
[25,] -9.647552e-17 3.583665e-17
[26,] 3.546522e-16 -9.647552e-17
[27,] -7.395779e-17 3.546522e-16
[28,] -2.266913e-17 -7.395779e-17
[29,] -2.575028e-17 -2.266913e-17
[30,] -5.083625e-17 -2.575028e-17
[31,] -5.261027e-17 -5.083625e-17
[32,] 1.225801e-18 -5.261027e-17
[33,] -9.286188e-17 1.225801e-18
[34,] -5.930945e-17 -9.286188e-17
[35,] -6.517936e-17 -5.930945e-17
[36,] -1.280720e-16 -6.517936e-17
[37,] -2.351986e-17 -1.280720e-16
[38,] 2.101752e-16 -2.351986e-17
[39,] -5.219528e-17 2.101752e-16
[40,] -1.314622e-16 -5.219528e-17
[41,] 2.309913e-17 -1.314622e-16
[42,] 1.214621e-16 2.309913e-17
[43,] 1.266266e-16 1.214621e-16
[44,] -4.418155e-17 1.266266e-16
[45,] -5.319796e-17 -4.418155e-17
[46,] 1.165003e-17 -5.319796e-17
[47,] -3.684866e-17 1.165003e-17
[48,] -9.627741e-17 -3.684866e-17
[49,] -1.139967e-16 -9.627741e-17
[50,] 2.517698e-16 -1.139967e-16
[51,] -4.556888e-18 2.517698e-16
[52,] 9.584887e-18 -4.556888e-18
[53,] 4.155360e-17 9.584887e-18
[54,] -1.013884e-16 4.155360e-17
[55,] -1.173651e-16 -1.013884e-16
[56,] 6.543406e-17 -1.173651e-16
[57,] 9.128916e-17 6.543406e-17
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.339375e-16 1.495525e-16
2 -1.087200e-15 3.339375e-16
3 1.394266e-16 -1.087200e-15
4 1.552651e-16 1.394266e-16
5 -5.195247e-17 1.552651e-16
6 6.004128e-17 -5.195247e-17
7 5.163499e-17 6.004128e-17
8 3.652647e-17 5.163499e-17
9 3.876108e-17 3.652647e-17
10 9.528622e-17 3.876108e-17
11 5.599780e-17 9.528622e-17
12 3.896025e-17 5.599780e-17
13 -9.994550e-17 3.896025e-17
14 2.706027e-16 -9.994550e-17
15 -8.716663e-18 2.706027e-16
16 -1.071866e-17 -8.716663e-18
17 1.305002e-17 -1.071866e-17
18 -2.927880e-17 1.305002e-17
19 -8.286291e-18 -2.927880e-17
20 -5.900478e-17 -8.286291e-18
21 1.600961e-17 -5.900478e-17
22 -4.762680e-17 1.600961e-17
23 4.603022e-17 -4.762680e-17
24 3.583665e-17 4.603022e-17
25 -9.647552e-17 3.583665e-17
26 3.546522e-16 -9.647552e-17
27 -7.395779e-17 3.546522e-16
28 -2.266913e-17 -7.395779e-17
29 -2.575028e-17 -2.266913e-17
30 -5.083625e-17 -2.575028e-17
31 -5.261027e-17 -5.083625e-17
32 1.225801e-18 -5.261027e-17
33 -9.286188e-17 1.225801e-18
34 -5.930945e-17 -9.286188e-17
35 -6.517936e-17 -5.930945e-17
36 -1.280720e-16 -6.517936e-17
37 -2.351986e-17 -1.280720e-16
38 2.101752e-16 -2.351986e-17
39 -5.219528e-17 2.101752e-16
40 -1.314622e-16 -5.219528e-17
41 2.309913e-17 -1.314622e-16
42 1.214621e-16 2.309913e-17
43 1.266266e-16 1.214621e-16
44 -4.418155e-17 1.266266e-16
45 -5.319796e-17 -4.418155e-17
46 1.165003e-17 -5.319796e-17
47 -3.684866e-17 1.165003e-17
48 -9.627741e-17 -3.684866e-17
49 -1.139967e-16 -9.627741e-17
50 2.517698e-16 -1.139967e-16
51 -4.556888e-18 2.517698e-16
52 9.584887e-18 -4.556888e-18
53 4.155360e-17 9.584887e-18
54 -1.013884e-16 4.155360e-17
55 -1.173651e-16 -1.013884e-16
56 6.543406e-17 -1.173651e-16
57 9.128916e-17 6.543406e-17
> 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/7yf4o1258724081.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/8pxgq1258724081.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/9u1s51258724081.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/10cr4y1258724081.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/11opsf1258724081.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/12svg71258724081.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/137yc21258724081.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/14lbr01258724081.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/15mo471258724081.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/16ypm41258724081.tab")
+ }
>
> system("convert tmp/1bvnr1258724081.ps tmp/1bvnr1258724081.png")
> system("convert tmp/2b2wz1258724081.ps tmp/2b2wz1258724081.png")
> system("convert tmp/3tl9o1258724081.ps tmp/3tl9o1258724081.png")
> system("convert tmp/4mr4z1258724081.ps tmp/4mr4z1258724081.png")
> system("convert tmp/575lh1258724081.ps tmp/575lh1258724081.png")
> system("convert tmp/6pwkf1258724081.ps tmp/6pwkf1258724081.png")
> system("convert tmp/7yf4o1258724081.ps tmp/7yf4o1258724081.png")
> system("convert tmp/8pxgq1258724081.ps tmp/8pxgq1258724081.png")
> system("convert tmp/9u1s51258724081.ps tmp/9u1s51258724081.png")
> system("convert tmp/10cr4y1258724081.ps tmp/10cr4y1258724081.png")
>
>
> proc.time()
user system elapsed
2.340 1.534 2.761