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(0,0,0,0,0,0,5,0,0,0,-3,0,0,0,0,0,0,0,-2,0,0,0,1,0,0,0,2,0,0,0,-2,0,0,0,-2,0,0,0,0,0,0,0,1,0,0,0,-2,0,173,183,2,0,70,118,0,-711,215,-110,-2,-2910,357,-41,1,8382,424,85,-1,-1743,-384,314,1,5212,123,242,0,1131,-138,-34,0,2046,195,164,-2,495,135,-160,2,-8138,19,-118,0,-8774,162,114,1,4445,244,-152,-3,611,242,-214,5,684,-227,223,-2,1554,555,124,0,-10927,-59,-410,0,1333,-18,356,2,-54,1155,-432,-2,-11544,-773,363,1,6842,192,-20,1,3572,66,-10,-3,11239,90,173,-1,963,54,44,2,-6157,-7,-328,1,-12126,-348,273,0,-15,-35,-188,-1,571,-6,1,-2,405,-38,238,1,1293,-89,-237,0,-4488,66,112,-1,899,106,-174,0,-9084,336,-18,-3,-2502,-143,-148,1,-14826,4,-65,1,444,10,-40,0,450,-74,30,1,856,-126,-219,1,-1850,289,103,-1,-5322,-92,-507,1,5734,8,74,-2,4214,700,-54,3,-1405,-212,-302,-2,-5082,197,-76,-2,-1907,859,-280,-2,-5241,-52,67,1,-16176,-80,-45,1,-5170,251,-452,1,-10205),dim=c(4,60),dimnames=list(c('h','e','f','w'),1:60))
> y <- array(NA,dim=c(4,60),dimnames=list(c('h','e','f','w'),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
h e f w
1 0 0 0 0
2 0 0 5 0
3 0 0 -3 0
4 0 0 0 0
5 0 0 -2 0
6 0 0 1 0
7 0 0 2 0
8 0 0 -2 0
9 0 0 -2 0
10 0 0 0 0
11 0 0 1 0
12 0 0 -2 0
13 173 183 2 0
14 70 118 0 -711
15 215 -110 -2 -2910
16 357 -41 1 8382
17 424 85 -1 -1743
18 -384 314 1 5212
19 123 242 0 1131
20 -138 -34 0 2046
21 195 164 -2 495
22 135 -160 2 -8138
23 19 -118 0 -8774
24 162 114 1 4445
25 244 -152 -3 611
26 242 -214 5 684
27 -227 223 -2 1554
28 555 124 0 -10927
29 -59 -410 0 1333
30 -18 356 2 -54
31 1155 -432 -2 -11544
32 -773 363 1 6842
33 192 -20 1 3572
34 66 -10 -3 11239
35 90 173 -1 963
36 54 44 2 -6157
37 -7 -328 1 -12126
38 -348 273 0 -15
39 -35 -188 -1 571
40 -6 1 -2 405
41 -38 238 1 1293
42 -89 -237 0 -4488
43 66 112 -1 899
44 106 -174 0 -9084
45 336 -18 -3 -2502
46 -143 -148 1 -14826
47 4 -65 1 444
48 10 -40 0 450
49 -74 30 1 856
50 -126 -219 1 -1850
51 289 103 -1 -5322
52 -92 -507 1 5734
53 8 74 -2 4214
54 700 -54 3 -1405
55 -212 -302 -2 -5082
56 197 -76 -2 -1907
57 859 -280 -2 -5241
58 -52 67 1 -16176
59 -80 -45 1 -5170
60 251 -452 1 -10205
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) e f w
47.8203 -0.3670 -21.0384 -0.0113
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-589.23 -167.63 -34.20 98.22 776.09
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 47.820313 34.921373 1.369 0.1764
e -0.367047 0.190277 -1.929 0.0588 .
f -21.038424 19.089932 -1.102 0.2751
w -0.011300 0.006886 -1.641 0.1064
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 261.2 on 56 degrees of freedom
Multiple R-Squared: 0.1646, Adjusted R-squared: 0.1199
F-statistic: 3.679 on 3 and 56 DF, p-value: 0.01723
> postscript(file="/var/www/html/rcomp/tmp/1sua01195839724.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/2lpk71195839724.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/355wj1195839724.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/44jpx1195839724.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/5mrqt1195839724.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
-47.820313 57.371807 -110.935585 -47.820313 -89.897161 -26.781889
7 8 9 10 11 12
-5.743465 -89.897161 -89.897161 -47.820313 -26.781889 -89.897161
13 14 15 16 17 18
234.426191 57.457049 51.844967 409.884818 366.644582 -236.634030
19 20 21 22 23 24
176.785305 -175.180354 170.892039 -21.429507 -171.277096 227.289497
25 26 27 28 29 30
84.177455 228.552807 -217.485590 429.219695 -242.246958 106.315180
31 32 33 34 35 36
776.092519 -589.229898 198.240359 78.393363 95.522236 -5.166783
37 38 39 40 41 42
-291.195824 -295.785899 -166.411395 -90.953660 37.186121 -274.524413
43 44 45 46 47 48
48.409158 -108.334710 190.185248 -391.637002 -41.622814 -47.417256
49 50 51 52 53 54
-80.097767 -254.070037 197.809140 -240.081322 -7.117941 679.598065
55 56 57 58 59 60
-470.171468 57.658361 607.106891 -236.976679 -181.719429 -57.002607
> postscript(file="/var/www/html/rcomp/tmp/677ns1195839724.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 -47.820313 NA
1 57.371807 -47.820313
2 -110.935585 57.371807
3 -47.820313 -110.935585
4 -89.897161 -47.820313
5 -26.781889 -89.897161
6 -5.743465 -26.781889
7 -89.897161 -5.743465
8 -89.897161 -89.897161
9 -47.820313 -89.897161
10 -26.781889 -47.820313
11 -89.897161 -26.781889
12 234.426191 -89.897161
13 57.457049 234.426191
14 51.844967 57.457049
15 409.884818 51.844967
16 366.644582 409.884818
17 -236.634030 366.644582
18 176.785305 -236.634030
19 -175.180354 176.785305
20 170.892039 -175.180354
21 -21.429507 170.892039
22 -171.277096 -21.429507
23 227.289497 -171.277096
24 84.177455 227.289497
25 228.552807 84.177455
26 -217.485590 228.552807
27 429.219695 -217.485590
28 -242.246958 429.219695
29 106.315180 -242.246958
30 776.092519 106.315180
31 -589.229898 776.092519
32 198.240359 -589.229898
33 78.393363 198.240359
34 95.522236 78.393363
35 -5.166783 95.522236
36 -291.195824 -5.166783
37 -295.785899 -291.195824
38 -166.411395 -295.785899
39 -90.953660 -166.411395
40 37.186121 -90.953660
41 -274.524413 37.186121
42 48.409158 -274.524413
43 -108.334710 48.409158
44 190.185248 -108.334710
45 -391.637002 190.185248
46 -41.622814 -391.637002
47 -47.417256 -41.622814
48 -80.097767 -47.417256
49 -254.070037 -80.097767
50 197.809140 -254.070037
51 -240.081322 197.809140
52 -7.117941 -240.081322
53 679.598065 -7.117941
54 -470.171468 679.598065
55 57.658361 -470.171468
56 607.106891 57.658361
57 -236.976679 607.106891
58 -181.719429 -236.976679
59 -57.002607 -181.719429
60 NA -57.002607
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 57.371807 -47.820313
[2,] -110.935585 57.371807
[3,] -47.820313 -110.935585
[4,] -89.897161 -47.820313
[5,] -26.781889 -89.897161
[6,] -5.743465 -26.781889
[7,] -89.897161 -5.743465
[8,] -89.897161 -89.897161
[9,] -47.820313 -89.897161
[10,] -26.781889 -47.820313
[11,] -89.897161 -26.781889
[12,] 234.426191 -89.897161
[13,] 57.457049 234.426191
[14,] 51.844967 57.457049
[15,] 409.884818 51.844967
[16,] 366.644582 409.884818
[17,] -236.634030 366.644582
[18,] 176.785305 -236.634030
[19,] -175.180354 176.785305
[20,] 170.892039 -175.180354
[21,] -21.429507 170.892039
[22,] -171.277096 -21.429507
[23,] 227.289497 -171.277096
[24,] 84.177455 227.289497
[25,] 228.552807 84.177455
[26,] -217.485590 228.552807
[27,] 429.219695 -217.485590
[28,] -242.246958 429.219695
[29,] 106.315180 -242.246958
[30,] 776.092519 106.315180
[31,] -589.229898 776.092519
[32,] 198.240359 -589.229898
[33,] 78.393363 198.240359
[34,] 95.522236 78.393363
[35,] -5.166783 95.522236
[36,] -291.195824 -5.166783
[37,] -295.785899 -291.195824
[38,] -166.411395 -295.785899
[39,] -90.953660 -166.411395
[40,] 37.186121 -90.953660
[41,] -274.524413 37.186121
[42,] 48.409158 -274.524413
[43,] -108.334710 48.409158
[44,] 190.185248 -108.334710
[45,] -391.637002 190.185248
[46,] -41.622814 -391.637002
[47,] -47.417256 -41.622814
[48,] -80.097767 -47.417256
[49,] -254.070037 -80.097767
[50,] 197.809140 -254.070037
[51,] -240.081322 197.809140
[52,] -7.117941 -240.081322
[53,] 679.598065 -7.117941
[54,] -470.171468 679.598065
[55,] 57.658361 -470.171468
[56,] 607.106891 57.658361
[57,] -236.976679 607.106891
[58,] -181.719429 -236.976679
[59,] -57.002607 -181.719429
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 57.371807 -47.820313
2 -110.935585 57.371807
3 -47.820313 -110.935585
4 -89.897161 -47.820313
5 -26.781889 -89.897161
6 -5.743465 -26.781889
7 -89.897161 -5.743465
8 -89.897161 -89.897161
9 -47.820313 -89.897161
10 -26.781889 -47.820313
11 -89.897161 -26.781889
12 234.426191 -89.897161
13 57.457049 234.426191
14 51.844967 57.457049
15 409.884818 51.844967
16 366.644582 409.884818
17 -236.634030 366.644582
18 176.785305 -236.634030
19 -175.180354 176.785305
20 170.892039 -175.180354
21 -21.429507 170.892039
22 -171.277096 -21.429507
23 227.289497 -171.277096
24 84.177455 227.289497
25 228.552807 84.177455
26 -217.485590 228.552807
27 429.219695 -217.485590
28 -242.246958 429.219695
29 106.315180 -242.246958
30 776.092519 106.315180
31 -589.229898 776.092519
32 198.240359 -589.229898
33 78.393363 198.240359
34 95.522236 78.393363
35 -5.166783 95.522236
36 -291.195824 -5.166783
37 -295.785899 -291.195824
38 -166.411395 -295.785899
39 -90.953660 -166.411395
40 37.186121 -90.953660
41 -274.524413 37.186121
42 48.409158 -274.524413
43 -108.334710 48.409158
44 190.185248 -108.334710
45 -391.637002 190.185248
46 -41.622814 -391.637002
47 -47.417256 -41.622814
48 -80.097767 -47.417256
49 -254.070037 -80.097767
50 197.809140 -254.070037
51 -240.081322 197.809140
52 -7.117941 -240.081322
53 679.598065 -7.117941
54 -470.171468 679.598065
55 57.658361 -470.171468
56 607.106891 57.658361
57 -236.976679 607.106891
58 -181.719429 -236.976679
59 -57.002607 -181.719429
> 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/7nsps1195839724.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/8mwsr1195839724.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/9xjcj1195839724.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/105i3p1195839724.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/11eg2e1195839724.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/120zq11195839725.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/13e0ml1195839725.tab")
>
> system("convert tmp/1sua01195839724.ps tmp/1sua01195839724.png")
> system("convert tmp/2lpk71195839724.ps tmp/2lpk71195839724.png")
> system("convert tmp/355wj1195839724.ps tmp/355wj1195839724.png")
> system("convert tmp/44jpx1195839724.ps tmp/44jpx1195839724.png")
> system("convert tmp/5mrqt1195839724.ps tmp/5mrqt1195839724.png")
> system("convert tmp/677ns1195839724.ps tmp/677ns1195839724.png")
> system("convert tmp/7nsps1195839724.ps tmp/7nsps1195839724.png")
> system("convert tmp/8mwsr1195839724.ps tmp/8mwsr1195839724.png")
> system("convert tmp/9xjcj1195839724.ps tmp/9xjcj1195839724.png")
>
>
> proc.time()
user system elapsed
2.291 1.468 2.704