R version 2.6.0 (2007-10-03)
Copyright (C) 2007 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,1,89.5,1,91.3,1,97.6,1,100.7,1,104.6,1,94.7,1,101.8,1,102.5,1,105.3,1,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('prijs/olie','war?
'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('prijs/olie','war?
'),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
prijs/olie war?\r
1 87.0 0
2 96.3 0
3 107.1 0
4 115.2 0
5 106.1 1
6 89.5 1
7 91.3 1
8 97.6 1
9 100.7 1
10 104.6 1
11 94.7 1
12 101.8 1
13 102.5 1
14 105.3 1
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) `war?\r`
101.40 76.17
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-88.071 -47.621 8.314 40.104 112.129
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 101.40 27.31 3.713 0.000462 ***
`war?\r` 76.17 28.27 2.694 0.009210 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 54.63 on 58 degrees of freedom
Multiple R-Squared: 0.1112, Adjusted R-squared: 0.09591
F-statistic: 7.259 on 1 and 58 DF, p-value: 0.00921
> postscript(file="/var/www/html/rcomp/tmp/1wx6b1197026060.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/23z2s1197026060.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/3l5jg1197026060.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/4ym8r1197026060.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/5gfhm1197026060.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
-14.400000 -5.100000 5.700000 13.800000 -71.471429 -88.071429 -86.271429
8 9 10 11 12 13 14
-79.971429 -76.871429 -72.971429 -82.871429 -75.771429 -75.071429 -72.271429
15 16 17 18 19 20 21
-67.271429 -67.771429 -60.271429 -58.771429 -46.271429 -51.671429 -44.471429
22 23 24 25 26 27 28
-30.571429 -31.771429 -13.171429 -27.771429 -39.871429 -25.871429 -20.771429
29 30 31 32 33 34 35
2.428571 2.828571 -7.171429 14.028571 21.928571 40.628571 39.928571
36 37 38 39 40 41 42
27.428571 16.428571 21.728571 41.728571 33.528571 37.628571 62.628571
43 44 45 46 47 48 49
64.628571 63.128571 77.828571 75.428571 40.628571 26.128571 28.028571
50 51 52 53 54 55 56
38.028571 10.928571 25.328571 36.428571 52.728571 52.428571 63.428571
57 58 59 60
82.028571 70.228571 92.728571 112.128571
> postscript(file="/var/www/html/rcomp/tmp/6q1z41197026060.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 -14.400000 NA
1 -5.100000 -14.400000
2 5.700000 -5.100000
3 13.800000 5.700000
4 -71.471429 13.800000
5 -88.071429 -71.471429
6 -86.271429 -88.071429
7 -79.971429 -86.271429
8 -76.871429 -79.971429
9 -72.971429 -76.871429
10 -82.871429 -72.971429
11 -75.771429 -82.871429
12 -75.071429 -75.771429
13 -72.271429 -75.071429
14 -67.271429 -72.271429
15 -67.771429 -67.271429
16 -60.271429 -67.771429
17 -58.771429 -60.271429
18 -46.271429 -58.771429
19 -51.671429 -46.271429
20 -44.471429 -51.671429
21 -30.571429 -44.471429
22 -31.771429 -30.571429
23 -13.171429 -31.771429
24 -27.771429 -13.171429
25 -39.871429 -27.771429
26 -25.871429 -39.871429
27 -20.771429 -25.871429
28 2.428571 -20.771429
29 2.828571 2.428571
30 -7.171429 2.828571
31 14.028571 -7.171429
32 21.928571 14.028571
33 40.628571 21.928571
34 39.928571 40.628571
35 27.428571 39.928571
36 16.428571 27.428571
37 21.728571 16.428571
38 41.728571 21.728571
39 33.528571 41.728571
40 37.628571 33.528571
41 62.628571 37.628571
42 64.628571 62.628571
43 63.128571 64.628571
44 77.828571 63.128571
45 75.428571 77.828571
46 40.628571 75.428571
47 26.128571 40.628571
48 28.028571 26.128571
49 38.028571 28.028571
50 10.928571 38.028571
51 25.328571 10.928571
52 36.428571 25.328571
53 52.728571 36.428571
54 52.428571 52.728571
55 63.428571 52.428571
56 82.028571 63.428571
57 70.228571 82.028571
58 92.728571 70.228571
59 112.128571 92.728571
60 NA 112.128571
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -5.100000 -14.400000
[2,] 5.700000 -5.100000
[3,] 13.800000 5.700000
[4,] -71.471429 13.800000
[5,] -88.071429 -71.471429
[6,] -86.271429 -88.071429
[7,] -79.971429 -86.271429
[8,] -76.871429 -79.971429
[9,] -72.971429 -76.871429
[10,] -82.871429 -72.971429
[11,] -75.771429 -82.871429
[12,] -75.071429 -75.771429
[13,] -72.271429 -75.071429
[14,] -67.271429 -72.271429
[15,] -67.771429 -67.271429
[16,] -60.271429 -67.771429
[17,] -58.771429 -60.271429
[18,] -46.271429 -58.771429
[19,] -51.671429 -46.271429
[20,] -44.471429 -51.671429
[21,] -30.571429 -44.471429
[22,] -31.771429 -30.571429
[23,] -13.171429 -31.771429
[24,] -27.771429 -13.171429
[25,] -39.871429 -27.771429
[26,] -25.871429 -39.871429
[27,] -20.771429 -25.871429
[28,] 2.428571 -20.771429
[29,] 2.828571 2.428571
[30,] -7.171429 2.828571
[31,] 14.028571 -7.171429
[32,] 21.928571 14.028571
[33,] 40.628571 21.928571
[34,] 39.928571 40.628571
[35,] 27.428571 39.928571
[36,] 16.428571 27.428571
[37,] 21.728571 16.428571
[38,] 41.728571 21.728571
[39,] 33.528571 41.728571
[40,] 37.628571 33.528571
[41,] 62.628571 37.628571
[42,] 64.628571 62.628571
[43,] 63.128571 64.628571
[44,] 77.828571 63.128571
[45,] 75.428571 77.828571
[46,] 40.628571 75.428571
[47,] 26.128571 40.628571
[48,] 28.028571 26.128571
[49,] 38.028571 28.028571
[50,] 10.928571 38.028571
[51,] 25.328571 10.928571
[52,] 36.428571 25.328571
[53,] 52.728571 36.428571
[54,] 52.428571 52.728571
[55,] 63.428571 52.428571
[56,] 82.028571 63.428571
[57,] 70.228571 82.028571
[58,] 92.728571 70.228571
[59,] 112.128571 92.728571
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -5.100000 -14.400000
2 5.700000 -5.100000
3 13.800000 5.700000
4 -71.471429 13.800000
5 -88.071429 -71.471429
6 -86.271429 -88.071429
7 -79.971429 -86.271429
8 -76.871429 -79.971429
9 -72.971429 -76.871429
10 -82.871429 -72.971429
11 -75.771429 -82.871429
12 -75.071429 -75.771429
13 -72.271429 -75.071429
14 -67.271429 -72.271429
15 -67.771429 -67.271429
16 -60.271429 -67.771429
17 -58.771429 -60.271429
18 -46.271429 -58.771429
19 -51.671429 -46.271429
20 -44.471429 -51.671429
21 -30.571429 -44.471429
22 -31.771429 -30.571429
23 -13.171429 -31.771429
24 -27.771429 -13.171429
25 -39.871429 -27.771429
26 -25.871429 -39.871429
27 -20.771429 -25.871429
28 2.428571 -20.771429
29 2.828571 2.428571
30 -7.171429 2.828571
31 14.028571 -7.171429
32 21.928571 14.028571
33 40.628571 21.928571
34 39.928571 40.628571
35 27.428571 39.928571
36 16.428571 27.428571
37 21.728571 16.428571
38 41.728571 21.728571
39 33.528571 41.728571
40 37.628571 33.528571
41 62.628571 37.628571
42 64.628571 62.628571
43 63.128571 64.628571
44 77.828571 63.128571
45 75.428571 77.828571
46 40.628571 75.428571
47 26.128571 40.628571
48 28.028571 26.128571
49 38.028571 28.028571
50 10.928571 38.028571
51 25.328571 10.928571
52 36.428571 25.328571
53 52.728571 36.428571
54 52.428571 52.728571
55 63.428571 52.428571
56 82.028571 63.428571
57 70.228571 82.028571
58 92.728571 70.228571
59 112.128571 92.728571
> 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/781c41197026060.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/8go2r1197026060.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/91tz51197026060.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
> 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/10ch4p1197026060.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/111fa81197026060.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/123has1197026061.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/13fbbl1197026061.tab")
>
> system("convert tmp/1wx6b1197026060.ps tmp/1wx6b1197026060.png")
> system("convert tmp/23z2s1197026060.ps tmp/23z2s1197026060.png")
> system("convert tmp/3l5jg1197026060.ps tmp/3l5jg1197026060.png")
> system("convert tmp/4ym8r1197026060.ps tmp/4ym8r1197026060.png")
> system("convert tmp/5gfhm1197026060.ps tmp/5gfhm1197026060.png")
> system("convert tmp/6q1z41197026060.ps tmp/6q1z41197026060.png")
> system("convert tmp/781c41197026060.ps tmp/781c41197026060.png")
> system("convert tmp/8go2r1197026060.ps tmp/8go2r1197026060.png")
> system("convert tmp/91tz51197026060.ps tmp/91tz51197026060.png")
>
>
> proc.time()
user system elapsed
2.238 1.459 2.683