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.
Natural language support but running in an English locale
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(595,0,597,0,593,0,590,0,580,0,574,0,573,0,573,0,620,0,626,0,620,0,588,0,566,0,557,0,561,0,549,0,532,0,526,0,511,0,499,0,555,0,565,0,542,0,527,0,510,0,514,0,517,0,508,0,493,0,490,0,469,0,478,0,528,0,534,0,518,0,506,0,502,1,516,1,528,1,533,1,536,1,537,1,524,1,536,1,587,1,597,1,581,1,564,1,558,1,575,1,580,1,575,1,563,1,552,1,537,1,545,1,601,1,604,1,586,1,564,1,549,1),dim=c(2,61),dimnames=list(c('werkloosheid','X'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('werkloosheid','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])
+ }
+ }
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = ''
> par9 = ''
> par8 = ''
> par7 = ''
> par6 = ''
> par5 = ''
> par4 = ''
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> ylab = ''
> xlab = ''
> main = ''
> #'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)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> 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 X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 595 0 1 0 0 0 0 0 0 0 0 0 0 1
2 597 0 0 1 0 0 0 0 0 0 0 0 0 2
3 593 0 0 0 1 0 0 0 0 0 0 0 0 3
4 590 0 0 0 0 1 0 0 0 0 0 0 0 4
5 580 0 0 0 0 0 1 0 0 0 0 0 0 5
6 574 0 0 0 0 0 0 1 0 0 0 0 0 6
7 573 0 0 0 0 0 0 0 1 0 0 0 0 7
8 573 0 0 0 0 0 0 0 0 1 0 0 0 8
9 620 0 0 0 0 0 0 0 0 0 1 0 0 9
10 626 0 0 0 0 0 0 0 0 0 0 1 0 10
11 620 0 0 0 0 0 0 0 0 0 0 0 1 11
12 588 0 0 0 0 0 0 0 0 0 0 0 0 12
13 566 0 1 0 0 0 0 0 0 0 0 0 0 13
14 557 0 0 1 0 0 0 0 0 0 0 0 0 14
15 561 0 0 0 1 0 0 0 0 0 0 0 0 15
16 549 0 0 0 0 1 0 0 0 0 0 0 0 16
17 532 0 0 0 0 0 1 0 0 0 0 0 0 17
18 526 0 0 0 0 0 0 1 0 0 0 0 0 18
19 511 0 0 0 0 0 0 0 1 0 0 0 0 19
20 499 0 0 0 0 0 0 0 0 1 0 0 0 20
21 555 0 0 0 0 0 0 0 0 0 1 0 0 21
22 565 0 0 0 0 0 0 0 0 0 0 1 0 22
23 542 0 0 0 0 0 0 0 0 0 0 0 1 23
24 527 0 0 0 0 0 0 0 0 0 0 0 0 24
25 510 0 1 0 0 0 0 0 0 0 0 0 0 25
26 514 0 0 1 0 0 0 0 0 0 0 0 0 26
27 517 0 0 0 1 0 0 0 0 0 0 0 0 27
28 508 0 0 0 0 1 0 0 0 0 0 0 0 28
29 493 0 0 0 0 0 1 0 0 0 0 0 0 29
30 490 0 0 0 0 0 0 1 0 0 0 0 0 30
31 469 0 0 0 0 0 0 0 1 0 0 0 0 31
32 478 0 0 0 0 0 0 0 0 1 0 0 0 32
33 528 0 0 0 0 0 0 0 0 0 1 0 0 33
34 534 0 0 0 0 0 0 0 0 0 0 1 0 34
35 518 0 0 0 0 0 0 0 0 0 0 0 1 35
36 506 0 0 0 0 0 0 0 0 0 0 0 0 36
37 502 1 1 0 0 0 0 0 0 0 0 0 0 37
38 516 1 0 1 0 0 0 0 0 0 0 0 0 38
39 528 1 0 0 1 0 0 0 0 0 0 0 0 39
40 533 1 0 0 0 1 0 0 0 0 0 0 0 40
41 536 1 0 0 0 0 1 0 0 0 0 0 0 41
42 537 1 0 0 0 0 0 1 0 0 0 0 0 42
43 524 1 0 0 0 0 0 0 1 0 0 0 0 43
44 536 1 0 0 0 0 0 0 0 1 0 0 0 44
45 587 1 0 0 0 0 0 0 0 0 1 0 0 45
46 597 1 0 0 0 0 0 0 0 0 0 1 0 46
47 581 1 0 0 0 0 0 0 0 0 0 0 1 47
48 564 1 0 0 0 0 0 0 0 0 0 0 0 48
49 558 1 1 0 0 0 0 0 0 0 0 0 0 49
50 575 1 0 1 0 0 0 0 0 0 0 0 0 50
51 580 1 0 0 1 0 0 0 0 0 0 0 0 51
52 575 1 0 0 0 1 0 0 0 0 0 0 0 52
53 563 1 0 0 0 0 1 0 0 0 0 0 0 53
54 552 1 0 0 0 0 0 1 0 0 0 0 0 54
55 537 1 0 0 0 0 0 0 1 0 0 0 0 55
56 545 1 0 0 0 0 0 0 0 1 0 0 0 56
57 601 1 0 0 0 0 0 0 0 0 1 0 0 57
58 604 1 0 0 0 0 0 0 0 0 0 1 0 58
59 586 1 0 0 0 0 0 0 0 0 0 0 1 59
60 564 1 0 0 0 0 0 0 0 0 0 0 0 60
61 549 1 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
602.686 84.276 -23.588 -22.054 -15.649 -18.044
M5 M6 M7 M8 M9 M10
-25.838 -28.433 -39.027 -33.222 21.184 30.589
M11 t
17.195 -2.405
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-72.37 -15.36 -4.09 21.18 32.36
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 602.6859 15.2911 39.414 < 2e-16 ***
X 84.2756 13.8107 6.102 1.88e-07 ***
M1 -23.5881 16.3434 -1.443 0.1556
M2 -22.0545 17.1847 -1.283 0.2057
M3 -15.6490 17.1001 -0.915 0.3648
M4 -18.0436 17.0240 -1.060 0.2946
M5 -25.8381 16.9566 -1.524 0.1343
M6 -28.4327 16.8979 -1.683 0.0991 .
M7 -39.0272 16.8482 -2.316 0.0249 *
M8 -33.2218 16.8073 -1.977 0.0540 .
M9 21.1837 16.7755 1.263 0.2129
M10 30.5891 16.7527 1.826 0.0742 .
M11 17.1946 16.7390 1.027 0.3096
t -2.4054 0.3908 -6.155 1.57e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 26.46 on 47 degrees of freedom
Multiple R-squared: 0.5953, Adjusted R-squared: 0.4833
F-statistic: 5.317 on 13 and 47 DF, p-value: 9.672e-06
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.04922898 0.09845795 0.95077102
[2,] 0.02776624 0.05553248 0.97223376
[3,] 0.07016244 0.14032488 0.92983756
[4,] 0.17593625 0.35187251 0.82406375
[5,] 0.16931489 0.33862979 0.83068511
[6,] 0.15539227 0.31078454 0.84460773
[7,] 0.24983475 0.49966951 0.75016525
[8,] 0.26606800 0.53213600 0.73393200
[9,] 0.34704703 0.69409406 0.65295297
[10,] 0.38560874 0.77121748 0.61439126
[11,] 0.41906819 0.83813638 0.58093181
[12,] 0.39037581 0.78075161 0.60962419
[13,] 0.31859411 0.63718823 0.68140589
[14,] 0.26767868 0.53535736 0.73232132
[15,] 0.19365217 0.38730434 0.80634783
[16,] 0.15233439 0.30466877 0.84766561
[17,] 0.11082780 0.22165560 0.88917220
[18,] 0.07472666 0.14945331 0.92527334
[19,] 0.04747625 0.09495249 0.95252375
[20,] 0.03738029 0.07476058 0.96261971
[21,] 0.03270852 0.06541703 0.96729148
[22,] 0.11201173 0.22402346 0.88798827
[23,] 0.38276430 0.76552860 0.61723570
[24,] 0.80127104 0.39745793 0.19872896
[25,] 0.93939961 0.12120078 0.06060039
[26,] 0.93954362 0.12091276 0.06045638
[27,] 0.91637153 0.16725693 0.08362847
[28,] 0.85994476 0.28011049 0.14005524
> postscript(file="/var/www/html/freestat/rcomp/tmp/1h6dx1293210298.ps",horizontal=F,onefile=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/freestat/rcomp/tmp/2h6dx1293210298.ps",horizontal=F,onefile=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/freestat/rcomp/tmp/3h6dx1293210298.ps",horizontal=F,onefile=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/freestat/rcomp/tmp/4axcz1293210298.ps",horizontal=F,onefile=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/freestat/rcomp/tmp/5axcz1293210298.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> 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
18.307692 21.179487 13.179487 14.979487 15.179487 14.179487 26.179487
8 9 10 11 12 13 14
22.779487 17.779487 16.779487 26.579487 14.179487 18.173077 10.044872
15 16 17 18 19 20 21
10.044872 2.844872 -3.955128 -4.955128 -6.955128 -22.355128 -18.355128
22 23 24 25 26 27 28
-15.355128 -22.555128 -17.955128 -8.961538 -4.089744 -5.089744 -9.289744
29 30 31 32 33 34 35
-14.089744 -12.089744 -20.089744 -14.489744 -16.489744 -17.489744 -17.689744
36 37 38 39 40 41 42
-10.089744 -72.371795 -57.500000 -49.500000 -39.700000 -26.500000 -20.500000
43 44 45 46 47 48 49
-20.500000 -11.900000 -12.900000 -9.900000 -10.100000 -7.500000 12.493590
50 51 52 53 54 55 56
30.365385 31.365385 31.165385 29.365385 23.365385 21.365385 25.965385
57 58 59 60 61
29.965385 25.965385 23.765385 21.365385 32.358974
> postscript(file="/var/www/html/freestat/rcomp/tmp/6axcz1293210298.ps",horizontal=F,onefile=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 18.307692 NA
1 21.179487 18.307692
2 13.179487 21.179487
3 14.979487 13.179487
4 15.179487 14.979487
5 14.179487 15.179487
6 26.179487 14.179487
7 22.779487 26.179487
8 17.779487 22.779487
9 16.779487 17.779487
10 26.579487 16.779487
11 14.179487 26.579487
12 18.173077 14.179487
13 10.044872 18.173077
14 10.044872 10.044872
15 2.844872 10.044872
16 -3.955128 2.844872
17 -4.955128 -3.955128
18 -6.955128 -4.955128
19 -22.355128 -6.955128
20 -18.355128 -22.355128
21 -15.355128 -18.355128
22 -22.555128 -15.355128
23 -17.955128 -22.555128
24 -8.961538 -17.955128
25 -4.089744 -8.961538
26 -5.089744 -4.089744
27 -9.289744 -5.089744
28 -14.089744 -9.289744
29 -12.089744 -14.089744
30 -20.089744 -12.089744
31 -14.489744 -20.089744
32 -16.489744 -14.489744
33 -17.489744 -16.489744
34 -17.689744 -17.489744
35 -10.089744 -17.689744
36 -72.371795 -10.089744
37 -57.500000 -72.371795
38 -49.500000 -57.500000
39 -39.700000 -49.500000
40 -26.500000 -39.700000
41 -20.500000 -26.500000
42 -20.500000 -20.500000
43 -11.900000 -20.500000
44 -12.900000 -11.900000
45 -9.900000 -12.900000
46 -10.100000 -9.900000
47 -7.500000 -10.100000
48 12.493590 -7.500000
49 30.365385 12.493590
50 31.365385 30.365385
51 31.165385 31.365385
52 29.365385 31.165385
53 23.365385 29.365385
54 21.365385 23.365385
55 25.965385 21.365385
56 29.965385 25.965385
57 25.965385 29.965385
58 23.765385 25.965385
59 21.365385 23.765385
60 32.358974 21.365385
61 NA 32.358974
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 21.179487 18.307692
[2,] 13.179487 21.179487
[3,] 14.979487 13.179487
[4,] 15.179487 14.979487
[5,] 14.179487 15.179487
[6,] 26.179487 14.179487
[7,] 22.779487 26.179487
[8,] 17.779487 22.779487
[9,] 16.779487 17.779487
[10,] 26.579487 16.779487
[11,] 14.179487 26.579487
[12,] 18.173077 14.179487
[13,] 10.044872 18.173077
[14,] 10.044872 10.044872
[15,] 2.844872 10.044872
[16,] -3.955128 2.844872
[17,] -4.955128 -3.955128
[18,] -6.955128 -4.955128
[19,] -22.355128 -6.955128
[20,] -18.355128 -22.355128
[21,] -15.355128 -18.355128
[22,] -22.555128 -15.355128
[23,] -17.955128 -22.555128
[24,] -8.961538 -17.955128
[25,] -4.089744 -8.961538
[26,] -5.089744 -4.089744
[27,] -9.289744 -5.089744
[28,] -14.089744 -9.289744
[29,] -12.089744 -14.089744
[30,] -20.089744 -12.089744
[31,] -14.489744 -20.089744
[32,] -16.489744 -14.489744
[33,] -17.489744 -16.489744
[34,] -17.689744 -17.489744
[35,] -10.089744 -17.689744
[36,] -72.371795 -10.089744
[37,] -57.500000 -72.371795
[38,] -49.500000 -57.500000
[39,] -39.700000 -49.500000
[40,] -26.500000 -39.700000
[41,] -20.500000 -26.500000
[42,] -20.500000 -20.500000
[43,] -11.900000 -20.500000
[44,] -12.900000 -11.900000
[45,] -9.900000 -12.900000
[46,] -10.100000 -9.900000
[47,] -7.500000 -10.100000
[48,] 12.493590 -7.500000
[49,] 30.365385 12.493590
[50,] 31.365385 30.365385
[51,] 31.165385 31.365385
[52,] 29.365385 31.165385
[53,] 23.365385 29.365385
[54,] 21.365385 23.365385
[55,] 25.965385 21.365385
[56,] 29.965385 25.965385
[57,] 25.965385 29.965385
[58,] 23.765385 25.965385
[59,] 21.365385 23.765385
[60,] 32.358974 21.365385
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 21.179487 18.307692
2 13.179487 21.179487
3 14.979487 13.179487
4 15.179487 14.979487
5 14.179487 15.179487
6 26.179487 14.179487
7 22.779487 26.179487
8 17.779487 22.779487
9 16.779487 17.779487
10 26.579487 16.779487
11 14.179487 26.579487
12 18.173077 14.179487
13 10.044872 18.173077
14 10.044872 10.044872
15 2.844872 10.044872
16 -3.955128 2.844872
17 -4.955128 -3.955128
18 -6.955128 -4.955128
19 -22.355128 -6.955128
20 -18.355128 -22.355128
21 -15.355128 -18.355128
22 -22.555128 -15.355128
23 -17.955128 -22.555128
24 -8.961538 -17.955128
25 -4.089744 -8.961538
26 -5.089744 -4.089744
27 -9.289744 -5.089744
28 -14.089744 -9.289744
29 -12.089744 -14.089744
30 -20.089744 -12.089744
31 -14.489744 -20.089744
32 -16.489744 -14.489744
33 -17.489744 -16.489744
34 -17.689744 -17.489744
35 -10.089744 -17.689744
36 -72.371795 -10.089744
37 -57.500000 -72.371795
38 -49.500000 -57.500000
39 -39.700000 -49.500000
40 -26.500000 -39.700000
41 -20.500000 -26.500000
42 -20.500000 -20.500000
43 -11.900000 -20.500000
44 -12.900000 -11.900000
45 -9.900000 -12.900000
46 -10.100000 -9.900000
47 -7.500000 -10.100000
48 12.493590 -7.500000
49 30.365385 12.493590
50 31.365385 30.365385
51 31.165385 31.365385
52 29.365385 31.165385
53 23.365385 29.365385
54 21.365385 23.365385
55 25.965385 21.365385
56 29.965385 25.965385
57 25.965385 29.965385
58 23.765385 25.965385
59 21.365385 23.765385
60 32.358974 21.365385
> 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/freestat/rcomp/tmp/73pbk1293210298.ps",horizontal=F,onefile=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/freestat/rcomp/tmp/8oqd01293210299.ps",horizontal=F,onefile=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/freestat/rcomp/tmp/9oqd01293210299.ps",horizontal=F,onefile=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
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10oqd01293210299.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/1120tr1293210299.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/freestat/rcomp/tmp/12v9au1293210299.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/freestat/rcomp/tmp/131s7o1293210299.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/freestat/rcomp/tmp/14u1or1293210299.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/15g25x1293210299.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/16cc2n1293210299.tab")
+ }
>
> try(system("convert tmp/1h6dx1293210298.ps tmp/1h6dx1293210298.png",intern=TRUE))
character(0)
> try(system("convert tmp/2h6dx1293210298.ps tmp/2h6dx1293210298.png",intern=TRUE))
character(0)
> try(system("convert tmp/3h6dx1293210298.ps tmp/3h6dx1293210298.png",intern=TRUE))
character(0)
> try(system("convert tmp/4axcz1293210298.ps tmp/4axcz1293210298.png",intern=TRUE))
character(0)
> try(system("convert tmp/5axcz1293210298.ps tmp/5axcz1293210298.png",intern=TRUE))
character(0)
> try(system("convert tmp/6axcz1293210298.ps tmp/6axcz1293210298.png",intern=TRUE))
character(0)
> try(system("convert tmp/73pbk1293210298.ps tmp/73pbk1293210298.png",intern=TRUE))
character(0)
> try(system("convert tmp/8oqd01293210299.ps tmp/8oqd01293210299.png",intern=TRUE))
character(0)
> try(system("convert tmp/9oqd01293210299.ps tmp/9oqd01293210299.png",intern=TRUE))
character(0)
> try(system("convert tmp/10oqd01293210299.ps tmp/10oqd01293210299.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.762 2.426 4.089