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(10511,0,10812,0,10738,0,10171,0,9721,0,9897,0,9828,0,9924,0,10371,0,10846,0,10413,0,10709,0,10662,0,10570,0,10297,0,10635,0,10872,0,10296,0,10383,0,10431,0,10574,0,10653,0,10805,0,10872,0,10625,0,10407,0,10463,0,10556,0,10646,0,10702,0,11353,1,11346,1,11451,1,11964,1,12574,1,13031,1,13812,1,14544,1,14931,1,14886,1,16005,1,17064,1,15168,1,16050,1,15839,1,15137,1,14954,1,15648,1,15305,1,15579,1,16348,1,15928,1,16171,1,15937,1,15713,1,15594,1,15683,1,16438,1,17032,1,17696,1,17745,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
1 10511 0
2 10812 0
3 10738 0
4 10171 0
5 9721 0
6 9897 0
7 9828 0
8 9924 0
9 10371 0
10 10846 0
11 10413 0
12 10709 0
13 10662 0
14 10570 0
15 10297 0
16 10635 0
17 10872 0
18 10296 0
19 10383 0
20 10431 0
21 10574 0
22 10653 0
23 10805 0
24 10872 0
25 10625 0
26 10407 0
27 10463 0
28 10556 0
29 10646 0
30 10702 0
31 11353 1
32 11346 1
33 11451 1
34 11964 1
35 12574 1
36 13031 1
37 13812 1
38 14544 1
39 14931 1
40 14886 1
41 16005 1
42 17064 1
43 15168 1
44 16050 1
45 15839 1
46 15137 1
47 14954 1
48 15648 1
49 15305 1
50 15579 1
51 16348 1
52 15928 1
53 16171 1
54 15937 1
55 15713 1
56 15594 1
57 15683 1
58 16438 1
59 17032 1
60 17696 1
61 17745 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x
10480 4582
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3716.1 -176.1 145.3 531.9 2682.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10479.7 235.1 44.58 <2e-16 ***
x 4582.5 329.7 13.90 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1288 on 59 degrees of freedom
Multiple R-squared: 0.766, Adjusted R-squared: 0.762
F-statistic: 193.1 on 1 and 59 DF, p-value: < 2.2e-16
> postscript(file="/var/www/html/rcomp/tmp/1alic1227555978.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/2nl651227555978.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/3ly471227555978.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/43fpc1227555978.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/5z6mf1227555978.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
31.33333 332.33333 258.33333 -308.66667 -758.66667 -582.66667
7 8 9 10 11 12
-651.66667 -555.66667 -108.66667 366.33333 -66.66667 229.33333
13 14 15 16 17 18
182.33333 90.33333 -182.66667 155.33333 392.33333 -183.66667
19 20 21 22 23 24
-96.66667 -48.66667 94.33333 173.33333 325.33333 392.33333
25 26 27 28 29 30
145.33333 -72.66667 -16.66667 76.33333 166.33333 222.33333
31 32 33 34 35 36
-3709.12903 -3716.12903 -3611.12903 -3098.12903 -2488.12903 -2031.12903
37 38 39 40 41 42
-1250.12903 -518.12903 -131.12903 -176.12903 942.87097 2001.87097
43 44 45 46 47 48
105.87097 987.87097 776.87097 74.87097 -108.12903 585.87097
49 50 51 52 53 54
242.87097 516.87097 1285.87097 865.87097 1108.87097 874.87097
55 56 57 58 59 60
650.87097 531.87097 620.87097 1375.87097 1969.87097 2633.87097
61
2682.87097
> postscript(file="/var/www/html/rcomp/tmp/6xhly1227555978.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 31.33333 NA
1 332.33333 31.33333
2 258.33333 332.33333
3 -308.66667 258.33333
4 -758.66667 -308.66667
5 -582.66667 -758.66667
6 -651.66667 -582.66667
7 -555.66667 -651.66667
8 -108.66667 -555.66667
9 366.33333 -108.66667
10 -66.66667 366.33333
11 229.33333 -66.66667
12 182.33333 229.33333
13 90.33333 182.33333
14 -182.66667 90.33333
15 155.33333 -182.66667
16 392.33333 155.33333
17 -183.66667 392.33333
18 -96.66667 -183.66667
19 -48.66667 -96.66667
20 94.33333 -48.66667
21 173.33333 94.33333
22 325.33333 173.33333
23 392.33333 325.33333
24 145.33333 392.33333
25 -72.66667 145.33333
26 -16.66667 -72.66667
27 76.33333 -16.66667
28 166.33333 76.33333
29 222.33333 166.33333
30 -3709.12903 222.33333
31 -3716.12903 -3709.12903
32 -3611.12903 -3716.12903
33 -3098.12903 -3611.12903
34 -2488.12903 -3098.12903
35 -2031.12903 -2488.12903
36 -1250.12903 -2031.12903
37 -518.12903 -1250.12903
38 -131.12903 -518.12903
39 -176.12903 -131.12903
40 942.87097 -176.12903
41 2001.87097 942.87097
42 105.87097 2001.87097
43 987.87097 105.87097
44 776.87097 987.87097
45 74.87097 776.87097
46 -108.12903 74.87097
47 585.87097 -108.12903
48 242.87097 585.87097
49 516.87097 242.87097
50 1285.87097 516.87097
51 865.87097 1285.87097
52 1108.87097 865.87097
53 874.87097 1108.87097
54 650.87097 874.87097
55 531.87097 650.87097
56 620.87097 531.87097
57 1375.87097 620.87097
58 1969.87097 1375.87097
59 2633.87097 1969.87097
60 2682.87097 2633.87097
61 NA 2682.87097
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 332.33333 31.33333
[2,] 258.33333 332.33333
[3,] -308.66667 258.33333
[4,] -758.66667 -308.66667
[5,] -582.66667 -758.66667
[6,] -651.66667 -582.66667
[7,] -555.66667 -651.66667
[8,] -108.66667 -555.66667
[9,] 366.33333 -108.66667
[10,] -66.66667 366.33333
[11,] 229.33333 -66.66667
[12,] 182.33333 229.33333
[13,] 90.33333 182.33333
[14,] -182.66667 90.33333
[15,] 155.33333 -182.66667
[16,] 392.33333 155.33333
[17,] -183.66667 392.33333
[18,] -96.66667 -183.66667
[19,] -48.66667 -96.66667
[20,] 94.33333 -48.66667
[21,] 173.33333 94.33333
[22,] 325.33333 173.33333
[23,] 392.33333 325.33333
[24,] 145.33333 392.33333
[25,] -72.66667 145.33333
[26,] -16.66667 -72.66667
[27,] 76.33333 -16.66667
[28,] 166.33333 76.33333
[29,] 222.33333 166.33333
[30,] -3709.12903 222.33333
[31,] -3716.12903 -3709.12903
[32,] -3611.12903 -3716.12903
[33,] -3098.12903 -3611.12903
[34,] -2488.12903 -3098.12903
[35,] -2031.12903 -2488.12903
[36,] -1250.12903 -2031.12903
[37,] -518.12903 -1250.12903
[38,] -131.12903 -518.12903
[39,] -176.12903 -131.12903
[40,] 942.87097 -176.12903
[41,] 2001.87097 942.87097
[42,] 105.87097 2001.87097
[43,] 987.87097 105.87097
[44,] 776.87097 987.87097
[45,] 74.87097 776.87097
[46,] -108.12903 74.87097
[47,] 585.87097 -108.12903
[48,] 242.87097 585.87097
[49,] 516.87097 242.87097
[50,] 1285.87097 516.87097
[51,] 865.87097 1285.87097
[52,] 1108.87097 865.87097
[53,] 874.87097 1108.87097
[54,] 650.87097 874.87097
[55,] 531.87097 650.87097
[56,] 620.87097 531.87097
[57,] 1375.87097 620.87097
[58,] 1969.87097 1375.87097
[59,] 2633.87097 1969.87097
[60,] 2682.87097 2633.87097
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 332.33333 31.33333
2 258.33333 332.33333
3 -308.66667 258.33333
4 -758.66667 -308.66667
5 -582.66667 -758.66667
6 -651.66667 -582.66667
7 -555.66667 -651.66667
8 -108.66667 -555.66667
9 366.33333 -108.66667
10 -66.66667 366.33333
11 229.33333 -66.66667
12 182.33333 229.33333
13 90.33333 182.33333
14 -182.66667 90.33333
15 155.33333 -182.66667
16 392.33333 155.33333
17 -183.66667 392.33333
18 -96.66667 -183.66667
19 -48.66667 -96.66667
20 94.33333 -48.66667
21 173.33333 94.33333
22 325.33333 173.33333
23 392.33333 325.33333
24 145.33333 392.33333
25 -72.66667 145.33333
26 -16.66667 -72.66667
27 76.33333 -16.66667
28 166.33333 76.33333
29 222.33333 166.33333
30 -3709.12903 222.33333
31 -3716.12903 -3709.12903
32 -3611.12903 -3716.12903
33 -3098.12903 -3611.12903
34 -2488.12903 -3098.12903
35 -2031.12903 -2488.12903
36 -1250.12903 -2031.12903
37 -518.12903 -1250.12903
38 -131.12903 -518.12903
39 -176.12903 -131.12903
40 942.87097 -176.12903
41 2001.87097 942.87097
42 105.87097 2001.87097
43 987.87097 105.87097
44 776.87097 987.87097
45 74.87097 776.87097
46 -108.12903 74.87097
47 585.87097 -108.12903
48 242.87097 585.87097
49 516.87097 242.87097
50 1285.87097 516.87097
51 865.87097 1285.87097
52 1108.87097 865.87097
53 874.87097 1108.87097
54 650.87097 874.87097
55 531.87097 650.87097
56 620.87097 531.87097
57 1375.87097 620.87097
58 1969.87097 1375.87097
59 2633.87097 1969.87097
60 2682.87097 2633.87097
> 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/7kgzw1227555978.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/8ag8b1227555978.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/98z3t1227555978.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/106g2w1227555978.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/114rzu1227555978.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/12g79i1227555978.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/1358gn1227555978.tab")
>
> system("convert tmp/1alic1227555978.ps tmp/1alic1227555978.png")
> system("convert tmp/2nl651227555978.ps tmp/2nl651227555978.png")
> system("convert tmp/3ly471227555978.ps tmp/3ly471227555978.png")
> system("convert tmp/43fpc1227555978.ps tmp/43fpc1227555978.png")
> system("convert tmp/5z6mf1227555978.ps tmp/5z6mf1227555978.png")
> system("convert tmp/6xhly1227555978.ps tmp/6xhly1227555978.png")
> system("convert tmp/7kgzw1227555978.ps tmp/7kgzw1227555978.png")
> system("convert tmp/8ag8b1227555978.ps tmp/8ag8b1227555978.png")
> system("convert tmp/98z3t1227555978.ps tmp/98z3t1227555978.png")
>
>
> proc.time()
user system elapsed
1.939 1.435 2.325