R version 2.6.1 (2007-11-26)
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 = '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 t
1 163414 0 1 0 0 0 0 0 0 0 0 0 0 1
2 163652 0 0 1 0 0 0 0 0 0 0 0 0 2
3 164603 0 0 0 1 0 0 0 0 0 0 0 0 3
4 165257 0 0 0 0 1 0 0 0 0 0 0 0 4
5 168731 0 0 0 0 0 1 0 0 0 0 0 0 5
6 171848 0 0 0 0 0 0 1 0 0 0 0 0 6
7 175032 0 0 0 0 0 0 0 1 0 0 0 0 7
8 179187 0 0 0 0 0 0 0 0 1 0 0 0 8
9 187369 0 0 0 0 0 0 0 0 0 1 0 0 9
10 194147 0 0 0 0 0 0 0 0 0 0 1 0 10
11 200145 0 0 0 0 0 0 0 0 0 0 0 1 11
12 203750 0 0 0 0 0 0 0 0 0 0 0 0 12
13 206464 0 1 0 0 0 0 0 0 0 0 0 0 13
14 205034 0 0 1 0 0 0 0 0 0 0 0 0 14
15 211782 0 0 0 1 0 0 0 0 0 0 0 0 15
16 244562 0 0 0 0 1 0 0 0 0 0 0 0 16
17 247059 0 0 0 0 0 1 0 0 0 0 0 0 17
18 255703 0 0 0 0 0 0 1 0 0 0 0 0 18
19 260218 0 0 0 0 0 0 0 1 0 0 0 0 19
20 268852 0 0 0 0 0 0 0 0 1 0 0 0 20
21 279436 0 0 0 0 0 0 0 0 0 1 0 0 21
22 281514 0 0 0 0 0 0 0 0 0 0 1 0 22
23 285458 1 0 0 0 0 0 0 0 0 0 0 1 23
24 288338 1 0 0 0 0 0 0 0 0 0 0 0 24
25 296369 1 1 0 0 0 0 0 0 0 0 0 0 25
26 302221 1 0 1 0 0 0 0 0 0 0 0 0 26
27 311016 1 0 0 1 0 0 0 0 0 0 0 0 27
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `ja/nee` M1 M2 M3 M4
132356 9923 7894 3407 2864 12150
M5 M6 M7 M8 M9 M10
9095 8935 6744 7098 10441 8829
M11 t
2798 6040
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-15294 -6018 -514 6018 17123
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 132356.5 10312.3 12.835 9.29e-09 ***
`ja/nee` 9922.9 10018.4 0.990 0.340
M1 7893.8 11221.7 0.703 0.494
M2 3406.8 11189.1 0.304 0.766
M3 2864.5 11174.2 0.256 0.802
M4 12149.6 12669.3 0.959 0.355
M5 9094.8 12677.1 0.717 0.486
M6 8935.0 12700.7 0.704 0.494
M7 6744.1 12739.8 0.529 0.605
M8 7098.3 12794.4 0.555 0.588
M9 10441.0 12864.3 0.812 0.432
M10 8828.6 12949.1 0.682 0.507
M11 2797.8 12180.4 0.230 0.822
t 6040.3 446.3 13.536 4.87e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12170 on 13 degrees of freedom
Multiple R-Squared: 0.9712, Adjusted R-squared: 0.9425
F-statistic: 33.77 on 13 and 13 DF, p-value: 7.431e-08
> postscript(file="/var/www/html/rcomp/tmp/14d251198924703.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/2rjld1198924703.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/30yz91198924703.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/4t4ta1198924703.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/5uq141198924703.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
17123.3387 15808.0054 11261.0054 -3410.4785 -2921.9785 -5685.4785
7 8 9 10 11 12
-6350.9785 -8590.4785 -9791.4785 -7441.4785 -1453.0349 -1090.5349
13 14 15 16 17 18
-12310.7043 -15294.0376 -14044.0376 3410.4785 2921.9785 5685.4785
19 20 21 22 23 24
6350.9785 8590.4785 9791.4785 7441.4785 1453.0349 1090.5349
25 26 27
-4812.6344 -513.9677 2783.0323
> postscript(file="/var/www/html/rcomp/tmp/6wyem1198924703.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 17123.3387 NA
1 15808.0054 17123.3387
2 11261.0054 15808.0054
3 -3410.4785 11261.0054
4 -2921.9785 -3410.4785
5 -5685.4785 -2921.9785
6 -6350.9785 -5685.4785
7 -8590.4785 -6350.9785
8 -9791.4785 -8590.4785
9 -7441.4785 -9791.4785
10 -1453.0349 -7441.4785
11 -1090.5349 -1453.0349
12 -12310.7043 -1090.5349
13 -15294.0376 -12310.7043
14 -14044.0376 -15294.0376
15 3410.4785 -14044.0376
16 2921.9785 3410.4785
17 5685.4785 2921.9785
18 6350.9785 5685.4785
19 8590.4785 6350.9785
20 9791.4785 8590.4785
21 7441.4785 9791.4785
22 1453.0349 7441.4785
23 1090.5349 1453.0349
24 -4812.6344 1090.5349
25 -513.9677 -4812.6344
26 2783.0323 -513.9677
27 NA 2783.0323
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 15808.0054 17123.3387
[2,] 11261.0054 15808.0054
[3,] -3410.4785 11261.0054
[4,] -2921.9785 -3410.4785
[5,] -5685.4785 -2921.9785
[6,] -6350.9785 -5685.4785
[7,] -8590.4785 -6350.9785
[8,] -9791.4785 -8590.4785
[9,] -7441.4785 -9791.4785
[10,] -1453.0349 -7441.4785
[11,] -1090.5349 -1453.0349
[12,] -12310.7043 -1090.5349
[13,] -15294.0376 -12310.7043
[14,] -14044.0376 -15294.0376
[15,] 3410.4785 -14044.0376
[16,] 2921.9785 3410.4785
[17,] 5685.4785 2921.9785
[18,] 6350.9785 5685.4785
[19,] 8590.4785 6350.9785
[20,] 9791.4785 8590.4785
[21,] 7441.4785 9791.4785
[22,] 1453.0349 7441.4785
[23,] 1090.5349 1453.0349
[24,] -4812.6344 1090.5349
[25,] -513.9677 -4812.6344
[26,] 2783.0323 -513.9677
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 15808.0054 17123.3387
2 11261.0054 15808.0054
3 -3410.4785 11261.0054
4 -2921.9785 -3410.4785
5 -5685.4785 -2921.9785
6 -6350.9785 -5685.4785
7 -8590.4785 -6350.9785
8 -9791.4785 -8590.4785
9 -7441.4785 -9791.4785
10 -1453.0349 -7441.4785
11 -1090.5349 -1453.0349
12 -12310.7043 -1090.5349
13 -15294.0376 -12310.7043
14 -14044.0376 -15294.0376
15 3410.4785 -14044.0376
16 2921.9785 3410.4785
17 5685.4785 2921.9785
18 6350.9785 5685.4785
19 8590.4785 6350.9785
20 9791.4785 8590.4785
21 7441.4785 9791.4785
22 1453.0349 7441.4785
23 1090.5349 1453.0349
24 -4812.6344 1090.5349
25 -513.9677 -4812.6344
26 2783.0323 -513.9677
> 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/78dms1198924704.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/8ddyb1198924704.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/9znty1198924704.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/10jqdn1198924704.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/11qwu31198924704.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/12akhc1198924704.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/132v481198924704.tab")
>
> system("convert tmp/14d251198924703.ps tmp/14d251198924703.png")
> system("convert tmp/2rjld1198924703.ps tmp/2rjld1198924703.png")
> system("convert tmp/30yz91198924703.ps tmp/30yz91198924703.png")
> system("convert tmp/4t4ta1198924703.ps tmp/4t4ta1198924703.png")
> system("convert tmp/5uq141198924703.ps tmp/5uq141198924703.png")
> system("convert tmp/6wyem1198924703.ps tmp/6wyem1198924703.png")
> system("convert tmp/78dms1198924704.ps tmp/78dms1198924704.png")
> system("convert tmp/8ddyb1198924704.ps tmp/8ddyb1198924704.png")
> system("convert tmp/9znty1198924704.ps tmp/9znty1198924704.png")
>
>
> proc.time()
user system elapsed
2.310 1.526 13.903