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(15859.4,0,15258.9,0,15498.6,0,15106.5,0,15023.6,0,12083.0,0,15761.3,0,16942.6,0,15070.3,0,13659.6,0,14768.9,0,14725.1,0,15998.1,0,15370.6,0,14956.9,0,15469.7,0,15101.8,0,11703.7,0,16283.6,0,16726.5,0,14968.9,0,14861.0,0,14583.3,0,15305.8,0,17903.9,0,16379.4,0,15420.3,0,17870.5,0,15912.8,0,13866.5,0,17823.2,0,17872.0,0,17422.0,0,16704.5,0,15991.2,0,16583.6,0,19123.5,0,17838.7,0,17209.4,0,18586.5,0,16258.1,0,15141.6,1,19202.1,1,17746.5,1,19090.1,1,18040.3,1,17515.5,1,17751.8,1,21072.4,1,17170.0,1,19439.5,1,19795.4,1,17574.9,1,16165.4,1,19464.6,1,19932.1,1,19961.2,1,17343.4,1,18924.2,1,18574.1,1,21350.6,1),dim=c(2,61),dimnames=list(c('y','x
'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('y','x
'),1:61))
> 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\r
1 15859.4 0
2 15258.9 0
3 15498.6 0
4 15106.5 0
5 15023.6 0
6 12083.0 0
7 15761.3 0
8 16942.6 0
9 15070.3 0
10 13659.6 0
11 14768.9 0
12 14725.1 0
13 15998.1 0
14 15370.6 0
15 14956.9 0
16 15469.7 0
17 15101.8 0
18 11703.7 0
19 16283.6 0
20 16726.5 0
21 14968.9 0
22 14861.0 0
23 14583.3 0
24 15305.8 0
25 17903.9 0
26 16379.4 0
27 15420.3 0
28 17870.5 0
29 15912.8 0
30 13866.5 0
31 17823.2 0
32 17872.0 0
33 17422.0 0
34 16704.5 0
35 15991.2 0
36 16583.6 0
37 19123.5 0
38 17838.7 0
39 17209.4 0
40 18586.5 0
41 16258.1 0
42 15141.6 1
43 19202.1 1
44 17746.5 1
45 19090.1 1
46 18040.3 1
47 17515.5 1
48 17751.8 1
49 21072.4 1
50 17170.0 1
51 19439.5 1
52 19795.4 1
53 17574.9 1
54 16165.4 1
55 19464.6 1
56 19932.1 1
57 19961.2 1
58 17343.4 1
59 18924.2 1
60 18574.1 1
61 21350.6 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `x\r`
15850 2713
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4146.393 -881.193 9.307 901.815 3273.407
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15850.1 243.5 65.080 < 2e-16 ***
`x\r` 2712.7 425.3 6.378 3.03e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1559 on 59 degrees of freedom
Multiple R-Squared: 0.4081, Adjusted R-squared: 0.398
F-statistic: 40.68 on 1 and 59 DF, p-value: 3.029e-08
> postscript(file="/var/www/html/rcomp/tmp/1eazd1196331917.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/20g1l1196331917.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/336hy1196331917.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/4mzia1196331917.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/59sjw1196331917.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 = 61
Frequency = 1
1 2 3 4 5 6
9.307317 -591.192683 -351.492683 -743.592683 -826.492683 -3767.092683
7 8 9 10 11 12
-88.792683 1092.507317 -779.792683 -2190.492683 -1081.192683 -1124.992683
13 14 15 16 17 18
148.007317 -479.492683 -893.192683 -380.392683 -748.292683 -4146.392683
19 20 21 22 23 24
433.507317 876.407317 -881.192683 -989.092683 -1266.792683 -544.292683
25 26 27 28 29 30
2053.807317 529.307317 -429.792683 2020.407317 62.707317 -1983.592683
31 32 33 34 35 36
1973.107317 2021.907317 1571.907317 854.407317 141.107317 733.507317
37 38 39 40 41 42
3273.407317 1988.607317 1359.307317 2736.407317 408.007317 -3421.185000
43 44 45 46 47 48
639.315000 -816.285000 527.315000 -522.485000 -1047.285000 -810.985000
49 50 51 52 53 54
2509.615000 -1392.785000 876.715000 1232.615000 -987.885000 -2397.385000
55 56 57 58 59 60
901.815000 1369.315000 1398.415000 -1219.385000 361.415000 11.315000
61
2787.815000
> postscript(file="/var/www/html/rcomp/tmp/6nz4p1196331917.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 9.307317 NA
1 -591.192683 9.307317
2 -351.492683 -591.192683
3 -743.592683 -351.492683
4 -826.492683 -743.592683
5 -3767.092683 -826.492683
6 -88.792683 -3767.092683
7 1092.507317 -88.792683
8 -779.792683 1092.507317
9 -2190.492683 -779.792683
10 -1081.192683 -2190.492683
11 -1124.992683 -1081.192683
12 148.007317 -1124.992683
13 -479.492683 148.007317
14 -893.192683 -479.492683
15 -380.392683 -893.192683
16 -748.292683 -380.392683
17 -4146.392683 -748.292683
18 433.507317 -4146.392683
19 876.407317 433.507317
20 -881.192683 876.407317
21 -989.092683 -881.192683
22 -1266.792683 -989.092683
23 -544.292683 -1266.792683
24 2053.807317 -544.292683
25 529.307317 2053.807317
26 -429.792683 529.307317
27 2020.407317 -429.792683
28 62.707317 2020.407317
29 -1983.592683 62.707317
30 1973.107317 -1983.592683
31 2021.907317 1973.107317
32 1571.907317 2021.907317
33 854.407317 1571.907317
34 141.107317 854.407317
35 733.507317 141.107317
36 3273.407317 733.507317
37 1988.607317 3273.407317
38 1359.307317 1988.607317
39 2736.407317 1359.307317
40 408.007317 2736.407317
41 -3421.185000 408.007317
42 639.315000 -3421.185000
43 -816.285000 639.315000
44 527.315000 -816.285000
45 -522.485000 527.315000
46 -1047.285000 -522.485000
47 -810.985000 -1047.285000
48 2509.615000 -810.985000
49 -1392.785000 2509.615000
50 876.715000 -1392.785000
51 1232.615000 876.715000
52 -987.885000 1232.615000
53 -2397.385000 -987.885000
54 901.815000 -2397.385000
55 1369.315000 901.815000
56 1398.415000 1369.315000
57 -1219.385000 1398.415000
58 361.415000 -1219.385000
59 11.315000 361.415000
60 2787.815000 11.315000
61 NA 2787.815000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -591.19268 9.307317
[2,] -351.49268 -591.192683
[3,] -743.59268 -351.492683
[4,] -826.49268 -743.592683
[5,] -3767.09268 -826.492683
[6,] -88.79268 -3767.092683
[7,] 1092.50732 -88.792683
[8,] -779.79268 1092.507317
[9,] -2190.49268 -779.792683
[10,] -1081.19268 -2190.492683
[11,] -1124.99268 -1081.192683
[12,] 148.00732 -1124.992683
[13,] -479.49268 148.007317
[14,] -893.19268 -479.492683
[15,] -380.39268 -893.192683
[16,] -748.29268 -380.392683
[17,] -4146.39268 -748.292683
[18,] 433.50732 -4146.392683
[19,] 876.40732 433.507317
[20,] -881.19268 876.407317
[21,] -989.09268 -881.192683
[22,] -1266.79268 -989.092683
[23,] -544.29268 -1266.792683
[24,] 2053.80732 -544.292683
[25,] 529.30732 2053.807317
[26,] -429.79268 529.307317
[27,] 2020.40732 -429.792683
[28,] 62.70732 2020.407317
[29,] -1983.59268 62.707317
[30,] 1973.10732 -1983.592683
[31,] 2021.90732 1973.107317
[32,] 1571.90732 2021.907317
[33,] 854.40732 1571.907317
[34,] 141.10732 854.407317
[35,] 733.50732 141.107317
[36,] 3273.40732 733.507317
[37,] 1988.60732 3273.407317
[38,] 1359.30732 1988.607317
[39,] 2736.40732 1359.307317
[40,] 408.00732 2736.407317
[41,] -3421.18500 408.007317
[42,] 639.31500 -3421.185000
[43,] -816.28500 639.315000
[44,] 527.31500 -816.285000
[45,] -522.48500 527.315000
[46,] -1047.28500 -522.485000
[47,] -810.98500 -1047.285000
[48,] 2509.61500 -810.985000
[49,] -1392.78500 2509.615000
[50,] 876.71500 -1392.785000
[51,] 1232.61500 876.715000
[52,] -987.88500 1232.615000
[53,] -2397.38500 -987.885000
[54,] 901.81500 -2397.385000
[55,] 1369.31500 901.815000
[56,] 1398.41500 1369.315000
[57,] -1219.38500 1398.415000
[58,] 361.41500 -1219.385000
[59,] 11.31500 361.415000
[60,] 2787.81500 11.315000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -591.19268 9.307317
2 -351.49268 -591.192683
3 -743.59268 -351.492683
4 -826.49268 -743.592683
5 -3767.09268 -826.492683
6 -88.79268 -3767.092683
7 1092.50732 -88.792683
8 -779.79268 1092.507317
9 -2190.49268 -779.792683
10 -1081.19268 -2190.492683
11 -1124.99268 -1081.192683
12 148.00732 -1124.992683
13 -479.49268 148.007317
14 -893.19268 -479.492683
15 -380.39268 -893.192683
16 -748.29268 -380.392683
17 -4146.39268 -748.292683
18 433.50732 -4146.392683
19 876.40732 433.507317
20 -881.19268 876.407317
21 -989.09268 -881.192683
22 -1266.79268 -989.092683
23 -544.29268 -1266.792683
24 2053.80732 -544.292683
25 529.30732 2053.807317
26 -429.79268 529.307317
27 2020.40732 -429.792683
28 62.70732 2020.407317
29 -1983.59268 62.707317
30 1973.10732 -1983.592683
31 2021.90732 1973.107317
32 1571.90732 2021.907317
33 854.40732 1571.907317
34 141.10732 854.407317
35 733.50732 141.107317
36 3273.40732 733.507317
37 1988.60732 3273.407317
38 1359.30732 1988.607317
39 2736.40732 1359.307317
40 408.00732 2736.407317
41 -3421.18500 408.007317
42 639.31500 -3421.185000
43 -816.28500 639.315000
44 527.31500 -816.285000
45 -522.48500 527.315000
46 -1047.28500 -522.485000
47 -810.98500 -1047.285000
48 2509.61500 -810.985000
49 -1392.78500 2509.615000
50 876.71500 -1392.785000
51 1232.61500 876.715000
52 -987.88500 1232.615000
53 -2397.38500 -987.885000
54 901.81500 -2397.385000
55 1369.31500 901.815000
56 1398.41500 1369.315000
57 -1219.38500 1398.415000
58 361.41500 -1219.385000
59 11.31500 361.415000
60 2787.81500 11.315000
> 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/73zsz1196331917.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/8helv1196331917.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/9glir1196331917.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/10vvok1196331917.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/11a2rn1196331917.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/124sdh1196331917.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/13t5b11196331917.tab")
>
> system("convert tmp/1eazd1196331917.ps tmp/1eazd1196331917.png")
> system("convert tmp/20g1l1196331917.ps tmp/20g1l1196331917.png")
> system("convert tmp/336hy1196331917.ps tmp/336hy1196331917.png")
> system("convert tmp/4mzia1196331917.ps tmp/4mzia1196331917.png")
> system("convert tmp/59sjw1196331917.ps tmp/59sjw1196331917.png")
> system("convert tmp/6nz4p1196331917.ps tmp/6nz4p1196331917.png")
> system("convert tmp/73zsz1196331917.ps tmp/73zsz1196331917.png")
> system("convert tmp/8helv1196331917.ps tmp/8helv1196331917.png")
> system("convert tmp/9glir1196331917.ps tmp/9glir1196331917.png")
>
>
> proc.time()
user system elapsed
2.247 1.474 2.656