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(163414,0,163652,0,164603,0,165257,0,168731,0,171848,0,175032,0,179187,0,187369,0,194147,0,200145,0,203750,0,206464,0,205034,0,211782,0,244562,0,247059,0,255703,0,260218,0,268852,0,279436,0,281514,0,285458,1,288338,1,296369,1,302221,1,311016,1),dim=c(2,27),dimnames=list(c('BBP','ja/nee'),1:27))
> y <- array(NA,dim=c(2,27),dimnames=list(c('BBP','ja/nee'),1:27))
> 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 = '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)
> 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
BBP ja/nee M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 163414 0 1 0 0 0 0 0 0 0 0 0 0
2 163652 0 0 1 0 0 0 0 0 0 0 0 0
3 164603 0 0 0 1 0 0 0 0 0 0 0 0
4 165257 0 0 0 0 1 0 0 0 0 0 0 0
5 168731 0 0 0 0 0 1 0 0 0 0 0 0
6 171848 0 0 0 0 0 0 1 0 0 0 0 0
7 175032 0 0 0 0 0 0 0 1 0 0 0 0
8 179187 0 0 0 0 0 0 0 0 1 0 0 0
9 187369 0 0 0 0 0 0 0 0 0 1 0 0
10 194147 0 0 0 0 0 0 0 0 0 0 1 0
11 200145 0 0 0 0 0 0 0 0 0 0 0 1
12 203750 0 0 0 0 0 0 0 0 0 0 0 0
13 206464 0 1 0 0 0 0 0 0 0 0 0 0
14 205034 0 0 1 0 0 0 0 0 0 0 0 0
15 211782 0 0 0 1 0 0 0 0 0 0 0 0
16 244562 0 0 0 0 1 0 0 0 0 0 0 0
17 247059 0 0 0 0 0 1 0 0 0 0 0 0
18 255703 0 0 0 0 0 0 1 0 0 0 0 0
19 260218 0 0 0 0 0 0 0 1 0 0 0 0
20 268852 0 0 0 0 0 0 0 0 1 0 0 0
21 279436 0 0 0 0 0 0 0 0 0 1 0 0
22 281514 0 0 0 0 0 0 0 0 0 0 1 0
23 285458 1 0 0 0 0 0 0 0 0 0 0 1
24 288338 1 0 0 0 0 0 0 0 0 0 0 0
25 296369 1 1 0 0 0 0 0 0 0 0 0 0
26 302221 1 0 1 0 0 0 0 0 0 0 0 0
27 311016 1 0 0 1 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `ja/nee` M1 M2 M3 M4
192760 106568 -6200 -4647 851 12150
M5 M6 M7 M8 M9 M10
15135 21016 24865 31260 40643 45071
M11
-3242
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-46033 -34086 7540 29534 46033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 192760 34804 5.538 7.3e-05 ***
`ja/nee` 106568 26309 4.051 0.00119 **
M1 -6200 41829 -0.148 0.88428
M2 -4647 41829 -0.111 0.91312
M3 851 41829 0.020 0.98405
M4 12150 47430 0.256 0.80155
M5 15135 47430 0.319 0.75436
M6 21016 47430 0.443 0.66447
M7 24865 47430 0.524 0.60830
M8 31260 47430 0.659 0.52054
M9 40643 47430 0.857 0.40593
M10 45071 47430 0.950 0.35810
M11 -3242 45569 -0.071 0.94428
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 45570 on 14 degrees of freedom
Multiple R-Squared: 0.5659, Adjusted R-squared: 0.1938
F-statistic: 1.521 on 12 and 14 DF, p-value: 0.2252
> postscript(file="/var/www/html/rcomp/tmp/1tb2e1198924218.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/26pfo1198924218.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/3j5l41198924218.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/47qat1198924218.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/5wwkx1198924218.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 = 27
Frequency = 1
1 2 3 4 5 6 7
-23145.574 -24460.907 -29007.907 -39652.500 -39164.000 -41927.500 -42593.000
8 9 10 11 12 13 14
-44832.500 -46033.500 -43683.500 10627.639 10990.139 19904.426 16921.093
15 16 17 18 19 20 21
18171.093 39652.500 39164.000 41927.500 42593.000 44832.500 46033.500
22 23 24 25 26 27
43683.500 -10627.639 -10990.139 3241.148 7539.815 10836.815
> postscript(file="/var/www/html/rcomp/tmp/63jth1198924218.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 = 27
Frequency = 1
lag(myerror, k = 1) myerror
0 -23145.574 NA
1 -24460.907 -23145.574
2 -29007.907 -24460.907
3 -39652.500 -29007.907
4 -39164.000 -39652.500
5 -41927.500 -39164.000
6 -42593.000 -41927.500
7 -44832.500 -42593.000
8 -46033.500 -44832.500
9 -43683.500 -46033.500
10 10627.639 -43683.500
11 10990.139 10627.639
12 19904.426 10990.139
13 16921.093 19904.426
14 18171.093 16921.093
15 39652.500 18171.093
16 39164.000 39652.500
17 41927.500 39164.000
18 42593.000 41927.500
19 44832.500 42593.000
20 46033.500 44832.500
21 43683.500 46033.500
22 -10627.639 43683.500
23 -10990.139 -10627.639
24 3241.148 -10990.139
25 7539.815 3241.148
26 10836.815 7539.815
27 NA 10836.815
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -24460.907 -23145.574
[2,] -29007.907 -24460.907
[3,] -39652.500 -29007.907
[4,] -39164.000 -39652.500
[5,] -41927.500 -39164.000
[6,] -42593.000 -41927.500
[7,] -44832.500 -42593.000
[8,] -46033.500 -44832.500
[9,] -43683.500 -46033.500
[10,] 10627.639 -43683.500
[11,] 10990.139 10627.639
[12,] 19904.426 10990.139
[13,] 16921.093 19904.426
[14,] 18171.093 16921.093
[15,] 39652.500 18171.093
[16,] 39164.000 39652.500
[17,] 41927.500 39164.000
[18,] 42593.000 41927.500
[19,] 44832.500 42593.000
[20,] 46033.500 44832.500
[21,] 43683.500 46033.500
[22,] -10627.639 43683.500
[23,] -10990.139 -10627.639
[24,] 3241.148 -10990.139
[25,] 7539.815 3241.148
[26,] 10836.815 7539.815
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -24460.907 -23145.574
2 -29007.907 -24460.907
3 -39652.500 -29007.907
4 -39164.000 -39652.500
5 -41927.500 -39164.000
6 -42593.000 -41927.500
7 -44832.500 -42593.000
8 -46033.500 -44832.500
9 -43683.500 -46033.500
10 10627.639 -43683.500
11 10990.139 10627.639
12 19904.426 10990.139
13 16921.093 19904.426
14 18171.093 16921.093
15 39652.500 18171.093
16 39164.000 39652.500
17 41927.500 39164.000
18 42593.000 41927.500
19 44832.500 42593.000
20 46033.500 44832.500
21 43683.500 46033.500
22 -10627.639 43683.500
23 -10990.139 -10627.639
24 3241.148 -10990.139
25 7539.815 3241.148
26 10836.815 7539.815
> 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/7rc251198924218.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/8e8i91198924218.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/9vx261198924218.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/10tj291198924218.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/11vzpl1198924218.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/12nggr1198924218.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/13ds9i1198924218.tab")
>
> system("convert tmp/1tb2e1198924218.ps tmp/1tb2e1198924218.png")
> system("convert tmp/26pfo1198924218.ps tmp/26pfo1198924218.png")
> system("convert tmp/3j5l41198924218.ps tmp/3j5l41198924218.png")
> system("convert tmp/47qat1198924218.ps tmp/47qat1198924218.png")
> system("convert tmp/5wwkx1198924218.ps tmp/5wwkx1198924218.png")
> system("convert tmp/63jth1198924218.ps tmp/63jth1198924218.png")
> system("convert tmp/7rc251198924218.ps tmp/7rc251198924218.png")
> system("convert tmp/8e8i91198924218.ps tmp/8e8i91198924218.png")
> system("convert tmp/9vx261198924218.ps tmp/9vx261198924218.png")
>
>
> proc.time()
user system elapsed
2.146 1.403 2.517