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(523000,519000,509000,512000,519000,517000,510000,509000,501000,507000,569000,580000,578000,565000,547000,555000,562000,561000,555000,544000,537000,543000,594000,611000,613000,611000,594000,595000,591000,589000,584000,573000,567000,569000,621000,629000,628000,612000,595000,597000,593000,590000,580000,574000,573000,573000,620000,626000,620000,588000,566000,557000,561000,549000,532000,526000,511000,499000,555000,565000,542000),dim=c(1,61),dimnames=list(c('Werkloosheid_(aantal)'),1:61))
> y <- array(NA,dim=c(1,61),dimnames=list(c('Werkloosheid_(aantal)'),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 = '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
Werkloosheid_(aantal) M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 523000 1 0 0 0 0 0 0 0 0 0 0
2 519000 0 1 0 0 0 0 0 0 0 0 0
3 509000 0 0 1 0 0 0 0 0 0 0 0
4 512000 0 0 0 1 0 0 0 0 0 0 0
5 519000 0 0 0 0 1 0 0 0 0 0 0
6 517000 0 0 0 0 0 1 0 0 0 0 0
7 510000 0 0 0 0 0 0 1 0 0 0 0
8 509000 0 0 0 0 0 0 0 1 0 0 0
9 501000 0 0 0 0 0 0 0 0 1 0 0
10 507000 0 0 0 0 0 0 0 0 0 1 0
11 569000 0 0 0 0 0 0 0 0 0 0 1
12 580000 0 0 0 0 0 0 0 0 0 0 0
13 578000 1 0 0 0 0 0 0 0 0 0 0
14 565000 0 1 0 0 0 0 0 0 0 0 0
15 547000 0 0 1 0 0 0 0 0 0 0 0
16 555000 0 0 0 1 0 0 0 0 0 0 0
17 562000 0 0 0 0 1 0 0 0 0 0 0
18 561000 0 0 0 0 0 1 0 0 0 0 0
19 555000 0 0 0 0 0 0 1 0 0 0 0
20 544000 0 0 0 0 0 0 0 1 0 0 0
21 537000 0 0 0 0 0 0 0 0 1 0 0
22 543000 0 0 0 0 0 0 0 0 0 1 0
23 594000 0 0 0 0 0 0 0 0 0 0 1
24 611000 0 0 0 0 0 0 0 0 0 0 0
25 613000 1 0 0 0 0 0 0 0 0 0 0
26 611000 0 1 0 0 0 0 0 0 0 0 0
27 594000 0 0 1 0 0 0 0 0 0 0 0
28 595000 0 0 0 1 0 0 0 0 0 0 0
29 591000 0 0 0 0 1 0 0 0 0 0 0
30 589000 0 0 0 0 0 1 0 0 0 0 0
31 584000 0 0 0 0 0 0 1 0 0 0 0
32 573000 0 0 0 0 0 0 0 1 0 0 0
33 567000 0 0 0 0 0 0 0 0 1 0 0
34 569000 0 0 0 0 0 0 0 0 0 1 0
35 621000 0 0 0 0 0 0 0 0 0 0 1
36 629000 0 0 0 0 0 0 0 0 0 0 0
37 628000 1 0 0 0 0 0 0 0 0 0 0
38 612000 0 1 0 0 0 0 0 0 0 0 0
39 595000 0 0 1 0 0 0 0 0 0 0 0
40 597000 0 0 0 1 0 0 0 0 0 0 0
41 593000 0 0 0 0 1 0 0 0 0 0 0
42 590000 0 0 0 0 0 1 0 0 0 0 0
43 580000 0 0 0 0 0 0 1 0 0 0 0
44 574000 0 0 0 0 0 0 0 1 0 0 0
45 573000 0 0 0 0 0 0 0 0 1 0 0
46 573000 0 0 0 0 0 0 0 0 0 1 0
47 620000 0 0 0 0 0 0 0 0 0 0 1
48 626000 0 0 0 0 0 0 0 0 0 0 0
49 620000 1 0 0 0 0 0 0 0 0 0 0
50 588000 0 1 0 0 0 0 0 0 0 0 0
51 566000 0 0 1 0 0 0 0 0 0 0 0
52 557000 0 0 0 1 0 0 0 0 0 0 0
53 561000 0 0 0 0 1 0 0 0 0 0 0
54 549000 0 0 0 0 0 1 0 0 0 0 0
55 532000 0 0 0 0 0 0 1 0 0 0 0
56 526000 0 0 0 0 0 0 0 1 0 0 0
57 511000 0 0 0 0 0 0 0 0 1 0 0
58 499000 0 0 0 0 0 0 0 0 0 1 0
59 555000 0 0 0 0 0 0 0 0 0 0 1
60 565000 0 0 0 0 0 0 0 0 0 0 0
61 542000 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
602200 -18200 -23200 -40000 -39000 -37000
M6 M7 M8 M9 M10 M11
-41000 -50000 -57000 -64400 -64000 -10400
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-61000 -22800 2200 28800 44000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 602200 15094 39.896 < 2e-16 ***
M1 -18200 20438 -0.891 0.37754
M2 -23200 21346 -1.087 0.28243
M3 -40000 21346 -1.874 0.06692 .
M4 -39000 21346 -1.827 0.07379 .
M5 -37000 21346 -1.733 0.08933 .
M6 -41000 21346 -1.921 0.06060 .
M7 -50000 21346 -2.342 0.02327 *
M8 -57000 21346 -2.670 0.01026 *
M9 -64400 21346 -3.017 0.00404 **
M10 -64000 21346 -2.998 0.00426 **
M11 -10400 21346 -0.487 0.62829
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 33750 on 49 degrees of freedom
Multiple R-Squared: 0.2993, Adjusted R-squared: 0.142
F-statistic: 1.902 on 11 and 49 DF, p-value: 0.06207
> postscript(file="/var/www/html/rcomp/tmp/193a11197454668.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/28hxw1197454669.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/344zd1197454669.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/42bpj1197454669.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/5widu1197454669.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 7 8 9 10 11
-61000 -60000 -53200 -51200 -46200 -44200 -42200 -36200 -36800 -31200 -22800
12 13 14 15 16 17 18 19 20 21 22
-22200 -6000 -14000 -15200 -8200 -3200 -200 2800 -1200 -800 4800
23 24 25 26 27 28 29 30 31 32 33
2200 8800 29000 32000 31800 31800 25800 27800 31800 27800 29200
34 35 36 37 38 39 40 41 42 43 44
30800 29200 26800 44000 33000 32800 33800 27800 28800 27800 28800
45 46 47 48 49 50 51 52 53 54 55
35200 34800 28200 23800 36000 9000 3800 -6200 -4200 -12200 -20200
56 57 58 59 60 61
-19200 -26800 -39200 -36800 -37200 -42000
> postscript(file="/var/www/html/rcomp/tmp/6wqd41197454669.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 -61000 NA
1 -60000 -61000
2 -53200 -60000
3 -51200 -53200
4 -46200 -51200
5 -44200 -46200
6 -42200 -44200
7 -36200 -42200
8 -36800 -36200
9 -31200 -36800
10 -22800 -31200
11 -22200 -22800
12 -6000 -22200
13 -14000 -6000
14 -15200 -14000
15 -8200 -15200
16 -3200 -8200
17 -200 -3200
18 2800 -200
19 -1200 2800
20 -800 -1200
21 4800 -800
22 2200 4800
23 8800 2200
24 29000 8800
25 32000 29000
26 31800 32000
27 31800 31800
28 25800 31800
29 27800 25800
30 31800 27800
31 27800 31800
32 29200 27800
33 30800 29200
34 29200 30800
35 26800 29200
36 44000 26800
37 33000 44000
38 32800 33000
39 33800 32800
40 27800 33800
41 28800 27800
42 27800 28800
43 28800 27800
44 35200 28800
45 34800 35200
46 28200 34800
47 23800 28200
48 36000 23800
49 9000 36000
50 3800 9000
51 -6200 3800
52 -4200 -6200
53 -12200 -4200
54 -20200 -12200
55 -19200 -20200
56 -26800 -19200
57 -39200 -26800
58 -36800 -39200
59 -37200 -36800
60 -42000 -37200
61 NA -42000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -60000 -61000
[2,] -53200 -60000
[3,] -51200 -53200
[4,] -46200 -51200
[5,] -44200 -46200
[6,] -42200 -44200
[7,] -36200 -42200
[8,] -36800 -36200
[9,] -31200 -36800
[10,] -22800 -31200
[11,] -22200 -22800
[12,] -6000 -22200
[13,] -14000 -6000
[14,] -15200 -14000
[15,] -8200 -15200
[16,] -3200 -8200
[17,] -200 -3200
[18,] 2800 -200
[19,] -1200 2800
[20,] -800 -1200
[21,] 4800 -800
[22,] 2200 4800
[23,] 8800 2200
[24,] 29000 8800
[25,] 32000 29000
[26,] 31800 32000
[27,] 31800 31800
[28,] 25800 31800
[29,] 27800 25800
[30,] 31800 27800
[31,] 27800 31800
[32,] 29200 27800
[33,] 30800 29200
[34,] 29200 30800
[35,] 26800 29200
[36,] 44000 26800
[37,] 33000 44000
[38,] 32800 33000
[39,] 33800 32800
[40,] 27800 33800
[41,] 28800 27800
[42,] 27800 28800
[43,] 28800 27800
[44,] 35200 28800
[45,] 34800 35200
[46,] 28200 34800
[47,] 23800 28200
[48,] 36000 23800
[49,] 9000 36000
[50,] 3800 9000
[51,] -6200 3800
[52,] -4200 -6200
[53,] -12200 -4200
[54,] -20200 -12200
[55,] -19200 -20200
[56,] -26800 -19200
[57,] -39200 -26800
[58,] -36800 -39200
[59,] -37200 -36800
[60,] -42000 -37200
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -60000 -61000
2 -53200 -60000
3 -51200 -53200
4 -46200 -51200
5 -44200 -46200
6 -42200 -44200
7 -36200 -42200
8 -36800 -36200
9 -31200 -36800
10 -22800 -31200
11 -22200 -22800
12 -6000 -22200
13 -14000 -6000
14 -15200 -14000
15 -8200 -15200
16 -3200 -8200
17 -200 -3200
18 2800 -200
19 -1200 2800
20 -800 -1200
21 4800 -800
22 2200 4800
23 8800 2200
24 29000 8800
25 32000 29000
26 31800 32000
27 31800 31800
28 25800 31800
29 27800 25800
30 31800 27800
31 27800 31800
32 29200 27800
33 30800 29200
34 29200 30800
35 26800 29200
36 44000 26800
37 33000 44000
38 32800 33000
39 33800 32800
40 27800 33800
41 28800 27800
42 27800 28800
43 28800 27800
44 35200 28800
45 34800 35200
46 28200 34800
47 23800 28200
48 36000 23800
49 9000 36000
50 3800 9000
51 -6200 3800
52 -4200 -6200
53 -12200 -4200
54 -20200 -12200
55 -19200 -20200
56 -26800 -19200
57 -39200 -26800
58 -36800 -39200
59 -37200 -36800
60 -42000 -37200
> 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/75l8u1197454669.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/8gt5k1197454669.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/9g6pe1197454669.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/105dm01197454669.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/11d6mp1197454669.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/12f81b1197454669.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/13527x1197454669.tab")
>
> system("convert tmp/193a11197454668.ps tmp/193a11197454668.png")
> system("convert tmp/28hxw1197454669.ps tmp/28hxw1197454669.png")
> system("convert tmp/344zd1197454669.ps tmp/344zd1197454669.png")
> system("convert tmp/42bpj1197454669.ps tmp/42bpj1197454669.png")
> system("convert tmp/5widu1197454669.ps tmp/5widu1197454669.png")
> system("convert tmp/6wqd41197454669.ps tmp/6wqd41197454669.png")
> system("convert tmp/75l8u1197454669.ps tmp/75l8u1197454669.png")
> system("convert tmp/8gt5k1197454669.ps tmp/8gt5k1197454669.png")
> system("convert tmp/9g6pe1197454669.ps tmp/9g6pe1197454669.png")
>
>
> proc.time()
user system elapsed
2.567 1.677 4.946