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(5.014,0,6.153,0,6.441,0,5.584,0,6.427,0,6.062,0,5.589,0,6.216,0,5.809,0,4.989,0,6.706,0,7.174,0,6.122,0,8.075,0,6.292,0,6.337,0,8.576,0,6.077,0,5.931,0,6.288,0,7.167,0,6.054,0,6.468,0,6.401,0,6.927,0,7.914,0,7.728,0,8.699,0,8.522,0,6.481,0,7.502,0,7.778,0,7.424,0,6.941,0,8.574,0,9.169,0,7.701,0,9.035,0,7.158,0,8.195,0,8.124,1,7.073,1,7.017,1,7.390,1,7.776,1,6.197,1,6.889,1,7.087,1,6.485,1,7.654,1,6.501,1,6.313,1,7.826,1,6.589,1,6.729,1,5.684,1,8.105,1,6.391,1,5.901,1,6.758,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 5.014 0
2 6.153 0
3 6.441 0
4 5.584 0
5 6.427 0
6 6.062 0
7 5.589 0
8 6.216 0
9 5.809 0
10 4.989 0
11 6.706 0
12 7.174 0
13 6.122 0
14 8.075 0
15 6.292 0
16 6.337 0
17 8.576 0
18 6.077 0
19 5.931 0
20 6.288 0
21 7.167 0
22 6.054 0
23 6.468 0
24 6.401 0
25 6.927 0
26 7.914 0
27 7.728 0
28 8.699 0
29 8.522 0
30 6.481 0
31 7.502 0
32 7.778 0
33 7.424 0
34 6.941 0
35 8.574 0
36 9.169 0
37 7.701 0
38 9.035 0
39 7.158 0
40 8.195 0
41 8.124 1
42 7.073 1
43 7.017 1
44 7.390 1
45 7.776 1
46 6.197 1
47 6.889 1
48 7.087 1
49 6.485 1
50 7.654 1
51 6.501 1
52 6.313 1
53 7.826 1
54 6.589 1
55 6.729 1
56 5.684 1
57 8.105 1
58 6.391 1
59 5.901 1
60 6.758 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x
6.94250 -0.01805
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.9535 -0.6725 -0.1809 0.7652 2.2265
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.94250 0.15579 44.562 <2e-16 ***
x -0.01805 0.26984 -0.067 0.947
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9853 on 58 degrees of freedom
Multiple R-squared: 7.714e-05, Adjusted R-squared: -0.01716
F-statistic: 0.004474 on 1 and 58 DF, p-value: 0.9469
> postscript(file="/var/www/html/rcomp/tmp/17igb1228995014.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/2a2no1228995014.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/3nq561228995014.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/40j4o1228995014.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/5tvz41228995014.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 7 8
-1.92850 -0.78950 -0.50150 -1.35850 -0.51550 -0.88050 -1.35350 -0.72650
9 10 11 12 13 14 15 16
-1.13350 -1.95350 -0.23650 0.23150 -0.82050 1.13250 -0.65050 -0.60550
17 18 19 20 21 22 23 24
1.63350 -0.86550 -1.01150 -0.65450 0.22450 -0.88850 -0.47450 -0.54150
25 26 27 28 29 30 31 32
-0.01550 0.97150 0.78550 1.75650 1.57950 -0.46150 0.55950 0.83550
33 34 35 36 37 38 39 40
0.48150 -0.00150 1.63150 2.22650 0.75850 2.09250 0.21550 1.25250
41 42 43 44 45 46 47 48
1.19955 0.14855 0.09255 0.46555 0.85155 -0.72745 -0.03545 0.16255
49 50 51 52 53 54 55 56
-0.43945 0.72955 -0.42345 -0.61145 0.90155 -0.33545 -0.19545 -1.24045
57 58 59 60
1.18055 -0.53345 -1.02345 -0.16645
> postscript(file="/var/www/html/rcomp/tmp/6to0g1228995014.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 -1.92850 NA
1 -0.78950 -1.92850
2 -0.50150 -0.78950
3 -1.35850 -0.50150
4 -0.51550 -1.35850
5 -0.88050 -0.51550
6 -1.35350 -0.88050
7 -0.72650 -1.35350
8 -1.13350 -0.72650
9 -1.95350 -1.13350
10 -0.23650 -1.95350
11 0.23150 -0.23650
12 -0.82050 0.23150
13 1.13250 -0.82050
14 -0.65050 1.13250
15 -0.60550 -0.65050
16 1.63350 -0.60550
17 -0.86550 1.63350
18 -1.01150 -0.86550
19 -0.65450 -1.01150
20 0.22450 -0.65450
21 -0.88850 0.22450
22 -0.47450 -0.88850
23 -0.54150 -0.47450
24 -0.01550 -0.54150
25 0.97150 -0.01550
26 0.78550 0.97150
27 1.75650 0.78550
28 1.57950 1.75650
29 -0.46150 1.57950
30 0.55950 -0.46150
31 0.83550 0.55950
32 0.48150 0.83550
33 -0.00150 0.48150
34 1.63150 -0.00150
35 2.22650 1.63150
36 0.75850 2.22650
37 2.09250 0.75850
38 0.21550 2.09250
39 1.25250 0.21550
40 1.19955 1.25250
41 0.14855 1.19955
42 0.09255 0.14855
43 0.46555 0.09255
44 0.85155 0.46555
45 -0.72745 0.85155
46 -0.03545 -0.72745
47 0.16255 -0.03545
48 -0.43945 0.16255
49 0.72955 -0.43945
50 -0.42345 0.72955
51 -0.61145 -0.42345
52 0.90155 -0.61145
53 -0.33545 0.90155
54 -0.19545 -0.33545
55 -1.24045 -0.19545
56 1.18055 -1.24045
57 -0.53345 1.18055
58 -1.02345 -0.53345
59 -0.16645 -1.02345
60 NA -0.16645
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.78950 -1.92850
[2,] -0.50150 -0.78950
[3,] -1.35850 -0.50150
[4,] -0.51550 -1.35850
[5,] -0.88050 -0.51550
[6,] -1.35350 -0.88050
[7,] -0.72650 -1.35350
[8,] -1.13350 -0.72650
[9,] -1.95350 -1.13350
[10,] -0.23650 -1.95350
[11,] 0.23150 -0.23650
[12,] -0.82050 0.23150
[13,] 1.13250 -0.82050
[14,] -0.65050 1.13250
[15,] -0.60550 -0.65050
[16,] 1.63350 -0.60550
[17,] -0.86550 1.63350
[18,] -1.01150 -0.86550
[19,] -0.65450 -1.01150
[20,] 0.22450 -0.65450
[21,] -0.88850 0.22450
[22,] -0.47450 -0.88850
[23,] -0.54150 -0.47450
[24,] -0.01550 -0.54150
[25,] 0.97150 -0.01550
[26,] 0.78550 0.97150
[27,] 1.75650 0.78550
[28,] 1.57950 1.75650
[29,] -0.46150 1.57950
[30,] 0.55950 -0.46150
[31,] 0.83550 0.55950
[32,] 0.48150 0.83550
[33,] -0.00150 0.48150
[34,] 1.63150 -0.00150
[35,] 2.22650 1.63150
[36,] 0.75850 2.22650
[37,] 2.09250 0.75850
[38,] 0.21550 2.09250
[39,] 1.25250 0.21550
[40,] 1.19955 1.25250
[41,] 0.14855 1.19955
[42,] 0.09255 0.14855
[43,] 0.46555 0.09255
[44,] 0.85155 0.46555
[45,] -0.72745 0.85155
[46,] -0.03545 -0.72745
[47,] 0.16255 -0.03545
[48,] -0.43945 0.16255
[49,] 0.72955 -0.43945
[50,] -0.42345 0.72955
[51,] -0.61145 -0.42345
[52,] 0.90155 -0.61145
[53,] -0.33545 0.90155
[54,] -0.19545 -0.33545
[55,] -1.24045 -0.19545
[56,] 1.18055 -1.24045
[57,] -0.53345 1.18055
[58,] -1.02345 -0.53345
[59,] -0.16645 -1.02345
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.78950 -1.92850
2 -0.50150 -0.78950
3 -1.35850 -0.50150
4 -0.51550 -1.35850
5 -0.88050 -0.51550
6 -1.35350 -0.88050
7 -0.72650 -1.35350
8 -1.13350 -0.72650
9 -1.95350 -1.13350
10 -0.23650 -1.95350
11 0.23150 -0.23650
12 -0.82050 0.23150
13 1.13250 -0.82050
14 -0.65050 1.13250
15 -0.60550 -0.65050
16 1.63350 -0.60550
17 -0.86550 1.63350
18 -1.01150 -0.86550
19 -0.65450 -1.01150
20 0.22450 -0.65450
21 -0.88850 0.22450
22 -0.47450 -0.88850
23 -0.54150 -0.47450
24 -0.01550 -0.54150
25 0.97150 -0.01550
26 0.78550 0.97150
27 1.75650 0.78550
28 1.57950 1.75650
29 -0.46150 1.57950
30 0.55950 -0.46150
31 0.83550 0.55950
32 0.48150 0.83550
33 -0.00150 0.48150
34 1.63150 -0.00150
35 2.22650 1.63150
36 0.75850 2.22650
37 2.09250 0.75850
38 0.21550 2.09250
39 1.25250 0.21550
40 1.19955 1.25250
41 0.14855 1.19955
42 0.09255 0.14855
43 0.46555 0.09255
44 0.85155 0.46555
45 -0.72745 0.85155
46 -0.03545 -0.72745
47 0.16255 -0.03545
48 -0.43945 0.16255
49 0.72955 -0.43945
50 -0.42345 0.72955
51 -0.61145 -0.42345
52 0.90155 -0.61145
53 -0.33545 0.90155
54 -0.19545 -0.33545
55 -1.24045 -0.19545
56 1.18055 -1.24045
57 -0.53345 1.18055
58 -1.02345 -0.53345
59 -0.16645 -1.02345
> 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/7fcdt1228995014.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/8a5e71228995014.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/9hhpj1228995014.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/100v7p1228995014.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/1125fy1228995014.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/12n47a1228995015.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/13yp8t1228995015.tab")
>
> system("convert tmp/17igb1228995014.ps tmp/17igb1228995014.png")
> system("convert tmp/2a2no1228995014.ps tmp/2a2no1228995014.png")
> system("convert tmp/3nq561228995014.ps tmp/3nq561228995014.png")
> system("convert tmp/40j4o1228995014.ps tmp/40j4o1228995014.png")
> system("convert tmp/5tvz41228995014.ps tmp/5tvz41228995014.png")
> system("convert tmp/6to0g1228995014.ps tmp/6to0g1228995014.png")
> system("convert tmp/7fcdt1228995014.ps tmp/7fcdt1228995014.png")
> system("convert tmp/8a5e71228995014.ps tmp/8a5e71228995014.png")
> system("convert tmp/9hhpj1228995014.ps tmp/9hhpj1228995014.png")
>
>
> proc.time()
user system elapsed
1.901 1.379 2.500