R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(87.0,0,96.3,0,107.1,0,115.2,0,106.1,0,89.5,0,91.3,0,97.6,0,100.7,0,104.6,0,94.7,0,101.8,0,102.5,0,105.3,0,110.3,1,109.8,1,117.3,1,118.8,1,131.3,1,125.9,1,133.1,1,147.0,1,145.8,1,164.4,1,149.8,1,137.7,1,151.7,1,156.8,1,180.0,1,180.4,1,170.4,1,191.6,1,199.5,1,218.2,1,217.5,1,205.0,1,194.0,1,199.3,1,219.3,1,211.1,1,215.2,1,240.2,1,242.2,1,240.7,1,255.4,1,253.0,1,218.2,1,203.7,1,205.6,1,215.6,1,188.5,1,202.9,1,214.0,1,230.3,1,230.0,1,241.0,1,259.6,1,247.8,1,270.3,1,289.7,1),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal 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)
> 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
1 87.0 0
2 96.3 0
3 107.1 0
4 115.2 0
5 106.1 0
6 89.5 0
7 91.3 0
8 97.6 0
9 100.7 0
10 104.6 0
11 94.7 0
12 101.8 0
13 102.5 0
14 105.3 0
15 110.3 1
16 109.8 1
17 117.3 1
18 118.8 1
19 131.3 1
20 125.9 1
21 133.1 1
22 147.0 1
23 145.8 1
24 164.4 1
25 149.8 1
26 137.7 1
27 151.7 1
28 156.8 1
29 180.0 1
30 180.4 1
31 170.4 1
32 191.6 1
33 199.5 1
34 218.2 1
35 217.5 1
36 205.0 1
37 194.0 1
38 199.3 1
39 219.3 1
40 211.1 1
41 215.2 1
42 240.2 1
43 242.2 1
44 240.7 1
45 255.4 1
46 253.0 1
47 218.2 1
48 203.7 1
49 205.6 1
50 215.6 1
51 188.5 1
52 202.9 1
53 214.0 1
54 230.3 1
55 230.0 1
56 241.0 1
57 259.6 1
58 247.8 1
59 270.3 1
60 289.7 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X
99.98 94.58
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-84.763 -16.963 4.679 23.112 95.137
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 99.98 11.03 9.068 1.02e-12 ***
X 94.58 12.59 7.512 4.05e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 41.25 on 58 degrees of freedom
Multiple R-squared: 0.4931, Adjusted R-squared: 0.4844
F-statistic: 56.42 on 1 and 58 DF, p-value: 4.047e-10
> postscript(file="/var/www/html/rcomp/tmp/1dfhs1229455478.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/2zxz41229455478.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/3drwg1229455478.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/4fvgk1229455478.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/54jrr1229455478.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 60
Frequency = 1
1 2 3 4 5 6
-12.9785714 -3.6785714 7.1214286 15.2214286 6.1214286 -10.4785714
7 8 9 10 11 12
-8.6785714 -2.3785714 0.7214286 4.6214286 -5.2785714 1.8214286
13 14 15 16 17 18
2.5214286 5.3214286 -84.2630435 -84.7630435 -77.2630435 -75.7630435
19 20 21 22 23 24
-63.2630435 -68.6630435 -61.4630435 -47.5630435 -48.7630435 -30.1630435
25 26 27 28 29 30
-44.7630435 -56.8630435 -42.8630435 -37.7630435 -14.5630435 -14.1630435
31 32 33 34 35 36
-24.1630435 -2.9630435 4.9369565 23.6369565 22.9369565 10.4369565
37 38 39 40 41 42
-0.5630435 4.7369565 24.7369565 16.5369565 20.6369565 45.6369565
43 44 45 46 47 48
47.6369565 46.1369565 60.8369565 58.4369565 23.6369565 9.1369565
49 50 51 52 53 54
11.0369565 21.0369565 -6.0630435 8.3369565 19.4369565 35.7369565
55 56 57 58 59 60
35.4369565 46.4369565 65.0369565 53.2369565 75.7369565 95.1369565
> postscript(file="/var/www/html/rcomp/tmp/6zz3w1229455478.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -12.9785714 NA
1 -3.6785714 -12.9785714
2 7.1214286 -3.6785714
3 15.2214286 7.1214286
4 6.1214286 15.2214286
5 -10.4785714 6.1214286
6 -8.6785714 -10.4785714
7 -2.3785714 -8.6785714
8 0.7214286 -2.3785714
9 4.6214286 0.7214286
10 -5.2785714 4.6214286
11 1.8214286 -5.2785714
12 2.5214286 1.8214286
13 5.3214286 2.5214286
14 -84.2630435 5.3214286
15 -84.7630435 -84.2630435
16 -77.2630435 -84.7630435
17 -75.7630435 -77.2630435
18 -63.2630435 -75.7630435
19 -68.6630435 -63.2630435
20 -61.4630435 -68.6630435
21 -47.5630435 -61.4630435
22 -48.7630435 -47.5630435
23 -30.1630435 -48.7630435
24 -44.7630435 -30.1630435
25 -56.8630435 -44.7630435
26 -42.8630435 -56.8630435
27 -37.7630435 -42.8630435
28 -14.5630435 -37.7630435
29 -14.1630435 -14.5630435
30 -24.1630435 -14.1630435
31 -2.9630435 -24.1630435
32 4.9369565 -2.9630435
33 23.6369565 4.9369565
34 22.9369565 23.6369565
35 10.4369565 22.9369565
36 -0.5630435 10.4369565
37 4.7369565 -0.5630435
38 24.7369565 4.7369565
39 16.5369565 24.7369565
40 20.6369565 16.5369565
41 45.6369565 20.6369565
42 47.6369565 45.6369565
43 46.1369565 47.6369565
44 60.8369565 46.1369565
45 58.4369565 60.8369565
46 23.6369565 58.4369565
47 9.1369565 23.6369565
48 11.0369565 9.1369565
49 21.0369565 11.0369565
50 -6.0630435 21.0369565
51 8.3369565 -6.0630435
52 19.4369565 8.3369565
53 35.7369565 19.4369565
54 35.4369565 35.7369565
55 46.4369565 35.4369565
56 65.0369565 46.4369565
57 53.2369565 65.0369565
58 75.7369565 53.2369565
59 95.1369565 75.7369565
60 NA 95.1369565
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.6785714 -12.9785714
[2,] 7.1214286 -3.6785714
[3,] 15.2214286 7.1214286
[4,] 6.1214286 15.2214286
[5,] -10.4785714 6.1214286
[6,] -8.6785714 -10.4785714
[7,] -2.3785714 -8.6785714
[8,] 0.7214286 -2.3785714
[9,] 4.6214286 0.7214286
[10,] -5.2785714 4.6214286
[11,] 1.8214286 -5.2785714
[12,] 2.5214286 1.8214286
[13,] 5.3214286 2.5214286
[14,] -84.2630435 5.3214286
[15,] -84.7630435 -84.2630435
[16,] -77.2630435 -84.7630435
[17,] -75.7630435 -77.2630435
[18,] -63.2630435 -75.7630435
[19,] -68.6630435 -63.2630435
[20,] -61.4630435 -68.6630435
[21,] -47.5630435 -61.4630435
[22,] -48.7630435 -47.5630435
[23,] -30.1630435 -48.7630435
[24,] -44.7630435 -30.1630435
[25,] -56.8630435 -44.7630435
[26,] -42.8630435 -56.8630435
[27,] -37.7630435 -42.8630435
[28,] -14.5630435 -37.7630435
[29,] -14.1630435 -14.5630435
[30,] -24.1630435 -14.1630435
[31,] -2.9630435 -24.1630435
[32,] 4.9369565 -2.9630435
[33,] 23.6369565 4.9369565
[34,] 22.9369565 23.6369565
[35,] 10.4369565 22.9369565
[36,] -0.5630435 10.4369565
[37,] 4.7369565 -0.5630435
[38,] 24.7369565 4.7369565
[39,] 16.5369565 24.7369565
[40,] 20.6369565 16.5369565
[41,] 45.6369565 20.6369565
[42,] 47.6369565 45.6369565
[43,] 46.1369565 47.6369565
[44,] 60.8369565 46.1369565
[45,] 58.4369565 60.8369565
[46,] 23.6369565 58.4369565
[47,] 9.1369565 23.6369565
[48,] 11.0369565 9.1369565
[49,] 21.0369565 11.0369565
[50,] -6.0630435 21.0369565
[51,] 8.3369565 -6.0630435
[52,] 19.4369565 8.3369565
[53,] 35.7369565 19.4369565
[54,] 35.4369565 35.7369565
[55,] 46.4369565 35.4369565
[56,] 65.0369565 46.4369565
[57,] 53.2369565 65.0369565
[58,] 75.7369565 53.2369565
[59,] 95.1369565 75.7369565
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.6785714 -12.9785714
2 7.1214286 -3.6785714
3 15.2214286 7.1214286
4 6.1214286 15.2214286
5 -10.4785714 6.1214286
6 -8.6785714 -10.4785714
7 -2.3785714 -8.6785714
8 0.7214286 -2.3785714
9 4.6214286 0.7214286
10 -5.2785714 4.6214286
11 1.8214286 -5.2785714
12 2.5214286 1.8214286
13 5.3214286 2.5214286
14 -84.2630435 5.3214286
15 -84.7630435 -84.2630435
16 -77.2630435 -84.7630435
17 -75.7630435 -77.2630435
18 -63.2630435 -75.7630435
19 -68.6630435 -63.2630435
20 -61.4630435 -68.6630435
21 -47.5630435 -61.4630435
22 -48.7630435 -47.5630435
23 -30.1630435 -48.7630435
24 -44.7630435 -30.1630435
25 -56.8630435 -44.7630435
26 -42.8630435 -56.8630435
27 -37.7630435 -42.8630435
28 -14.5630435 -37.7630435
29 -14.1630435 -14.5630435
30 -24.1630435 -14.1630435
31 -2.9630435 -24.1630435
32 4.9369565 -2.9630435
33 23.6369565 4.9369565
34 22.9369565 23.6369565
35 10.4369565 22.9369565
36 -0.5630435 10.4369565
37 4.7369565 -0.5630435
38 24.7369565 4.7369565
39 16.5369565 24.7369565
40 20.6369565 16.5369565
41 45.6369565 20.6369565
42 47.6369565 45.6369565
43 46.1369565 47.6369565
44 60.8369565 46.1369565
45 58.4369565 60.8369565
46 23.6369565 58.4369565
47 9.1369565 23.6369565
48 11.0369565 9.1369565
49 21.0369565 11.0369565
50 -6.0630435 21.0369565
51 8.3369565 -6.0630435
52 19.4369565 8.3369565
53 35.7369565 19.4369565
54 35.4369565 35.7369565
55 46.4369565 35.4369565
56 65.0369565 46.4369565
57 53.2369565 65.0369565
58 75.7369565 53.2369565
59 95.1369565 75.7369565
> 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/7vhs31229455478.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/82weo1229455478.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/9b79g1229455478.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
>
> #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/10vntf1229455478.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/11rpti1229455478.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/12g74w1229455478.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/13jbvh1229455478.tab")
>
> system("convert tmp/1dfhs1229455478.ps tmp/1dfhs1229455478.png")
> system("convert tmp/2zxz41229455478.ps tmp/2zxz41229455478.png")
> system("convert tmp/3drwg1229455478.ps tmp/3drwg1229455478.png")
> system("convert tmp/4fvgk1229455478.ps tmp/4fvgk1229455478.png")
> system("convert tmp/54jrr1229455478.ps tmp/54jrr1229455478.png")
> system("convert tmp/6zz3w1229455478.ps tmp/6zz3w1229455478.png")
> system("convert tmp/7vhs31229455478.ps tmp/7vhs31229455478.png")
> system("convert tmp/82weo1229455478.ps tmp/82weo1229455478.png")
> system("convert tmp/9b79g1229455478.ps tmp/9b79g1229455478.png")
>
>
> proc.time()
user system elapsed
1.821 1.326 2.360