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(281,0.88,295,0.87,294,0.88,302,0.89,314,0.92,321,0.96,313,0.99,310,0.98,319,0.98,316,0.98,319,1.00,333,1.02,356,1.06,358,1.08,340,1.08,328,1.08,355,1.16,356,1.17,351,1.14,359,1.11,378,1.12,378,1.17,389.,1.17,407,1.23,413,1.26,404,1.26,406,1.23,402,1.20,383,1.20,392,1.21,398,1.23,400,1.22,405,1.22,420,1.25,439,1.30,441,1.34,424,1.31,423,1.30,434,1.32,429,1.29,421,1.27,430,1.22,424,1.20,437,1.23,456,1.23,469,1.20,476,1.18,510,1.19,549,1.21,554,1.19,557,1.20,610,1.23,675,1.28,596,1.27,633,1.27,632,1.28,596,1.27,585,1.26,627,1.29,629,1.32),dim=c(2,60),dimnames=list(c('Goud','Dollar'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Goud','Dollar'),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
Goud Dollar
1 281 0.88
2 295 0.87
3 294 0.88
4 302 0.89
5 314 0.92
6 321 0.96
7 313 0.99
8 310 0.98
9 319 0.98
10 316 0.98
11 319 1.00
12 333 1.02
13 356 1.06
14 358 1.08
15 340 1.08
16 328 1.08
17 355 1.16
18 356 1.17
19 351 1.14
20 359 1.11
21 378 1.12
22 378 1.17
23 389 1.17
24 407 1.23
25 413 1.26
26 404 1.26
27 406 1.23
28 402 1.20
29 383 1.20
30 392 1.21
31 398 1.23
32 400 1.22
33 405 1.22
34 420 1.25
35 439 1.30
36 441 1.34
37 424 1.31
38 423 1.30
39 434 1.32
40 429 1.29
41 421 1.27
42 430 1.22
43 424 1.20
44 437 1.23
45 456 1.23
46 469 1.20
47 476 1.18
48 510 1.19
49 549 1.21
50 554 1.19
51 557 1.20
52 610 1.23
53 675 1.28
54 596 1.27
55 633 1.27
56 632 1.28
57 596 1.27
58 585 1.26
59 627 1.29
60 629 1.32
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dollar
-241.3 573.1
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-85.61 -58.41 -21.00 34.38 182.77
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -241.29 87.74 -2.750 0.00793 **
Dollar 573.06 74.92 7.649 2.38e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 73.62 on 58 degrees of freedom
Multiple R-Squared: 0.5022, Adjusted R-squared: 0.4936
F-statistic: 58.51 on 1 and 58 DF, p-value: 2.379e-10
> postscript(file="/var/www/html/rcomp/tmp/1eyde1200517975.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/2ixuw1200517975.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/3od8b1200517975.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/4nmyi1200517975.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/5sucb1200517975.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 7
17.995123 37.725743 30.995123 33.264503 28.072643 12.150162 -13.041699
8 9 10 11 12 13 14
-10.311079 -1.311079 -4.311079 -12.772319 -10.233559 -10.156040 -19.617280
15 16 17 18 19 20 21
-37.617280 -49.617280 -68.462242 -73.192862 -61.001001 -35.809141 -22.539761
22 23 24 25 26 27 28
-51.192862 -40.192862 -56.576583 -67.768443 -76.768443 -57.576583 -44.384722
29 30 31 32 33 34 35
-63.384722 -60.115342 -65.576583 -57.845963 -52.845963 -55.037823 -64.690924
36 37 38 39 40 41 42
-85.613405 -85.421544 -80.690924 -81.152164 -68.960304 -65.499063 -27.845963
43 44 45 46 47 48 49
-22.384722 -26.576583 -7.576583 22.615278 41.076518 69.345898 96.884658
50 51 52 53 54 55 56
113.345898 110.615278 146.423417 182.770316 109.500937 146.500937 139.770316
57 58 59 60
109.500937 104.231557 129.039696 113.847836
> postscript(file="/var/www/html/rcomp/tmp/6puh41200517975.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 17.995123 NA
1 37.725743 17.995123
2 30.995123 37.725743
3 33.264503 30.995123
4 28.072643 33.264503
5 12.150162 28.072643
6 -13.041699 12.150162
7 -10.311079 -13.041699
8 -1.311079 -10.311079
9 -4.311079 -1.311079
10 -12.772319 -4.311079
11 -10.233559 -12.772319
12 -10.156040 -10.233559
13 -19.617280 -10.156040
14 -37.617280 -19.617280
15 -49.617280 -37.617280
16 -68.462242 -49.617280
17 -73.192862 -68.462242
18 -61.001001 -73.192862
19 -35.809141 -61.001001
20 -22.539761 -35.809141
21 -51.192862 -22.539761
22 -40.192862 -51.192862
23 -56.576583 -40.192862
24 -67.768443 -56.576583
25 -76.768443 -67.768443
26 -57.576583 -76.768443
27 -44.384722 -57.576583
28 -63.384722 -44.384722
29 -60.115342 -63.384722
30 -65.576583 -60.115342
31 -57.845963 -65.576583
32 -52.845963 -57.845963
33 -55.037823 -52.845963
34 -64.690924 -55.037823
35 -85.613405 -64.690924
36 -85.421544 -85.613405
37 -80.690924 -85.421544
38 -81.152164 -80.690924
39 -68.960304 -81.152164
40 -65.499063 -68.960304
41 -27.845963 -65.499063
42 -22.384722 -27.845963
43 -26.576583 -22.384722
44 -7.576583 -26.576583
45 22.615278 -7.576583
46 41.076518 22.615278
47 69.345898 41.076518
48 96.884658 69.345898
49 113.345898 96.884658
50 110.615278 113.345898
51 146.423417 110.615278
52 182.770316 146.423417
53 109.500937 182.770316
54 146.500937 109.500937
55 139.770316 146.500937
56 109.500937 139.770316
57 104.231557 109.500937
58 129.039696 104.231557
59 113.847836 129.039696
60 NA 113.847836
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 37.725743 17.995123
[2,] 30.995123 37.725743
[3,] 33.264503 30.995123
[4,] 28.072643 33.264503
[5,] 12.150162 28.072643
[6,] -13.041699 12.150162
[7,] -10.311079 -13.041699
[8,] -1.311079 -10.311079
[9,] -4.311079 -1.311079
[10,] -12.772319 -4.311079
[11,] -10.233559 -12.772319
[12,] -10.156040 -10.233559
[13,] -19.617280 -10.156040
[14,] -37.617280 -19.617280
[15,] -49.617280 -37.617280
[16,] -68.462242 -49.617280
[17,] -73.192862 -68.462242
[18,] -61.001001 -73.192862
[19,] -35.809141 -61.001001
[20,] -22.539761 -35.809141
[21,] -51.192862 -22.539761
[22,] -40.192862 -51.192862
[23,] -56.576583 -40.192862
[24,] -67.768443 -56.576583
[25,] -76.768443 -67.768443
[26,] -57.576583 -76.768443
[27,] -44.384722 -57.576583
[28,] -63.384722 -44.384722
[29,] -60.115342 -63.384722
[30,] -65.576583 -60.115342
[31,] -57.845963 -65.576583
[32,] -52.845963 -57.845963
[33,] -55.037823 -52.845963
[34,] -64.690924 -55.037823
[35,] -85.613405 -64.690924
[36,] -85.421544 -85.613405
[37,] -80.690924 -85.421544
[38,] -81.152164 -80.690924
[39,] -68.960304 -81.152164
[40,] -65.499063 -68.960304
[41,] -27.845963 -65.499063
[42,] -22.384722 -27.845963
[43,] -26.576583 -22.384722
[44,] -7.576583 -26.576583
[45,] 22.615278 -7.576583
[46,] 41.076518 22.615278
[47,] 69.345898 41.076518
[48,] 96.884658 69.345898
[49,] 113.345898 96.884658
[50,] 110.615278 113.345898
[51,] 146.423417 110.615278
[52,] 182.770316 146.423417
[53,] 109.500937 182.770316
[54,] 146.500937 109.500937
[55,] 139.770316 146.500937
[56,] 109.500937 139.770316
[57,] 104.231557 109.500937
[58,] 129.039696 104.231557
[59,] 113.847836 129.039696
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 37.725743 17.995123
2 30.995123 37.725743
3 33.264503 30.995123
4 28.072643 33.264503
5 12.150162 28.072643
6 -13.041699 12.150162
7 -10.311079 -13.041699
8 -1.311079 -10.311079
9 -4.311079 -1.311079
10 -12.772319 -4.311079
11 -10.233559 -12.772319
12 -10.156040 -10.233559
13 -19.617280 -10.156040
14 -37.617280 -19.617280
15 -49.617280 -37.617280
16 -68.462242 -49.617280
17 -73.192862 -68.462242
18 -61.001001 -73.192862
19 -35.809141 -61.001001
20 -22.539761 -35.809141
21 -51.192862 -22.539761
22 -40.192862 -51.192862
23 -56.576583 -40.192862
24 -67.768443 -56.576583
25 -76.768443 -67.768443
26 -57.576583 -76.768443
27 -44.384722 -57.576583
28 -63.384722 -44.384722
29 -60.115342 -63.384722
30 -65.576583 -60.115342
31 -57.845963 -65.576583
32 -52.845963 -57.845963
33 -55.037823 -52.845963
34 -64.690924 -55.037823
35 -85.613405 -64.690924
36 -85.421544 -85.613405
37 -80.690924 -85.421544
38 -81.152164 -80.690924
39 -68.960304 -81.152164
40 -65.499063 -68.960304
41 -27.845963 -65.499063
42 -22.384722 -27.845963
43 -26.576583 -22.384722
44 -7.576583 -26.576583
45 22.615278 -7.576583
46 41.076518 22.615278
47 69.345898 41.076518
48 96.884658 69.345898
49 113.345898 96.884658
50 110.615278 113.345898
51 146.423417 110.615278
52 182.770316 146.423417
53 109.500937 182.770316
54 146.500937 109.500937
55 139.770316 146.500937
56 109.500937 139.770316
57 104.231557 109.500937
58 129.039696 104.231557
59 113.847836 129.039696
> 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/7109l1200517975.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/86h6y1200517975.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/9jx5k1200517975.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/1025nf1200517975.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/11kfdb1200517975.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/12cisj1200517975.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/13k8wf1200517975.tab")
>
> system("convert tmp/1eyde1200517975.ps tmp/1eyde1200517975.png")
> system("convert tmp/2ixuw1200517975.ps tmp/2ixuw1200517975.png")
> system("convert tmp/3od8b1200517975.ps tmp/3od8b1200517975.png")
> system("convert tmp/4nmyi1200517975.ps tmp/4nmyi1200517975.png")
> system("convert tmp/5sucb1200517975.ps tmp/5sucb1200517975.png")
> system("convert tmp/6puh41200517975.ps tmp/6puh41200517975.png")
> system("convert tmp/7109l1200517975.ps tmp/7109l1200517975.png")
> system("convert tmp/86h6y1200517975.ps tmp/86h6y1200517975.png")
> system("convert tmp/9jx5k1200517975.ps tmp/9jx5k1200517975.png")
>
>
> proc.time()
user system elapsed
2.254 1.513 5.022