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.08
+ ,1.00
+ ,2.05
+ ,2.09
+ ,2.06
+ ,1.00
+ ,2.08
+ ,2.05
+ ,2.06
+ ,1.00
+ ,2.06
+ ,2.08
+ ,2.08
+ ,1.00
+ ,2.06
+ ,2.06
+ ,2.07
+ ,1.00
+ ,2.08
+ ,2.06
+ ,2.06
+ ,1.00
+ ,2.07
+ ,2.08
+ ,2.07
+ ,1.00
+ ,2.06
+ ,2.07
+ ,2.06
+ ,1.00
+ ,2.07
+ ,2.06
+ ,2.09
+ ,1.00
+ ,2.06
+ ,2.07
+ ,2.07
+ ,1.00
+ ,2.09
+ ,2.06
+ ,2.09
+ ,1.00
+ ,2.07
+ ,2.09
+ ,2.28
+ ,1.25
+ ,2.09
+ ,2.07
+ ,2.33
+ ,1.25
+ ,2.28
+ ,2.09
+ ,2.35
+ ,1.25
+ ,2.33
+ ,2.28
+ ,2.52
+ ,1.50
+ ,2.35
+ ,2.33
+ ,2.63
+ ,1.50
+ ,2.52
+ ,2.35
+ ,2.58
+ ,1.50
+ ,2.63
+ ,2.52
+ ,2.70
+ ,1.75
+ ,2.58
+ ,2.63
+ ,2.81
+ ,1.75
+ ,2.70
+ ,2.58
+ ,2.97
+ ,2.00
+ ,2.81
+ ,2.70
+ ,3.04
+ ,2.00
+ ,2.97
+ ,2.81
+ ,3.28
+ ,2.25
+ ,3.04
+ ,2.97
+ ,3.33
+ ,2.25
+ ,3.28
+ ,3.04
+ ,3.50
+ ,2.50
+ ,3.33
+ ,3.28
+ ,3.56
+ ,2.50
+ ,3.50
+ ,3.33
+ ,3.57
+ ,2.50
+ ,3.56
+ ,3.50
+ ,3.69
+ ,2.75
+ ,3.57
+ ,3.56
+ ,3.82
+ ,2.75
+ ,3.69
+ ,3.57
+ ,3.79
+ ,2.75
+ ,3.82
+ ,3.69
+ ,3.96
+ ,3.00
+ ,3.79
+ ,3.82
+ ,4.06
+ ,3.00
+ ,3.96
+ ,3.79
+ ,4.05
+ ,3.00
+ ,4.06
+ ,3.96
+ ,4.03
+ ,3.00
+ ,4.05
+ ,4.06
+ ,3.94
+ ,3.00
+ ,4.03
+ ,4.05
+ ,4.02
+ ,3.00
+ ,3.94
+ ,4.03
+ ,3.88
+ ,3.00
+ ,4.02
+ ,3.94
+ ,4.02
+ ,3.00
+ ,3.88
+ ,4.02
+ ,4.03
+ ,3.00
+ ,4.02
+ ,3.88
+ ,4.09
+ ,3.00
+ ,4.03
+ ,4.02
+ ,3.99
+ ,3.00
+ ,4.09
+ ,4.03
+ ,4.01
+ ,3.00
+ ,3.99
+ ,4.09
+ ,4.01
+ ,3.00
+ ,4.01
+ ,3.99
+ ,4.19
+ ,3.25
+ ,4.01
+ ,4.01
+ ,4.30
+ ,3.25
+ ,4.19
+ ,4.01
+ ,4.27
+ ,3.25
+ ,4.30
+ ,4.19
+ ,3.82
+ ,3.25
+ ,4.27
+ ,4.30
+ ,3.15
+ ,2.75
+ ,3.82
+ ,4.27
+ ,2.49
+ ,2.00
+ ,3.15
+ ,3.82
+ ,1.81
+ ,1.00
+ ,2.49
+ ,3.15
+ ,1.26
+ ,1.00
+ ,1.81
+ ,2.49
+ ,1.06
+ ,0.50
+ ,1.26
+ ,1.81
+ ,0.84
+ ,0.25
+ ,1.06
+ ,1.26
+ ,0.78
+ ,0.25
+ ,0.84
+ ,1.06
+ ,0.70
+ ,0.25
+ ,0.78
+ ,0.84
+ ,0.36
+ ,0.25
+ ,0.70
+ ,0.78
+ ,0.35
+ ,0.25
+ ,0.36
+ ,0.70)
+ ,dim=c(4
+ ,56)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y-1'
+ ,'Y-2')
+ ,1:56))
> y <- array(NA,dim=c(4,56),dimnames=list(c('Y','X','Y-1','Y-2'),1:56))
> 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 Y-1 Y-2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 2.08 1.00 2.05 2.09 1 0 0 0 0 0 0 0 0 0 0 1
2 2.06 1.00 2.08 2.05 0 1 0 0 0 0 0 0 0 0 0 2
3 2.06 1.00 2.06 2.08 0 0 1 0 0 0 0 0 0 0 0 3
4 2.08 1.00 2.06 2.06 0 0 0 1 0 0 0 0 0 0 0 4
5 2.07 1.00 2.08 2.06 0 0 0 0 1 0 0 0 0 0 0 5
6 2.06 1.00 2.07 2.08 0 0 0 0 0 1 0 0 0 0 0 6
7 2.07 1.00 2.06 2.07 0 0 0 0 0 0 1 0 0 0 0 7
8 2.06 1.00 2.07 2.06 0 0 0 0 0 0 0 1 0 0 0 8
9 2.09 1.00 2.06 2.07 0 0 0 0 0 0 0 0 1 0 0 9
10 2.07 1.00 2.09 2.06 0 0 0 0 0 0 0 0 0 1 0 10
11 2.09 1.00 2.07 2.09 0 0 0 0 0 0 0 0 0 0 1 11
12 2.28 1.25 2.09 2.07 0 0 0 0 0 0 0 0 0 0 0 12
13 2.33 1.25 2.28 2.09 1 0 0 0 0 0 0 0 0 0 0 13
14 2.35 1.25 2.33 2.28 0 1 0 0 0 0 0 0 0 0 0 14
15 2.52 1.50 2.35 2.33 0 0 1 0 0 0 0 0 0 0 0 15
16 2.63 1.50 2.52 2.35 0 0 0 1 0 0 0 0 0 0 0 16
17 2.58 1.50 2.63 2.52 0 0 0 0 1 0 0 0 0 0 0 17
18 2.70 1.75 2.58 2.63 0 0 0 0 0 1 0 0 0 0 0 18
19 2.81 1.75 2.70 2.58 0 0 0 0 0 0 1 0 0 0 0 19
20 2.97 2.00 2.81 2.70 0 0 0 0 0 0 0 1 0 0 0 20
21 3.04 2.00 2.97 2.81 0 0 0 0 0 0 0 0 1 0 0 21
22 3.28 2.25 3.04 2.97 0 0 0 0 0 0 0 0 0 1 0 22
23 3.33 2.25 3.28 3.04 0 0 0 0 0 0 0 0 0 0 1 23
24 3.50 2.50 3.33 3.28 0 0 0 0 0 0 0 0 0 0 0 24
25 3.56 2.50 3.50 3.33 1 0 0 0 0 0 0 0 0 0 0 25
26 3.57 2.50 3.56 3.50 0 1 0 0 0 0 0 0 0 0 0 26
27 3.69 2.75 3.57 3.56 0 0 1 0 0 0 0 0 0 0 0 27
28 3.82 2.75 3.69 3.57 0 0 0 1 0 0 0 0 0 0 0 28
29 3.79 2.75 3.82 3.69 0 0 0 0 1 0 0 0 0 0 0 29
30 3.96 3.00 3.79 3.82 0 0 0 0 0 1 0 0 0 0 0 30
31 4.06 3.00 3.96 3.79 0 0 0 0 0 0 1 0 0 0 0 31
32 4.05 3.00 4.06 3.96 0 0 0 0 0 0 0 1 0 0 0 32
33 4.03 3.00 4.05 4.06 0 0 0 0 0 0 0 0 1 0 0 33
34 3.94 3.00 4.03 4.05 0 0 0 0 0 0 0 0 0 1 0 34
35 4.02 3.00 3.94 4.03 0 0 0 0 0 0 0 0 0 0 1 35
36 3.88 3.00 4.02 3.94 0 0 0 0 0 0 0 0 0 0 0 36
37 4.02 3.00 3.88 4.02 1 0 0 0 0 0 0 0 0 0 0 37
38 4.03 3.00 4.02 3.88 0 1 0 0 0 0 0 0 0 0 0 38
39 4.09 3.00 4.03 4.02 0 0 1 0 0 0 0 0 0 0 0 39
40 3.99 3.00 4.09 4.03 0 0 0 1 0 0 0 0 0 0 0 40
41 4.01 3.00 3.99 4.09 0 0 0 0 1 0 0 0 0 0 0 41
42 4.01 3.00 4.01 3.99 0 0 0 0 0 1 0 0 0 0 0 42
43 4.19 3.25 4.01 4.01 0 0 0 0 0 0 1 0 0 0 0 43
44 4.30 3.25 4.19 4.01 0 0 0 0 0 0 0 1 0 0 0 44
45 4.27 3.25 4.30 4.19 0 0 0 0 0 0 0 0 1 0 0 45
46 3.82 3.25 4.27 4.30 0 0 0 0 0 0 0 0 0 1 0 46
47 3.15 2.75 3.82 4.27 0 0 0 0 0 0 0 0 0 0 1 47
48 2.49 2.00 3.15 3.82 0 0 0 0 0 0 0 0 0 0 0 48
49 1.81 1.00 2.49 3.15 1 0 0 0 0 0 0 0 0 0 0 49
50 1.26 1.00 1.81 2.49 0 1 0 0 0 0 0 0 0 0 0 50
51 1.06 0.50 1.26 1.81 0 0 1 0 0 0 0 0 0 0 0 51
52 0.84 0.25 1.06 1.26 0 0 0 1 0 0 0 0 0 0 0 52
53 0.78 0.25 0.84 1.06 0 0 0 0 1 0 0 0 0 0 0 53
54 0.70 0.25 0.78 0.84 0 0 0 0 0 1 0 0 0 0 0 54
55 0.36 0.25 0.70 0.78 0 0 0 0 0 0 1 0 0 0 0 55
56 0.35 0.25 0.36 0.70 0 0 0 0 0 0 0 1 0 0 0 56
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X `Y-1` `Y-2` M1 M2
0.599181 0.690171 0.799965 -0.415248 0.090575 0.016319
M3 M4 M5 M6 M7 M8
0.105507 0.067612 0.071280 0.065691 0.007998 0.036112
M9 M10 M11 t
0.054978 -0.044593 -0.011523 -0.007611
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.260647 -0.046232 0.002773 0.058112 0.149818
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.599181 0.116311 5.152 7.30e-06 ***
X 0.690171 0.117726 5.863 7.38e-07 ***
`Y-1` 0.799965 0.179596 4.454 6.61e-05 ***
`Y-2` -0.415248 0.110474 -3.759 0.000546 ***
M1 0.090575 0.073231 1.237 0.223363
M2 0.016319 0.072225 0.226 0.822390
M3 0.105507 0.071539 1.475 0.148091
M4 0.067612 0.075284 0.898 0.374511
M5 0.071280 0.074097 0.962 0.341838
M6 0.065691 0.072807 0.902 0.372322
M7 0.007998 0.074430 0.107 0.914959
M8 0.036112 0.073909 0.489 0.627798
M9 0.054978 0.077543 0.709 0.482437
M10 -0.044593 0.075920 -0.587 0.560257
M11 -0.011523 0.075387 -0.153 0.879288
t -0.007611 0.001516 -5.020 1.11e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1059 on 40 degrees of freedom
Multiple R-squared: 0.9935, Adjusted R-squared: 0.991
F-statistic: 406.8 on 15 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,] 1.433843e-02 2.867687e-02 0.9856616
[2,] 2.529300e-03 5.058600e-03 0.9974707
[3,] 4.688397e-04 9.376794e-04 0.9995312
[4,] 5.129493e-04 1.025899e-03 0.9994871
[5,] 9.588193e-05 1.917639e-04 0.9999041
[6,] 2.085093e-05 4.170186e-05 0.9999791
[7,] 4.625651e-06 9.251302e-06 0.9999954
[8,] 7.492885e-07 1.498577e-06 0.9999993
[9,] 6.841263e-07 1.368253e-06 0.9999993
[10,] 1.295713e-07 2.591426e-07 0.9999999
[11,] 6.430246e-08 1.286049e-07 0.9999999
[12,] 2.038335e-08 4.076671e-08 1.0000000
[13,] 4.249733e-09 8.499466e-09 1.0000000
[14,] 2.802443e-09 5.604887e-09 1.0000000
[15,] 7.761248e-09 1.552250e-08 1.0000000
[16,] 2.443151e-07 4.886302e-07 0.9999998
[17,] 5.043404e-07 1.008681e-06 0.9999995
[18,] 1.259466e-05 2.518933e-05 0.9999874
[19,] 3.738952e-06 7.477904e-06 0.9999963
> postscript(file="/var/www/html/rcomp/tmp/17pqo1258739549.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/2zzhz1258739549.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/3afz01258739549.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/46g7c1258739549.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/5t4il1258739549.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 = 56
Frequency = 1
1 2 3 4 5 6
-0.064373440 -0.043115648 -0.096235294 -0.039033537 -0.061089856 -0.041584268
7 8 9 10 11 12
0.037566477 -0.005087472 0.025809995 0.084840783 0.107838403 0.097080240
13 14 15 16 17 18
-0.079571243 0.061194369 -0.018161588 0.009656108 -0.053804817 -0.007470984
19 20 21 22 23 24
0.051074421 -0.020136403 -0.043708097 0.141373628 0.002990367 0.056197826
25 26 27 28 29 30
-0.081996914 0.032464085 -0.084739742 -0.001076295 -0.081298930 0.007340576
31 32 33 34 35 36
0.024192714 -0.015713357 0.002556459 0.031585480 0.149818211 -0.095462525
37 38 39 40 41 42
0.106789232 0.028526079 0.057084807 -0.041253866 0.087600474 0.043277324
43 44 45 46 47 48
0.124343181 0.069847722 0.015341643 -0.257799891 -0.260646981 -0.057815540
49 50 51 52 53 54
0.119152366 -0.079068885 0.142051817 0.071707590 0.108593127 -0.001562648
55 56
-0.237176793 -0.028910490
> postscript(file="/var/www/html/rcomp/tmp/6n7w21258739549.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 = 56
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.064373440 NA
1 -0.043115648 -0.064373440
2 -0.096235294 -0.043115648
3 -0.039033537 -0.096235294
4 -0.061089856 -0.039033537
5 -0.041584268 -0.061089856
6 0.037566477 -0.041584268
7 -0.005087472 0.037566477
8 0.025809995 -0.005087472
9 0.084840783 0.025809995
10 0.107838403 0.084840783
11 0.097080240 0.107838403
12 -0.079571243 0.097080240
13 0.061194369 -0.079571243
14 -0.018161588 0.061194369
15 0.009656108 -0.018161588
16 -0.053804817 0.009656108
17 -0.007470984 -0.053804817
18 0.051074421 -0.007470984
19 -0.020136403 0.051074421
20 -0.043708097 -0.020136403
21 0.141373628 -0.043708097
22 0.002990367 0.141373628
23 0.056197826 0.002990367
24 -0.081996914 0.056197826
25 0.032464085 -0.081996914
26 -0.084739742 0.032464085
27 -0.001076295 -0.084739742
28 -0.081298930 -0.001076295
29 0.007340576 -0.081298930
30 0.024192714 0.007340576
31 -0.015713357 0.024192714
32 0.002556459 -0.015713357
33 0.031585480 0.002556459
34 0.149818211 0.031585480
35 -0.095462525 0.149818211
36 0.106789232 -0.095462525
37 0.028526079 0.106789232
38 0.057084807 0.028526079
39 -0.041253866 0.057084807
40 0.087600474 -0.041253866
41 0.043277324 0.087600474
42 0.124343181 0.043277324
43 0.069847722 0.124343181
44 0.015341643 0.069847722
45 -0.257799891 0.015341643
46 -0.260646981 -0.257799891
47 -0.057815540 -0.260646981
48 0.119152366 -0.057815540
49 -0.079068885 0.119152366
50 0.142051817 -0.079068885
51 0.071707590 0.142051817
52 0.108593127 0.071707590
53 -0.001562648 0.108593127
54 -0.237176793 -0.001562648
55 -0.028910490 -0.237176793
56 NA -0.028910490
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.043115648 -0.064373440
[2,] -0.096235294 -0.043115648
[3,] -0.039033537 -0.096235294
[4,] -0.061089856 -0.039033537
[5,] -0.041584268 -0.061089856
[6,] 0.037566477 -0.041584268
[7,] -0.005087472 0.037566477
[8,] 0.025809995 -0.005087472
[9,] 0.084840783 0.025809995
[10,] 0.107838403 0.084840783
[11,] 0.097080240 0.107838403
[12,] -0.079571243 0.097080240
[13,] 0.061194369 -0.079571243
[14,] -0.018161588 0.061194369
[15,] 0.009656108 -0.018161588
[16,] -0.053804817 0.009656108
[17,] -0.007470984 -0.053804817
[18,] 0.051074421 -0.007470984
[19,] -0.020136403 0.051074421
[20,] -0.043708097 -0.020136403
[21,] 0.141373628 -0.043708097
[22,] 0.002990367 0.141373628
[23,] 0.056197826 0.002990367
[24,] -0.081996914 0.056197826
[25,] 0.032464085 -0.081996914
[26,] -0.084739742 0.032464085
[27,] -0.001076295 -0.084739742
[28,] -0.081298930 -0.001076295
[29,] 0.007340576 -0.081298930
[30,] 0.024192714 0.007340576
[31,] -0.015713357 0.024192714
[32,] 0.002556459 -0.015713357
[33,] 0.031585480 0.002556459
[34,] 0.149818211 0.031585480
[35,] -0.095462525 0.149818211
[36,] 0.106789232 -0.095462525
[37,] 0.028526079 0.106789232
[38,] 0.057084807 0.028526079
[39,] -0.041253866 0.057084807
[40,] 0.087600474 -0.041253866
[41,] 0.043277324 0.087600474
[42,] 0.124343181 0.043277324
[43,] 0.069847722 0.124343181
[44,] 0.015341643 0.069847722
[45,] -0.257799891 0.015341643
[46,] -0.260646981 -0.257799891
[47,] -0.057815540 -0.260646981
[48,] 0.119152366 -0.057815540
[49,] -0.079068885 0.119152366
[50,] 0.142051817 -0.079068885
[51,] 0.071707590 0.142051817
[52,] 0.108593127 0.071707590
[53,] -0.001562648 0.108593127
[54,] -0.237176793 -0.001562648
[55,] -0.028910490 -0.237176793
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.043115648 -0.064373440
2 -0.096235294 -0.043115648
3 -0.039033537 -0.096235294
4 -0.061089856 -0.039033537
5 -0.041584268 -0.061089856
6 0.037566477 -0.041584268
7 -0.005087472 0.037566477
8 0.025809995 -0.005087472
9 0.084840783 0.025809995
10 0.107838403 0.084840783
11 0.097080240 0.107838403
12 -0.079571243 0.097080240
13 0.061194369 -0.079571243
14 -0.018161588 0.061194369
15 0.009656108 -0.018161588
16 -0.053804817 0.009656108
17 -0.007470984 -0.053804817
18 0.051074421 -0.007470984
19 -0.020136403 0.051074421
20 -0.043708097 -0.020136403
21 0.141373628 -0.043708097
22 0.002990367 0.141373628
23 0.056197826 0.002990367
24 -0.081996914 0.056197826
25 0.032464085 -0.081996914
26 -0.084739742 0.032464085
27 -0.001076295 -0.084739742
28 -0.081298930 -0.001076295
29 0.007340576 -0.081298930
30 0.024192714 0.007340576
31 -0.015713357 0.024192714
32 0.002556459 -0.015713357
33 0.031585480 0.002556459
34 0.149818211 0.031585480
35 -0.095462525 0.149818211
36 0.106789232 -0.095462525
37 0.028526079 0.106789232
38 0.057084807 0.028526079
39 -0.041253866 0.057084807
40 0.087600474 -0.041253866
41 0.043277324 0.087600474
42 0.124343181 0.043277324
43 0.069847722 0.124343181
44 0.015341643 0.069847722
45 -0.257799891 0.015341643
46 -0.260646981 -0.257799891
47 -0.057815540 -0.260646981
48 0.119152366 -0.057815540
49 -0.079068885 0.119152366
50 0.142051817 -0.079068885
51 0.071707590 0.142051817
52 0.108593127 0.071707590
53 -0.001562648 0.108593127
54 -0.237176793 -0.001562648
55 -0.028910490 -0.237176793
> 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/78ka71258739549.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/88cwu1258739549.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/9l2sf1258739549.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/1003gn1258739549.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/1154bo1258739549.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/12hdh81258739549.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/13ibh01258739549.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/14mfyt1258739549.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/15t94e1258739549.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/16ol4d1258739549.tab")
+ }
>
> system("convert tmp/17pqo1258739549.ps tmp/17pqo1258739549.png")
> system("convert tmp/2zzhz1258739549.ps tmp/2zzhz1258739549.png")
> system("convert tmp/3afz01258739549.ps tmp/3afz01258739549.png")
> system("convert tmp/46g7c1258739549.ps tmp/46g7c1258739549.png")
> system("convert tmp/5t4il1258739549.ps tmp/5t4il1258739549.png")
> system("convert tmp/6n7w21258739549.ps tmp/6n7w21258739549.png")
> system("convert tmp/78ka71258739549.ps tmp/78ka71258739549.png")
> system("convert tmp/88cwu1258739549.ps tmp/88cwu1258739549.png")
> system("convert tmp/9l2sf1258739549.ps tmp/9l2sf1258739549.png")
> system("convert tmp/1003gn1258739549.ps tmp/1003gn1258739549.png")
>
>
> proc.time()
user system elapsed
2.332 1.563 2.784