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.
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(8.4
+ ,7.6
+ ,9.5
+ ,25
+ ,6.6
+ ,8.4
+ ,7.9
+ ,9.1
+ ,23.6
+ ,6.7
+ ,8.4
+ ,7.9
+ ,9
+ ,22.3
+ ,6.8
+ ,8.6
+ ,8.1
+ ,9.3
+ ,21.8
+ ,7.2
+ ,8.9
+ ,8.2
+ ,9.9
+ ,20.8
+ ,7.6
+ ,8.8
+ ,8
+ ,9.8
+ ,19.7
+ ,7.6
+ ,8.3
+ ,7.5
+ ,9.4
+ ,18.3
+ ,7.3
+ ,7.5
+ ,6.8
+ ,8.3
+ ,17.4
+ ,6.4
+ ,7.2
+ ,6.5
+ ,8
+ ,17
+ ,6.1
+ ,7.5
+ ,6.6
+ ,8.5
+ ,18.1
+ ,6.3
+ ,8.8
+ ,7.6
+ ,10.4
+ ,23.9
+ ,7.1
+ ,9.3
+ ,8
+ ,11.1
+ ,25.6
+ ,7.5
+ ,9.3
+ ,8
+ ,10.9
+ ,25.3
+ ,7.4
+ ,8.7
+ ,7.7
+ ,9.9
+ ,23.6
+ ,7.1
+ ,8.2
+ ,7.5
+ ,9.2
+ ,21.9
+ ,6.8
+ ,8.3
+ ,7.6
+ ,9.2
+ ,21.4
+ ,6.9
+ ,8.5
+ ,7.7
+ ,9.5
+ ,20.6
+ ,7.2
+ ,8.6
+ ,7.9
+ ,9.6
+ ,20.5
+ ,7.4
+ ,8.6
+ ,7.8
+ ,9.5
+ ,20.2
+ ,7.3
+ ,8.2
+ ,7.5
+ ,9.1
+ ,20.6
+ ,6.9
+ ,8.1
+ ,7.5
+ ,8.9
+ ,19.7
+ ,6.9
+ ,8
+ ,7.1
+ ,9
+ ,19.3
+ ,6.8
+ ,8.6
+ ,7.5
+ ,10.1
+ ,22.8
+ ,7.1
+ ,8.7
+ ,7.5
+ ,10.3
+ ,23.5
+ ,7.2
+ ,8.8
+ ,7.6
+ ,10.2
+ ,23.8
+ ,7.1
+ ,8.5
+ ,7.7
+ ,9.6
+ ,22.6
+ ,7
+ ,8.4
+ ,7.7
+ ,9.2
+ ,22
+ ,6.9
+ ,8.5
+ ,7.9
+ ,9.3
+ ,21.7
+ ,7
+ ,8.7
+ ,8.1
+ ,9.4
+ ,20.7
+ ,7.4
+ ,8.7
+ ,8.2
+ ,9.4
+ ,20.2
+ ,7.5
+ ,8.6
+ ,8.2
+ ,9.2
+ ,19.1
+ ,7.5
+ ,8.5
+ ,8.1
+ ,9
+ ,19.5
+ ,7.4
+ ,8.3
+ ,7.9
+ ,9
+ ,18.7
+ ,7.3
+ ,8.1
+ ,7.3
+ ,9
+ ,18.6
+ ,7
+ ,8.2
+ ,6.9
+ ,9.8
+ ,22.2
+ ,6.7
+ ,8.1
+ ,6.6
+ ,10
+ ,23.2
+ ,6.5
+ ,8.1
+ ,6.7
+ ,9.9
+ ,23.5
+ ,6.5
+ ,7.9
+ ,6.9
+ ,9.3
+ ,21.3
+ ,6.5
+ ,7.9
+ ,7
+ ,9
+ ,20
+ ,6.6
+ ,7.9
+ ,7.1
+ ,9
+ ,18.7
+ ,6.8
+ ,8
+ ,7.2
+ ,9.1
+ ,18.9
+ ,6.9
+ ,8
+ ,7.1
+ ,9.1
+ ,18.3
+ ,6.9
+ ,7.9
+ ,6.9
+ ,9.1
+ ,18.4
+ ,6.8
+ ,8
+ ,7
+ ,9.2
+ ,19.9
+ ,6.8
+ ,7.7
+ ,6.8
+ ,8.8
+ ,19.2
+ ,6.5
+ ,7.2
+ ,6.4
+ ,8.3
+ ,18.5
+ ,6.1
+ ,7.5
+ ,6.7
+ ,8.4
+ ,20.9
+ ,6
+ ,7.3
+ ,6.7
+ ,8.1
+ ,20.5
+ ,5.9
+ ,7
+ ,6.4
+ ,7.8
+ ,19.4
+ ,5.8
+ ,7
+ ,6.3
+ ,7.9
+ ,18.1
+ ,5.9
+ ,7
+ ,6.2
+ ,7.9
+ ,17
+ ,5.9
+ ,7.2
+ ,6.5
+ ,8
+ ,17
+ ,6.2
+ ,7.3
+ ,6.8
+ ,7.9
+ ,17.3
+ ,6.3
+ ,7.1
+ ,6.8
+ ,7.5
+ ,16.7
+ ,6.2
+ ,6.8
+ ,6.5
+ ,7.2
+ ,15.5
+ ,6
+ ,6.6
+ ,6.3
+ ,6.9
+ ,15.3
+ ,5.8
+ ,6.2
+ ,5.9
+ ,6.6
+ ,13.7
+ ,5.5
+ ,6.2
+ ,5.9
+ ,6.7
+ ,14.1
+ ,5.5
+ ,6.8
+ ,6.4
+ ,7.3
+ ,17.3
+ ,5.7
+ ,6.9
+ ,6.4
+ ,7.5
+ ,18.1
+ ,5.8
+ ,6.8
+ ,6.5
+ ,7.2
+ ,18.1
+ ,5.7)
+ ,dim=c(5
+ ,61)
+ ,dimnames=list(c('Totaal'
+ ,'Mannen'
+ ,'Vrouwen'
+ ,'<25j'
+ ,'>25j')
+ ,1:61))
> y <- array(NA,dim=c(5,61),dimnames=list(c('Totaal','Mannen','Vrouwen','<25j','>25j'),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 = '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)
> 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
Totaal Mannen Vrouwen <25j >25j M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 8.4 7.6 9.5 25.0 6.6 1 0 0 0 0 0 0 0 0 0 0 1
2 8.4 7.9 9.1 23.6 6.7 0 1 0 0 0 0 0 0 0 0 0 2
3 8.4 7.9 9.0 22.3 6.8 0 0 1 0 0 0 0 0 0 0 0 3
4 8.6 8.1 9.3 21.8 7.2 0 0 0 1 0 0 0 0 0 0 0 4
5 8.9 8.2 9.9 20.8 7.6 0 0 0 0 1 0 0 0 0 0 0 5
6 8.8 8.0 9.8 19.7 7.6 0 0 0 0 0 1 0 0 0 0 0 6
7 8.3 7.5 9.4 18.3 7.3 0 0 0 0 0 0 1 0 0 0 0 7
8 7.5 6.8 8.3 17.4 6.4 0 0 0 0 0 0 0 1 0 0 0 8
9 7.2 6.5 8.0 17.0 6.1 0 0 0 0 0 0 0 0 1 0 0 9
10 7.5 6.6 8.5 18.1 6.3 0 0 0 0 0 0 0 0 0 1 0 10
11 8.8 7.6 10.4 23.9 7.1 0 0 0 0 0 0 0 0 0 0 1 11
12 9.3 8.0 11.1 25.6 7.5 0 0 0 0 0 0 0 0 0 0 0 12
13 9.3 8.0 10.9 25.3 7.4 1 0 0 0 0 0 0 0 0 0 0 13
14 8.7 7.7 9.9 23.6 7.1 0 1 0 0 0 0 0 0 0 0 0 14
15 8.2 7.5 9.2 21.9 6.8 0 0 1 0 0 0 0 0 0 0 0 15
16 8.3 7.6 9.2 21.4 6.9 0 0 0 1 0 0 0 0 0 0 0 16
17 8.5 7.7 9.5 20.6 7.2 0 0 0 0 1 0 0 0 0 0 0 17
18 8.6 7.9 9.6 20.5 7.4 0 0 0 0 0 1 0 0 0 0 0 18
19 8.6 7.8 9.5 20.2 7.3 0 0 0 0 0 0 1 0 0 0 0 19
20 8.2 7.5 9.1 20.6 6.9 0 0 0 0 0 0 0 1 0 0 0 20
21 8.1 7.5 8.9 19.7 6.9 0 0 0 0 0 0 0 0 1 0 0 21
22 8.0 7.1 9.0 19.3 6.8 0 0 0 0 0 0 0 0 0 1 0 22
23 8.6 7.5 10.1 22.8 7.1 0 0 0 0 0 0 0 0 0 0 1 23
24 8.7 7.5 10.3 23.5 7.2 0 0 0 0 0 0 0 0 0 0 0 24
25 8.8 7.6 10.2 23.8 7.1 1 0 0 0 0 0 0 0 0 0 0 25
26 8.5 7.7 9.6 22.6 7.0 0 1 0 0 0 0 0 0 0 0 0 26
27 8.4 7.7 9.2 22.0 6.9 0 0 1 0 0 0 0 0 0 0 0 27
28 8.5 7.9 9.3 21.7 7.0 0 0 0 1 0 0 0 0 0 0 0 28
29 8.7 8.1 9.4 20.7 7.4 0 0 0 0 1 0 0 0 0 0 0 29
30 8.7 8.2 9.4 20.2 7.5 0 0 0 0 0 1 0 0 0 0 0 30
31 8.6 8.2 9.2 19.1 7.5 0 0 0 0 0 0 1 0 0 0 0 31
32 8.5 8.1 9.0 19.5 7.4 0 0 0 0 0 0 0 1 0 0 0 32
33 8.3 7.9 9.0 18.7 7.3 0 0 0 0 0 0 0 0 1 0 0 33
34 8.1 7.3 9.0 18.6 7.0 0 0 0 0 0 0 0 0 0 1 0 34
35 8.2 6.9 9.8 22.2 6.7 0 0 0 0 0 0 0 0 0 0 1 35
36 8.1 6.6 10.0 23.2 6.5 0 0 0 0 0 0 0 0 0 0 0 36
37 8.1 6.7 9.9 23.5 6.5 1 0 0 0 0 0 0 0 0 0 0 37
38 7.9 6.9 9.3 21.3 6.5 0 1 0 0 0 0 0 0 0 0 0 38
39 7.9 7.0 9.0 20.0 6.6 0 0 1 0 0 0 0 0 0 0 0 39
40 7.9 7.1 9.0 18.7 6.8 0 0 0 1 0 0 0 0 0 0 0 40
41 8.0 7.2 9.1 18.9 6.9 0 0 0 0 1 0 0 0 0 0 0 41
42 8.0 7.1 9.1 18.3 6.9 0 0 0 0 0 1 0 0 0 0 0 42
43 7.9 6.9 9.1 18.4 6.8 0 0 0 0 0 0 1 0 0 0 0 43
44 8.0 7.0 9.2 19.9 6.8 0 0 0 0 0 0 0 1 0 0 0 44
45 7.7 6.8 8.8 19.2 6.5 0 0 0 0 0 0 0 0 1 0 0 45
46 7.2 6.4 8.3 18.5 6.1 0 0 0 0 0 0 0 0 0 1 0 46
47 7.5 6.7 8.4 20.9 6.0 0 0 0 0 0 0 0 0 0 0 1 47
48 7.3 6.7 8.1 20.5 5.9 0 0 0 0 0 0 0 0 0 0 0 48
49 7.0 6.4 7.8 19.4 5.8 1 0 0 0 0 0 0 0 0 0 0 49
50 7.0 6.3 7.9 18.1 5.9 0 1 0 0 0 0 0 0 0 0 0 50
51 7.0 6.2 7.9 17.0 5.9 0 0 1 0 0 0 0 0 0 0 0 51
52 7.2 6.5 8.0 17.0 6.2 0 0 0 1 0 0 0 0 0 0 0 52
53 7.3 6.8 7.9 17.3 6.3 0 0 0 0 1 0 0 0 0 0 0 53
54 7.1 6.8 7.5 16.7 6.2 0 0 0 0 0 1 0 0 0 0 0 54
55 6.8 6.5 7.2 15.5 6.0 0 0 0 0 0 0 1 0 0 0 0 55
56 6.6 6.3 6.9 15.3 5.8 0 0 0 0 0 0 0 1 0 0 0 56
57 6.2 5.9 6.6 13.7 5.5 0 0 0 0 0 0 0 0 1 0 0 57
58 6.2 5.9 6.7 14.1 5.5 0 0 0 0 0 0 0 0 0 1 0 58
59 6.8 6.4 7.3 17.3 5.7 0 0 0 0 0 0 0 0 0 0 1 59
60 6.9 6.4 7.5 18.1 5.8 0 0 0 0 0 0 0 0 0 0 0 60
61 6.8 6.5 7.2 18.1 5.7 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) Mannen Vrouwen `<25j` `>25j` M1
0.0915443 0.4657883 0.3928321 0.0140979 0.1060079 0.0222052
M2 M3 M4 M5 M6 M7
0.0078470 0.0458619 0.0271220 0.0333285 0.0290676 0.0363209
M8 M9 M10 M11 t
0.0483014 0.0190383 0.0366811 0.0245810 -0.0003761
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.068830 -0.027028 0.000954 0.022327 0.060588
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0915443 0.1115857 0.820 0.4164
Mannen 0.4657883 0.0814618 5.718 8.75e-07 ***
Vrouwen 0.3928321 0.0656696 5.982 3.59e-07 ***
`<25j` 0.0140979 0.0163584 0.862 0.3935
`>25j` 0.1060079 0.1332133 0.796 0.4304
M1 0.0222052 0.0229028 0.970 0.3376
M2 0.0078470 0.0245430 0.320 0.7507
M3 0.0458619 0.0264334 1.735 0.0897 .
M4 0.0271220 0.0292726 0.927 0.3592
M5 0.0333285 0.0333358 1.000 0.3229
M6 0.0290676 0.0360310 0.807 0.4242
M7 0.0363209 0.0373102 0.973 0.3356
M8 0.0483014 0.0332356 1.453 0.1532
M9 0.0190383 0.0338650 0.562 0.5768
M10 0.0366811 0.0325349 1.127 0.2657
M11 0.0245810 0.0242353 1.014 0.3160
t -0.0003761 0.0005187 -0.725 0.4723
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.03713 on 44 degrees of freedom
Multiple R-squared: 0.9982, Adjusted R-squared: 0.9975
F-statistic: 1511 on 16 and 44 DF, p-value: < 2.2e-16
> 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.9246275 0.1507450 0.07537249
[2,] 0.8502794 0.2994413 0.14972064
[3,] 0.8311556 0.3376888 0.16884441
[4,] 0.7814725 0.4370550 0.21852751
[5,] 0.7523276 0.4953448 0.24767238
[6,] 0.8060799 0.3878402 0.19392008
[7,] 0.8673765 0.2652469 0.13262347
[8,] 0.8243589 0.3512823 0.17564113
[9,] 0.7815513 0.4368973 0.21844867
[10,] 0.8039541 0.3920918 0.19604592
[11,] 0.7248301 0.5503398 0.27516990
[12,] 0.6630070 0.6739861 0.33699303
[13,] 0.5662933 0.8674134 0.43370672
[14,] 0.7421493 0.5157015 0.25785074
[15,] 0.9612778 0.0774445 0.03872225
[16,] 0.9287001 0.1425999 0.07129994
[17,] 0.8996818 0.2006365 0.10031825
[18,] 0.9486352 0.1027296 0.05136479
[19,] 0.9436805 0.1126390 0.05631948
[20,] 0.9357467 0.1285067 0.06425333
[21,] 0.8665999 0.2668002 0.13340010
[22,] 0.7575895 0.4848210 0.24241049
> postscript(file="/var/www/html/freestat/rcomp/tmp/1jiq51227282650.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/freestat/rcomp/tmp/2mwlg1227282650.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/freestat/rcomp/tmp/3vnu91227282650.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/freestat/rcomp/tmp/43wir1227282650.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/freestat/rcomp/tmp/5l5jg1227282650.ps",horizontal=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
-0.0373679383 0.0038988666 0.0132696052 -0.0139759497 -0.0303897272
6 7 8 9 10
0.0221957856 -0.0431151326 0.0115427210 0.0362095219 0.0392386730
11 12 13 14 15
-0.0270283230 -0.0297385914 0.0418288558 0.0449004023 -0.0688296588
16 17 18 19 20
0.0001556018 0.0093726555 -0.0382230180 0.0555918855 -0.0223792019
21 22 23 24 25
-0.0014855356 0.0445197680 -0.0425795486 -0.0166582257 0.0605884627
26 27 28 29 30
-0.0080386531 0.0305147715 0.0108184210 0.0442418378 -0.0012518616
31 32 33 34 35
-0.0140550790 0.0044473913 -0.0508766865 0.0445417043 0.0101176085
36 37 38 39 40
0.0033484345 -0.0300056662 -0.0417145343 -0.0003562172 -0.0306934507
41 42 43 44 45
-0.0358062769 0.0238682640 0.0193396467 0.0007263919 0.0223269267
46 47 48 49 50
-0.0599368097 0.0502855712 0.0093321547 -0.0288024015 0.0009539185
51 52 53 54 55
0.0254014994 0.0336953776 0.0125815109 -0.0065891700 -0.0177613206
56 57 58 59 60
0.0056626977 -0.0061742266 -0.0683633356 0.0092046919 0.0337162279
61
-0.0062413124
> postscript(file="/var/www/html/freestat/rcomp/tmp/6av621227282650.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 -0.0373679383 NA
1 0.0038988666 -0.0373679383
2 0.0132696052 0.0038988666
3 -0.0139759497 0.0132696052
4 -0.0303897272 -0.0139759497
5 0.0221957856 -0.0303897272
6 -0.0431151326 0.0221957856
7 0.0115427210 -0.0431151326
8 0.0362095219 0.0115427210
9 0.0392386730 0.0362095219
10 -0.0270283230 0.0392386730
11 -0.0297385914 -0.0270283230
12 0.0418288558 -0.0297385914
13 0.0449004023 0.0418288558
14 -0.0688296588 0.0449004023
15 0.0001556018 -0.0688296588
16 0.0093726555 0.0001556018
17 -0.0382230180 0.0093726555
18 0.0555918855 -0.0382230180
19 -0.0223792019 0.0555918855
20 -0.0014855356 -0.0223792019
21 0.0445197680 -0.0014855356
22 -0.0425795486 0.0445197680
23 -0.0166582257 -0.0425795486
24 0.0605884627 -0.0166582257
25 -0.0080386531 0.0605884627
26 0.0305147715 -0.0080386531
27 0.0108184210 0.0305147715
28 0.0442418378 0.0108184210
29 -0.0012518616 0.0442418378
30 -0.0140550790 -0.0012518616
31 0.0044473913 -0.0140550790
32 -0.0508766865 0.0044473913
33 0.0445417043 -0.0508766865
34 0.0101176085 0.0445417043
35 0.0033484345 0.0101176085
36 -0.0300056662 0.0033484345
37 -0.0417145343 -0.0300056662
38 -0.0003562172 -0.0417145343
39 -0.0306934507 -0.0003562172
40 -0.0358062769 -0.0306934507
41 0.0238682640 -0.0358062769
42 0.0193396467 0.0238682640
43 0.0007263919 0.0193396467
44 0.0223269267 0.0007263919
45 -0.0599368097 0.0223269267
46 0.0502855712 -0.0599368097
47 0.0093321547 0.0502855712
48 -0.0288024015 0.0093321547
49 0.0009539185 -0.0288024015
50 0.0254014994 0.0009539185
51 0.0336953776 0.0254014994
52 0.0125815109 0.0336953776
53 -0.0065891700 0.0125815109
54 -0.0177613206 -0.0065891700
55 0.0056626977 -0.0177613206
56 -0.0061742266 0.0056626977
57 -0.0683633356 -0.0061742266
58 0.0092046919 -0.0683633356
59 0.0337162279 0.0092046919
60 -0.0062413124 0.0337162279
61 NA -0.0062413124
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.0038988666 -0.0373679383
[2,] 0.0132696052 0.0038988666
[3,] -0.0139759497 0.0132696052
[4,] -0.0303897272 -0.0139759497
[5,] 0.0221957856 -0.0303897272
[6,] -0.0431151326 0.0221957856
[7,] 0.0115427210 -0.0431151326
[8,] 0.0362095219 0.0115427210
[9,] 0.0392386730 0.0362095219
[10,] -0.0270283230 0.0392386730
[11,] -0.0297385914 -0.0270283230
[12,] 0.0418288558 -0.0297385914
[13,] 0.0449004023 0.0418288558
[14,] -0.0688296588 0.0449004023
[15,] 0.0001556018 -0.0688296588
[16,] 0.0093726555 0.0001556018
[17,] -0.0382230180 0.0093726555
[18,] 0.0555918855 -0.0382230180
[19,] -0.0223792019 0.0555918855
[20,] -0.0014855356 -0.0223792019
[21,] 0.0445197680 -0.0014855356
[22,] -0.0425795486 0.0445197680
[23,] -0.0166582257 -0.0425795486
[24,] 0.0605884627 -0.0166582257
[25,] -0.0080386531 0.0605884627
[26,] 0.0305147715 -0.0080386531
[27,] 0.0108184210 0.0305147715
[28,] 0.0442418378 0.0108184210
[29,] -0.0012518616 0.0442418378
[30,] -0.0140550790 -0.0012518616
[31,] 0.0044473913 -0.0140550790
[32,] -0.0508766865 0.0044473913
[33,] 0.0445417043 -0.0508766865
[34,] 0.0101176085 0.0445417043
[35,] 0.0033484345 0.0101176085
[36,] -0.0300056662 0.0033484345
[37,] -0.0417145343 -0.0300056662
[38,] -0.0003562172 -0.0417145343
[39,] -0.0306934507 -0.0003562172
[40,] -0.0358062769 -0.0306934507
[41,] 0.0238682640 -0.0358062769
[42,] 0.0193396467 0.0238682640
[43,] 0.0007263919 0.0193396467
[44,] 0.0223269267 0.0007263919
[45,] -0.0599368097 0.0223269267
[46,] 0.0502855712 -0.0599368097
[47,] 0.0093321547 0.0502855712
[48,] -0.0288024015 0.0093321547
[49,] 0.0009539185 -0.0288024015
[50,] 0.0254014994 0.0009539185
[51,] 0.0336953776 0.0254014994
[52,] 0.0125815109 0.0336953776
[53,] -0.0065891700 0.0125815109
[54,] -0.0177613206 -0.0065891700
[55,] 0.0056626977 -0.0177613206
[56,] -0.0061742266 0.0056626977
[57,] -0.0683633356 -0.0061742266
[58,] 0.0092046919 -0.0683633356
[59,] 0.0337162279 0.0092046919
[60,] -0.0062413124 0.0337162279
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.0038988666 -0.0373679383
2 0.0132696052 0.0038988666
3 -0.0139759497 0.0132696052
4 -0.0303897272 -0.0139759497
5 0.0221957856 -0.0303897272
6 -0.0431151326 0.0221957856
7 0.0115427210 -0.0431151326
8 0.0362095219 0.0115427210
9 0.0392386730 0.0362095219
10 -0.0270283230 0.0392386730
11 -0.0297385914 -0.0270283230
12 0.0418288558 -0.0297385914
13 0.0449004023 0.0418288558
14 -0.0688296588 0.0449004023
15 0.0001556018 -0.0688296588
16 0.0093726555 0.0001556018
17 -0.0382230180 0.0093726555
18 0.0555918855 -0.0382230180
19 -0.0223792019 0.0555918855
20 -0.0014855356 -0.0223792019
21 0.0445197680 -0.0014855356
22 -0.0425795486 0.0445197680
23 -0.0166582257 -0.0425795486
24 0.0605884627 -0.0166582257
25 -0.0080386531 0.0605884627
26 0.0305147715 -0.0080386531
27 0.0108184210 0.0305147715
28 0.0442418378 0.0108184210
29 -0.0012518616 0.0442418378
30 -0.0140550790 -0.0012518616
31 0.0044473913 -0.0140550790
32 -0.0508766865 0.0044473913
33 0.0445417043 -0.0508766865
34 0.0101176085 0.0445417043
35 0.0033484345 0.0101176085
36 -0.0300056662 0.0033484345
37 -0.0417145343 -0.0300056662
38 -0.0003562172 -0.0417145343
39 -0.0306934507 -0.0003562172
40 -0.0358062769 -0.0306934507
41 0.0238682640 -0.0358062769
42 0.0193396467 0.0238682640
43 0.0007263919 0.0193396467
44 0.0223269267 0.0007263919
45 -0.0599368097 0.0223269267
46 0.0502855712 -0.0599368097
47 0.0093321547 0.0502855712
48 -0.0288024015 0.0093321547
49 0.0009539185 -0.0288024015
50 0.0254014994 0.0009539185
51 0.0336953776 0.0254014994
52 0.0125815109 0.0336953776
53 -0.0065891700 0.0125815109
54 -0.0177613206 -0.0065891700
55 0.0056626977 -0.0177613206
56 -0.0061742266 0.0056626977
57 -0.0683633356 -0.0061742266
58 0.0092046919 -0.0683633356
59 0.0337162279 0.0092046919
60 -0.0062413124 0.0337162279
> 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/7pl5z1227282650.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/freestat/rcomp/tmp/8s6kg1227282650.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/freestat/rcomp/tmp/9g40n1227282650.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
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10t7nm1227282650.ps",horizontal=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/113uvn1227282650.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/128zla1227282650.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/1382ar1227282650.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/14n3nz1227282651.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/15b7vt1227282651.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/168zsr1227282651.tab")
+ }
>
> system("convert tmp/1jiq51227282650.ps tmp/1jiq51227282650.png")
> system("convert tmp/2mwlg1227282650.ps tmp/2mwlg1227282650.png")
> system("convert tmp/3vnu91227282650.ps tmp/3vnu91227282650.png")
> system("convert tmp/43wir1227282650.ps tmp/43wir1227282650.png")
> system("convert tmp/5l5jg1227282650.ps tmp/5l5jg1227282650.png")
> system("convert tmp/6av621227282650.ps tmp/6av621227282650.png")
> system("convert tmp/7pl5z1227282650.ps tmp/7pl5z1227282650.png")
> system("convert tmp/8s6kg1227282650.ps tmp/8s6kg1227282650.png")
> system("convert tmp/9g40n1227282650.ps tmp/9g40n1227282650.png")
> system("convert tmp/10t7nm1227282650.ps tmp/10t7nm1227282650.png")
>
>
> proc.time()
user system elapsed
3.638 2.503 3.968