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(476,2.9,475,2.6,470,2.7,461,1.8,455,1.3,456,0.9,517,1.3,525,1.3,523,1.3,519,1.3,509,1.1,512,1.4,519,1.2,517,1.7,510,1.8,509,1.5,501,1,507,1.6,569,1.5,580,1.8,578,1.8,565,1.6,547,1.9,555,1.7,562,1.6,561,1.3,555,1.1,544,1.9,537,2.6,543,2.3,594,2.4,611,2.2,613,2,611,2.9,594,2.6,595,2.3,591,2.3,589,2.6,584,3.1,573,2.8,567,2.5,569,2.9,621,3.1,629,3.1,628,3.2,612,2.5,595,2.6,597,2.9,593,2.6,590,2.4,580,1.7,574,2,573,2.2,573,1.9,620,1.6,626,1.6,620,1.2,588,1.2,566,1.5,557,1.6),dim=c(2,60),dimnames=list(c('y','x'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('y','x'),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
y x
1 476 2.9
2 475 2.6
3 470 2.7
4 461 1.8
5 455 1.3
6 456 0.9
7 517 1.3
8 525 1.3
9 523 1.3
10 519 1.3
11 509 1.1
12 512 1.4
13 519 1.2
14 517 1.7
15 510 1.8
16 509 1.5
17 501 1.0
18 507 1.6
19 569 1.5
20 580 1.8
21 578 1.8
22 565 1.6
23 547 1.9
24 555 1.7
25 562 1.6
26 561 1.3
27 555 1.1
28 544 1.9
29 537 2.6
30 543 2.3
31 594 2.4
32 611 2.2
33 613 2.0
34 611 2.9
35 594 2.6
36 595 2.3
37 591 2.3
38 589 2.6
39 584 3.1
40 573 2.8
41 567 2.5
42 569 2.9
43 621 3.1
44 629 3.1
45 628 3.2
46 612 2.5
47 595 2.6
48 597 2.9
49 593 2.6
50 590 2.4
51 580 1.7
52 574 2.0
53 573 2.2
54 573 1.9
55 620 1.6
56 626 1.6
57 620 1.2
58 588 1.2
59 566 1.5
60 557 1.6
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x
498.34 30.09
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-109.60 -20.95 11.91 25.43 85.55
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 498.337 18.395 27.091 < 2e-16 ***
x 30.090 8.817 3.413 0.00118 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 43.18 on 58 degrees of freedom
Multiple R-Squared: 0.1672, Adjusted R-squared: 0.1529
F-statistic: 11.65 on 1 and 58 DF, p-value: 0.001178
> postscript(file="/var/www/html/rcomp/tmp/13vfe1195321373.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/26qhq1195321373.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/3rrvt1195321373.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/4dfq41195321373.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/5rlcn1195321373.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
-109.598842 -101.571801 -109.580815 -91.499691 -82.454622 -69.418567
7 8 9 10 11 12
-20.454622 -12.454622 -14.454622 -18.454622 -22.436594 -28.463636
13 14 15 16 17 18
-15.445608 -32.490677 -42.499691 -34.472649 -27.427581 -39.481663
19 20 21 22 23 24
25.527351 27.500309 25.500309 18.518337 -8.508705 5.509323
25 26 27 28 29 30
15.518337 23.545378 23.563406 -11.508705 -39.571801 -24.544760
31 32 33 34 35 36
23.446227 46.464254 54.482282 25.401158 17.428199 27.455240
37 38 39 40 41 42
23.455240 12.428199 -7.616870 -9.589828 -6.562787 -16.598842
43 44 45 46 47 48
29.383130 37.383130 33.374117 38.437213 18.428199 11.401158
49 50 51 52 53 54
16.428199 19.446227 30.509323 15.482282 8.464254 17.491295
55 56 57 58 59 60
73.518337 79.518337 85.554392 53.554392 22.527351 10.518337
> postscript(file="/var/www/html/rcomp/tmp/6nogw1195321373.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 -109.598842 NA
1 -101.571801 -109.598842
2 -109.580815 -101.571801
3 -91.499691 -109.580815
4 -82.454622 -91.499691
5 -69.418567 -82.454622
6 -20.454622 -69.418567
7 -12.454622 -20.454622
8 -14.454622 -12.454622
9 -18.454622 -14.454622
10 -22.436594 -18.454622
11 -28.463636 -22.436594
12 -15.445608 -28.463636
13 -32.490677 -15.445608
14 -42.499691 -32.490677
15 -34.472649 -42.499691
16 -27.427581 -34.472649
17 -39.481663 -27.427581
18 25.527351 -39.481663
19 27.500309 25.527351
20 25.500309 27.500309
21 18.518337 25.500309
22 -8.508705 18.518337
23 5.509323 -8.508705
24 15.518337 5.509323
25 23.545378 15.518337
26 23.563406 23.545378
27 -11.508705 23.563406
28 -39.571801 -11.508705
29 -24.544760 -39.571801
30 23.446227 -24.544760
31 46.464254 23.446227
32 54.482282 46.464254
33 25.401158 54.482282
34 17.428199 25.401158
35 27.455240 17.428199
36 23.455240 27.455240
37 12.428199 23.455240
38 -7.616870 12.428199
39 -9.589828 -7.616870
40 -6.562787 -9.589828
41 -16.598842 -6.562787
42 29.383130 -16.598842
43 37.383130 29.383130
44 33.374117 37.383130
45 38.437213 33.374117
46 18.428199 38.437213
47 11.401158 18.428199
48 16.428199 11.401158
49 19.446227 16.428199
50 30.509323 19.446227
51 15.482282 30.509323
52 8.464254 15.482282
53 17.491295 8.464254
54 73.518337 17.491295
55 79.518337 73.518337
56 85.554392 79.518337
57 53.554392 85.554392
58 22.527351 53.554392
59 10.518337 22.527351
60 NA 10.518337
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -101.571801 -109.598842
[2,] -109.580815 -101.571801
[3,] -91.499691 -109.580815
[4,] -82.454622 -91.499691
[5,] -69.418567 -82.454622
[6,] -20.454622 -69.418567
[7,] -12.454622 -20.454622
[8,] -14.454622 -12.454622
[9,] -18.454622 -14.454622
[10,] -22.436594 -18.454622
[11,] -28.463636 -22.436594
[12,] -15.445608 -28.463636
[13,] -32.490677 -15.445608
[14,] -42.499691 -32.490677
[15,] -34.472649 -42.499691
[16,] -27.427581 -34.472649
[17,] -39.481663 -27.427581
[18,] 25.527351 -39.481663
[19,] 27.500309 25.527351
[20,] 25.500309 27.500309
[21,] 18.518337 25.500309
[22,] -8.508705 18.518337
[23,] 5.509323 -8.508705
[24,] 15.518337 5.509323
[25,] 23.545378 15.518337
[26,] 23.563406 23.545378
[27,] -11.508705 23.563406
[28,] -39.571801 -11.508705
[29,] -24.544760 -39.571801
[30,] 23.446227 -24.544760
[31,] 46.464254 23.446227
[32,] 54.482282 46.464254
[33,] 25.401158 54.482282
[34,] 17.428199 25.401158
[35,] 27.455240 17.428199
[36,] 23.455240 27.455240
[37,] 12.428199 23.455240
[38,] -7.616870 12.428199
[39,] -9.589828 -7.616870
[40,] -6.562787 -9.589828
[41,] -16.598842 -6.562787
[42,] 29.383130 -16.598842
[43,] 37.383130 29.383130
[44,] 33.374117 37.383130
[45,] 38.437213 33.374117
[46,] 18.428199 38.437213
[47,] 11.401158 18.428199
[48,] 16.428199 11.401158
[49,] 19.446227 16.428199
[50,] 30.509323 19.446227
[51,] 15.482282 30.509323
[52,] 8.464254 15.482282
[53,] 17.491295 8.464254
[54,] 73.518337 17.491295
[55,] 79.518337 73.518337
[56,] 85.554392 79.518337
[57,] 53.554392 85.554392
[58,] 22.527351 53.554392
[59,] 10.518337 22.527351
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -101.571801 -109.598842
2 -109.580815 -101.571801
3 -91.499691 -109.580815
4 -82.454622 -91.499691
5 -69.418567 -82.454622
6 -20.454622 -69.418567
7 -12.454622 -20.454622
8 -14.454622 -12.454622
9 -18.454622 -14.454622
10 -22.436594 -18.454622
11 -28.463636 -22.436594
12 -15.445608 -28.463636
13 -32.490677 -15.445608
14 -42.499691 -32.490677
15 -34.472649 -42.499691
16 -27.427581 -34.472649
17 -39.481663 -27.427581
18 25.527351 -39.481663
19 27.500309 25.527351
20 25.500309 27.500309
21 18.518337 25.500309
22 -8.508705 18.518337
23 5.509323 -8.508705
24 15.518337 5.509323
25 23.545378 15.518337
26 23.563406 23.545378
27 -11.508705 23.563406
28 -39.571801 -11.508705
29 -24.544760 -39.571801
30 23.446227 -24.544760
31 46.464254 23.446227
32 54.482282 46.464254
33 25.401158 54.482282
34 17.428199 25.401158
35 27.455240 17.428199
36 23.455240 27.455240
37 12.428199 23.455240
38 -7.616870 12.428199
39 -9.589828 -7.616870
40 -6.562787 -9.589828
41 -16.598842 -6.562787
42 29.383130 -16.598842
43 37.383130 29.383130
44 33.374117 37.383130
45 38.437213 33.374117
46 18.428199 38.437213
47 11.401158 18.428199
48 16.428199 11.401158
49 19.446227 16.428199
50 30.509323 19.446227
51 15.482282 30.509323
52 8.464254 15.482282
53 17.491295 8.464254
54 73.518337 17.491295
55 79.518337 73.518337
56 85.554392 79.518337
57 53.554392 85.554392
58 22.527351 53.554392
59 10.518337 22.527351
> 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/76opw1195321373.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/8vupu1195321373.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/9viv51195321373.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/10punr1195321373.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/115aex1195321373.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/12t3ox1195321373.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/13yajb1195321374.tab")
>
> system("convert tmp/13vfe1195321373.ps tmp/13vfe1195321373.png")
> system("convert tmp/26qhq1195321373.ps tmp/26qhq1195321373.png")
> system("convert tmp/3rrvt1195321373.ps tmp/3rrvt1195321373.png")
> system("convert tmp/4dfq41195321373.ps tmp/4dfq41195321373.png")
> system("convert tmp/5rlcn1195321373.ps tmp/5rlcn1195321373.png")
> system("convert tmp/6nogw1195321373.ps tmp/6nogw1195321373.png")
> system("convert tmp/76opw1195321373.ps tmp/76opw1195321373.png")
> system("convert tmp/8vupu1195321373.ps tmp/8vupu1195321373.png")
> system("convert tmp/9viv51195321373.ps tmp/9viv51195321373.png")
>
>
> proc.time()
user system elapsed
2.212 1.424 2.626