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(12.8
+ ,23
+ ,20.3
+ ,13.2
+ ,15.7
+ ,12.6
+ ,8
+ ,20
+ ,12.8
+ ,20.3
+ ,13.2
+ ,15.7
+ ,0.9
+ ,20
+ ,8
+ ,12.8
+ ,20.3
+ ,13.2
+ ,3.6
+ ,15
+ ,0.9
+ ,8
+ ,12.8
+ ,20.3
+ ,14.1
+ ,17
+ ,3.6
+ ,0.9
+ ,8
+ ,12.8
+ ,21.7
+ ,16
+ ,14.1
+ ,3.6
+ ,0.9
+ ,8
+ ,24.5
+ ,15
+ ,21.7
+ ,14.1
+ ,3.6
+ ,0.9
+ ,18.9
+ ,10
+ ,24.5
+ ,21.7
+ ,14.1
+ ,3.6
+ ,13.9
+ ,13
+ ,18.9
+ ,24.5
+ ,21.7
+ ,14.1
+ ,11
+ ,10
+ ,13.9
+ ,18.9
+ ,24.5
+ ,21.7
+ ,5.8
+ ,19
+ ,11
+ ,13.9
+ ,18.9
+ ,24.5
+ ,15.5
+ ,21
+ ,5.8
+ ,11
+ ,13.9
+ ,18.9
+ ,22.4
+ ,17
+ ,15.5
+ ,5.8
+ ,11
+ ,13.9
+ ,31.7
+ ,16
+ ,22.4
+ ,15.5
+ ,5.8
+ ,11
+ ,30.3
+ ,17
+ ,31.7
+ ,22.4
+ ,15.5
+ ,5.8
+ ,31.4
+ ,14
+ ,30.3
+ ,31.7
+ ,22.4
+ ,15.5
+ ,20.2
+ ,18
+ ,31.4
+ ,30.3
+ ,31.7
+ ,22.4
+ ,19.7
+ ,17
+ ,20.2
+ ,31.4
+ ,30.3
+ ,31.7
+ ,10.8
+ ,14
+ ,19.7
+ ,20.2
+ ,31.4
+ ,30.3
+ ,13.2
+ ,15
+ ,10.8
+ ,19.7
+ ,20.2
+ ,31.4
+ ,15.1
+ ,16
+ ,13.2
+ ,10.8
+ ,19.7
+ ,20.2
+ ,15.6
+ ,11
+ ,15.1
+ ,13.2
+ ,10.8
+ ,19.7
+ ,15.5
+ ,15
+ ,15.6
+ ,15.1
+ ,13.2
+ ,10.8
+ ,12.7
+ ,13
+ ,15.5
+ ,15.6
+ ,15.1
+ ,13.2
+ ,10.9
+ ,17
+ ,12.7
+ ,15.5
+ ,15.6
+ ,15.1
+ ,10
+ ,16
+ ,10.9
+ ,12.7
+ ,15.5
+ ,15.6
+ ,9.1
+ ,9
+ ,10
+ ,10.9
+ ,12.7
+ ,15.5
+ ,10.3
+ ,17
+ ,9.1
+ ,10
+ ,10.9
+ ,12.7
+ ,16.9
+ ,15
+ ,10.3
+ ,9.1
+ ,10
+ ,10.9
+ ,22
+ ,12
+ ,16.9
+ ,10.3
+ ,9.1
+ ,10
+ ,27.6
+ ,12
+ ,22
+ ,16.9
+ ,10.3
+ ,9.1
+ ,28.9
+ ,12
+ ,27.6
+ ,22
+ ,16.9
+ ,10.3
+ ,31
+ ,12
+ ,28.9
+ ,27.6
+ ,22
+ ,16.9
+ ,32.9
+ ,4
+ ,31
+ ,28.9
+ ,27.6
+ ,22
+ ,38.1
+ ,7
+ ,32.9
+ ,31
+ ,28.9
+ ,27.6
+ ,28.8
+ ,4
+ ,38.1
+ ,32.9
+ ,31
+ ,28.9
+ ,29
+ ,3
+ ,28.8
+ ,38.1
+ ,32.9
+ ,31
+ ,21.8
+ ,3
+ ,29
+ ,28.8
+ ,38.1
+ ,32.9
+ ,28.8
+ ,0
+ ,21.8
+ ,29
+ ,28.8
+ ,38.1
+ ,25.6
+ ,5
+ ,28.8
+ ,21.8
+ ,29
+ ,28.8
+ ,28.2
+ ,3
+ ,25.6
+ ,28.8
+ ,21.8
+ ,29
+ ,20.2
+ ,4
+ ,28.2
+ ,25.6
+ ,28.8
+ ,21.8
+ ,17.9
+ ,3
+ ,20.2
+ ,28.2
+ ,25.6
+ ,28.8
+ ,16.3
+ ,10
+ ,17.9
+ ,20.2
+ ,28.2
+ ,25.6
+ ,13.2
+ ,4
+ ,16.3
+ ,17.9
+ ,20.2
+ ,28.2
+ ,8.1
+ ,1
+ ,13.2
+ ,16.3
+ ,17.9
+ ,20.2
+ ,4.5
+ ,1
+ ,8.1
+ ,13.2
+ ,16.3
+ ,17.9
+ ,-0.1
+ ,8
+ ,4.5
+ ,8.1
+ ,13.2
+ ,16.3
+ ,0
+ ,5
+ ,-0.1
+ ,4.5
+ ,8.1
+ ,13.2
+ ,2.3
+ ,4
+ ,0
+ ,-0.1
+ ,4.5
+ ,8.1
+ ,2.8
+ ,0
+ ,2.3
+ ,0
+ ,-0.1
+ ,4.5
+ ,2.9
+ ,2
+ ,2.8
+ ,2.3
+ ,0
+ ,-0.1
+ ,0.1
+ ,7
+ ,2.9
+ ,2.8
+ ,2.3
+ ,0
+ ,3.5
+ ,6
+ ,0.1
+ ,2.9
+ ,2.8
+ ,2.3
+ ,8.6
+ ,9
+ ,3.5
+ ,0.1
+ ,2.9
+ ,2.8
+ ,13.8
+ ,10
+ ,8.6
+ ,3.5
+ ,0.1
+ ,2.9)
+ ,dim=c(6
+ ,56)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3'
+ ,'Y4')
+ ,1:56))
> y <- array(NA,dim=c(6,56),dimnames=list(c('Y','X','Y1','Y2','Y3','Y4'),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 Y1 Y2 Y3 Y4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 12.8 23 20.3 13.2 15.7 12.6 1 0 0 0 0 0 0 0 0 0 0 1
2 8.0 20 12.8 20.3 13.2 15.7 0 1 0 0 0 0 0 0 0 0 0 2
3 0.9 20 8.0 12.8 20.3 13.2 0 0 1 0 0 0 0 0 0 0 0 3
4 3.6 15 0.9 8.0 12.8 20.3 0 0 0 1 0 0 0 0 0 0 0 4
5 14.1 17 3.6 0.9 8.0 12.8 0 0 0 0 1 0 0 0 0 0 0 5
6 21.7 16 14.1 3.6 0.9 8.0 0 0 0 0 0 1 0 0 0 0 0 6
7 24.5 15 21.7 14.1 3.6 0.9 0 0 0 0 0 0 1 0 0 0 0 7
8 18.9 10 24.5 21.7 14.1 3.6 0 0 0 0 0 0 0 1 0 0 0 8
9 13.9 13 18.9 24.5 21.7 14.1 0 0 0 0 0 0 0 0 1 0 0 9
10 11.0 10 13.9 18.9 24.5 21.7 0 0 0 0 0 0 0 0 0 1 0 10
11 5.8 19 11.0 13.9 18.9 24.5 0 0 0 0 0 0 0 0 0 0 1 11
12 15.5 21 5.8 11.0 13.9 18.9 0 0 0 0 0 0 0 0 0 0 0 12
13 22.4 17 15.5 5.8 11.0 13.9 1 0 0 0 0 0 0 0 0 0 0 13
14 31.7 16 22.4 15.5 5.8 11.0 0 1 0 0 0 0 0 0 0 0 0 14
15 30.3 17 31.7 22.4 15.5 5.8 0 0 1 0 0 0 0 0 0 0 0 15
16 31.4 14 30.3 31.7 22.4 15.5 0 0 0 1 0 0 0 0 0 0 0 16
17 20.2 18 31.4 30.3 31.7 22.4 0 0 0 0 1 0 0 0 0 0 0 17
18 19.7 17 20.2 31.4 30.3 31.7 0 0 0 0 0 1 0 0 0 0 0 18
19 10.8 14 19.7 20.2 31.4 30.3 0 0 0 0 0 0 1 0 0 0 0 19
20 13.2 15 10.8 19.7 20.2 31.4 0 0 0 0 0 0 0 1 0 0 0 20
21 15.1 16 13.2 10.8 19.7 20.2 0 0 0 0 0 0 0 0 1 0 0 21
22 15.6 11 15.1 13.2 10.8 19.7 0 0 0 0 0 0 0 0 0 1 0 22
23 15.5 15 15.6 15.1 13.2 10.8 0 0 0 0 0 0 0 0 0 0 1 23
24 12.7 13 15.5 15.6 15.1 13.2 0 0 0 0 0 0 0 0 0 0 0 24
25 10.9 17 12.7 15.5 15.6 15.1 1 0 0 0 0 0 0 0 0 0 0 25
26 10.0 16 10.9 12.7 15.5 15.6 0 1 0 0 0 0 0 0 0 0 0 26
27 9.1 9 10.0 10.9 12.7 15.5 0 0 1 0 0 0 0 0 0 0 0 27
28 10.3 17 9.1 10.0 10.9 12.7 0 0 0 1 0 0 0 0 0 0 0 28
29 16.9 15 10.3 9.1 10.0 10.9 0 0 0 0 1 0 0 0 0 0 0 29
30 22.0 12 16.9 10.3 9.1 10.0 0 0 0 0 0 1 0 0 0 0 0 30
31 27.6 12 22.0 16.9 10.3 9.1 0 0 0 0 0 0 1 0 0 0 0 31
32 28.9 12 27.6 22.0 16.9 10.3 0 0 0 0 0 0 0 1 0 0 0 32
33 31.0 12 28.9 27.6 22.0 16.9 0 0 0 0 0 0 0 0 1 0 0 33
34 32.9 4 31.0 28.9 27.6 22.0 0 0 0 0 0 0 0 0 0 1 0 34
35 38.1 7 32.9 31.0 28.9 27.6 0 0 0 0 0 0 0 0 0 0 1 35
36 28.8 4 38.1 32.9 31.0 28.9 0 0 0 0 0 0 0 0 0 0 0 36
37 29.0 3 28.8 38.1 32.9 31.0 1 0 0 0 0 0 0 0 0 0 0 37
38 21.8 3 29.0 28.8 38.1 32.9 0 1 0 0 0 0 0 0 0 0 0 38
39 28.8 0 21.8 29.0 28.8 38.1 0 0 1 0 0 0 0 0 0 0 0 39
40 25.6 5 28.8 21.8 29.0 28.8 0 0 0 1 0 0 0 0 0 0 0 40
41 28.2 3 25.6 28.8 21.8 29.0 0 0 0 0 1 0 0 0 0 0 0 41
42 20.2 4 28.2 25.6 28.8 21.8 0 0 0 0 0 1 0 0 0 0 0 42
43 17.9 3 20.2 28.2 25.6 28.8 0 0 0 0 0 0 1 0 0 0 0 43
44 16.3 10 17.9 20.2 28.2 25.6 0 0 0 0 0 0 0 1 0 0 0 44
45 13.2 4 16.3 17.9 20.2 28.2 0 0 0 0 0 0 0 0 1 0 0 45
46 8.1 1 13.2 16.3 17.9 20.2 0 0 0 0 0 0 0 0 0 1 0 46
47 4.5 1 8.1 13.2 16.3 17.9 0 0 0 0 0 0 0 0 0 0 1 47
48 -0.1 8 4.5 8.1 13.2 16.3 0 0 0 0 0 0 0 0 0 0 0 48
49 0.0 5 -0.1 4.5 8.1 13.2 1 0 0 0 0 0 0 0 0 0 0 49
50 2.3 4 0.0 -0.1 4.5 8.1 0 1 0 0 0 0 0 0 0 0 0 50
51 2.8 0 2.3 0.0 -0.1 4.5 0 0 1 0 0 0 0 0 0 0 0 51
52 2.9 2 2.8 2.3 0.0 -0.1 0 0 0 1 0 0 0 0 0 0 0 52
53 0.1 7 2.9 2.8 2.3 0.0 0 0 0 0 1 0 0 0 0 0 0 53
54 3.5 6 0.1 2.9 2.8 2.3 0 0 0 0 0 1 0 0 0 0 0 54
55 8.6 9 3.5 0.1 2.9 2.8 0 0 0 0 0 0 1 0 0 0 0 55
56 13.8 10 8.6 3.5 0.1 2.9 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 Y1 Y2 Y3 Y4
-3.14673 0.15927 1.19239 0.05475 -0.84696 0.51658
M1 M2 M3 M4 M5 M6
1.18489 0.78770 1.77624 1.99845 2.44383 2.51934
M7 M8 M9 M10 M11 t
1.65025 1.96425 1.97879 1.43407 1.18188 0.03609
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-7.0775 -1.8773 -0.4816 2.6805 9.3601
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.14673 5.71413 -0.551 0.585067
X 0.15927 0.20409 0.780 0.439991
Y1 1.19239 0.13805 8.637 1.70e-10 ***
Y2 0.05475 0.19001 0.288 0.774812
Y3 -0.84696 0.18777 -4.511 6.04e-05 ***
Y4 0.51658 0.14344 3.601 0.000903 ***
M1 1.18489 3.04402 0.389 0.699262
M2 0.78770 3.05494 0.258 0.797917
M3 1.77624 3.13858 0.566 0.574764
M4 1.99845 3.07205 0.651 0.519266
M5 2.44383 3.04758 0.802 0.427602
M6 2.51934 3.05860 0.824 0.415255
M7 1.65025 3.06921 0.538 0.593934
M8 1.96425 3.06992 0.640 0.526117
M9 1.97879 3.21532 0.615 0.541943
M10 1.43407 3.37683 0.425 0.673464
M11 1.18188 3.20688 0.369 0.714513
t 0.03609 0.07551 0.478 0.635398
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.522 on 38 degrees of freedom
Multiple R-squared: 0.8589, Adjusted R-squared: 0.7957
F-statistic: 13.6 on 17 and 38 DF, p-value: 2.898e-11
> 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.4583713 0.91674259 0.541628704
[2,] 0.8445518 0.31089635 0.155448174
[3,] 0.7842150 0.43157000 0.215784999
[4,] 0.9491377 0.10172450 0.050862250
[5,] 0.9390260 0.12194799 0.060973996
[6,] 0.9290574 0.14188519 0.070942597
[7,] 0.9249936 0.15001278 0.075006389
[8,] 0.9915682 0.01686366 0.008431832
[9,] 0.9807057 0.03858856 0.019294282
[10,] 0.9573804 0.08523926 0.042619628
[11,] 0.9342191 0.13156184 0.065780919
[12,] 0.8893077 0.22138465 0.110692326
[13,] 0.9030290 0.19394207 0.096971035
[14,] 0.8960187 0.20796268 0.103981338
[15,] 0.8748256 0.25034873 0.125174363
> postscript(file="/var/www/html/rcomp/tmp/1xfqf1258746591.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/24xx61258746591.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/33bmi1258746591.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/4z3f61258746591.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/5hk5i1258746591.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 7
-7.0774890 -6.2031238 -0.8887942 1.0581585 7.7363503 -0.8177568 -0.7080414
8 9 10 11 12 13 14
-2.1182755 -0.1098245 2.6906565 -6.1843939 9.3601232 4.5214238 2.6771160
15 16 17 18 19 20 21
-0.4721117 2.8407700 -6.4004694 0.4518545 -4.2730776 -1.7969246 2.8809059
22 23 24 25 26 27 28
-4.9906834 0.4183662 -0.4560101 -1.3279339 0.2490589 -1.7087746 -0.9969398
29 30 31 32 33 34 35
4.2261223 1.4594968 2.9312148 1.8945198 2.9972421 6.2132543 6.9792180
36 37 38 39 40 41 42
-5.8945923 4.5726652 1.4271279 5.8916904 -1.3419668 -1.6738790 -3.2217107
43 44 45 46 47 48 49
-1.4589554 2.5116810 -5.7683235 -3.9132274 -1.2131903 -3.0095207 -0.6886661
50 51 52 53 54 55 56
1.8498210 -2.8220100 -1.5600220 -3.8881242 2.1281162 3.5088596 -0.4910008
> postscript(file="/var/www/html/rcomp/tmp/6hihn1258746591.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 -7.0774890 NA
1 -6.2031238 -7.0774890
2 -0.8887942 -6.2031238
3 1.0581585 -0.8887942
4 7.7363503 1.0581585
5 -0.8177568 7.7363503
6 -0.7080414 -0.8177568
7 -2.1182755 -0.7080414
8 -0.1098245 -2.1182755
9 2.6906565 -0.1098245
10 -6.1843939 2.6906565
11 9.3601232 -6.1843939
12 4.5214238 9.3601232
13 2.6771160 4.5214238
14 -0.4721117 2.6771160
15 2.8407700 -0.4721117
16 -6.4004694 2.8407700
17 0.4518545 -6.4004694
18 -4.2730776 0.4518545
19 -1.7969246 -4.2730776
20 2.8809059 -1.7969246
21 -4.9906834 2.8809059
22 0.4183662 -4.9906834
23 -0.4560101 0.4183662
24 -1.3279339 -0.4560101
25 0.2490589 -1.3279339
26 -1.7087746 0.2490589
27 -0.9969398 -1.7087746
28 4.2261223 -0.9969398
29 1.4594968 4.2261223
30 2.9312148 1.4594968
31 1.8945198 2.9312148
32 2.9972421 1.8945198
33 6.2132543 2.9972421
34 6.9792180 6.2132543
35 -5.8945923 6.9792180
36 4.5726652 -5.8945923
37 1.4271279 4.5726652
38 5.8916904 1.4271279
39 -1.3419668 5.8916904
40 -1.6738790 -1.3419668
41 -3.2217107 -1.6738790
42 -1.4589554 -3.2217107
43 2.5116810 -1.4589554
44 -5.7683235 2.5116810
45 -3.9132274 -5.7683235
46 -1.2131903 -3.9132274
47 -3.0095207 -1.2131903
48 -0.6886661 -3.0095207
49 1.8498210 -0.6886661
50 -2.8220100 1.8498210
51 -1.5600220 -2.8220100
52 -3.8881242 -1.5600220
53 2.1281162 -3.8881242
54 3.5088596 2.1281162
55 -0.4910008 3.5088596
56 NA -0.4910008
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -6.2031238 -7.0774890
[2,] -0.8887942 -6.2031238
[3,] 1.0581585 -0.8887942
[4,] 7.7363503 1.0581585
[5,] -0.8177568 7.7363503
[6,] -0.7080414 -0.8177568
[7,] -2.1182755 -0.7080414
[8,] -0.1098245 -2.1182755
[9,] 2.6906565 -0.1098245
[10,] -6.1843939 2.6906565
[11,] 9.3601232 -6.1843939
[12,] 4.5214238 9.3601232
[13,] 2.6771160 4.5214238
[14,] -0.4721117 2.6771160
[15,] 2.8407700 -0.4721117
[16,] -6.4004694 2.8407700
[17,] 0.4518545 -6.4004694
[18,] -4.2730776 0.4518545
[19,] -1.7969246 -4.2730776
[20,] 2.8809059 -1.7969246
[21,] -4.9906834 2.8809059
[22,] 0.4183662 -4.9906834
[23,] -0.4560101 0.4183662
[24,] -1.3279339 -0.4560101
[25,] 0.2490589 -1.3279339
[26,] -1.7087746 0.2490589
[27,] -0.9969398 -1.7087746
[28,] 4.2261223 -0.9969398
[29,] 1.4594968 4.2261223
[30,] 2.9312148 1.4594968
[31,] 1.8945198 2.9312148
[32,] 2.9972421 1.8945198
[33,] 6.2132543 2.9972421
[34,] 6.9792180 6.2132543
[35,] -5.8945923 6.9792180
[36,] 4.5726652 -5.8945923
[37,] 1.4271279 4.5726652
[38,] 5.8916904 1.4271279
[39,] -1.3419668 5.8916904
[40,] -1.6738790 -1.3419668
[41,] -3.2217107 -1.6738790
[42,] -1.4589554 -3.2217107
[43,] 2.5116810 -1.4589554
[44,] -5.7683235 2.5116810
[45,] -3.9132274 -5.7683235
[46,] -1.2131903 -3.9132274
[47,] -3.0095207 -1.2131903
[48,] -0.6886661 -3.0095207
[49,] 1.8498210 -0.6886661
[50,] -2.8220100 1.8498210
[51,] -1.5600220 -2.8220100
[52,] -3.8881242 -1.5600220
[53,] 2.1281162 -3.8881242
[54,] 3.5088596 2.1281162
[55,] -0.4910008 3.5088596
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -6.2031238 -7.0774890
2 -0.8887942 -6.2031238
3 1.0581585 -0.8887942
4 7.7363503 1.0581585
5 -0.8177568 7.7363503
6 -0.7080414 -0.8177568
7 -2.1182755 -0.7080414
8 -0.1098245 -2.1182755
9 2.6906565 -0.1098245
10 -6.1843939 2.6906565
11 9.3601232 -6.1843939
12 4.5214238 9.3601232
13 2.6771160 4.5214238
14 -0.4721117 2.6771160
15 2.8407700 -0.4721117
16 -6.4004694 2.8407700
17 0.4518545 -6.4004694
18 -4.2730776 0.4518545
19 -1.7969246 -4.2730776
20 2.8809059 -1.7969246
21 -4.9906834 2.8809059
22 0.4183662 -4.9906834
23 -0.4560101 0.4183662
24 -1.3279339 -0.4560101
25 0.2490589 -1.3279339
26 -1.7087746 0.2490589
27 -0.9969398 -1.7087746
28 4.2261223 -0.9969398
29 1.4594968 4.2261223
30 2.9312148 1.4594968
31 1.8945198 2.9312148
32 2.9972421 1.8945198
33 6.2132543 2.9972421
34 6.9792180 6.2132543
35 -5.8945923 6.9792180
36 4.5726652 -5.8945923
37 1.4271279 4.5726652
38 5.8916904 1.4271279
39 -1.3419668 5.8916904
40 -1.6738790 -1.3419668
41 -3.2217107 -1.6738790
42 -1.4589554 -3.2217107
43 2.5116810 -1.4589554
44 -5.7683235 2.5116810
45 -3.9132274 -5.7683235
46 -1.2131903 -3.9132274
47 -3.0095207 -1.2131903
48 -0.6886661 -3.0095207
49 1.8498210 -0.6886661
50 -2.8220100 1.8498210
51 -1.5600220 -2.8220100
52 -3.8881242 -1.5600220
53 2.1281162 -3.8881242
54 3.5088596 2.1281162
55 -0.4910008 3.5088596
> 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/7ufs31258746591.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/88fgk1258746591.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/97ny71258746591.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/10v0b61258746592.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/11ustl1258746592.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/12zobp1258746592.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/13iqc11258746592.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/14kuj91258746592.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/15q8cm1258746592.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/16gc061258746592.tab")
+ }
> system("convert tmp/1xfqf1258746591.ps tmp/1xfqf1258746591.png")
> system("convert tmp/24xx61258746591.ps tmp/24xx61258746591.png")
> system("convert tmp/33bmi1258746591.ps tmp/33bmi1258746591.png")
> system("convert tmp/4z3f61258746591.ps tmp/4z3f61258746591.png")
> system("convert tmp/5hk5i1258746591.ps tmp/5hk5i1258746591.png")
> system("convert tmp/6hihn1258746591.ps tmp/6hihn1258746591.png")
> system("convert tmp/7ufs31258746591.ps tmp/7ufs31258746591.png")
> system("convert tmp/88fgk1258746591.ps tmp/88fgk1258746591.png")
> system("convert tmp/97ny71258746591.ps tmp/97ny71258746591.png")
> system("convert tmp/10v0b61258746592.ps tmp/10v0b61258746592.png")
>
>
> proc.time()
user system elapsed
2.363 1.567 2.728