R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(325412,285351,326011,286602,328282,283042,317480,276687,317539,277915,313737,277128,312276,277103,309391,275037,302950,270150,300316,267140,304035,264993,333476,287259,337698,291186,335932,292300,323931,288186,313927,281477,314485,282656,313218,280190,309664,280408,302963,276836,298989,275216,298423,274352,310631,271311,329765,289802,335083,290726,327616,292300,309119,278506,295916,269826,291413,265861,291542,269034,284678,264176,276475,255198,272566,253353,264981,246057,263290,235372,296806,258556,303598,260993,286994,254663,276427,250643,266424,243422,267153,247105,268381,248541,262522,245039,255542,237080,253158,237085,243803,225554,250741,226839,280445,247934,285257,248333,270976,246969,261076,245098,255603,246263,260376,255765,263903,264319,264291,268347,263276,273046,262572,273963,256167,267430,264221,271993,293860,292710,300713,295881,287224,293299),dim=c(2,62),dimnames=list(c('Werkl_vrouwen','Werkl_mannen'),1:62))
> y <- array(NA,dim=c(2,62),dimnames=list(c('Werkl_vrouwen','Werkl_mannen'),1:62))
> 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
Werkl_vrouwen Werkl_mannen M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 325412 285351 1 0 0 0 0 0 0 0 0 0 0 1
2 326011 286602 0 1 0 0 0 0 0 0 0 0 0 2
3 328282 283042 0 0 1 0 0 0 0 0 0 0 0 3
4 317480 276687 0 0 0 1 0 0 0 0 0 0 0 4
5 317539 277915 0 0 0 0 1 0 0 0 0 0 0 5
6 313737 277128 0 0 0 0 0 1 0 0 0 0 0 6
7 312276 277103 0 0 0 0 0 0 1 0 0 0 0 7
8 309391 275037 0 0 0 0 0 0 0 1 0 0 0 8
9 302950 270150 0 0 0 0 0 0 0 0 1 0 0 9
10 300316 267140 0 0 0 0 0 0 0 0 0 1 0 10
11 304035 264993 0 0 0 0 0 0 0 0 0 0 1 11
12 333476 287259 0 0 0 0 0 0 0 0 0 0 0 12
13 337698 291186 1 0 0 0 0 0 0 0 0 0 0 13
14 335932 292300 0 1 0 0 0 0 0 0 0 0 0 14
15 323931 288186 0 0 1 0 0 0 0 0 0 0 0 15
16 313927 281477 0 0 0 1 0 0 0 0 0 0 0 16
17 314485 282656 0 0 0 0 1 0 0 0 0 0 0 17
18 313218 280190 0 0 0 0 0 1 0 0 0 0 0 18
19 309664 280408 0 0 0 0 0 0 1 0 0 0 0 19
20 302963 276836 0 0 0 0 0 0 0 1 0 0 0 20
21 298989 275216 0 0 0 0 0 0 0 0 1 0 0 21
22 298423 274352 0 0 0 0 0 0 0 0 0 1 0 22
23 310631 271311 0 0 0 0 0 0 0 0 0 0 1 23
24 329765 289802 0 0 0 0 0 0 0 0 0 0 0 24
25 335083 290726 1 0 0 0 0 0 0 0 0 0 0 25
26 327616 292300 0 1 0 0 0 0 0 0 0 0 0 26
27 309119 278506 0 0 1 0 0 0 0 0 0 0 0 27
28 295916 269826 0 0 0 1 0 0 0 0 0 0 0 28
29 291413 265861 0 0 0 0 1 0 0 0 0 0 0 29
30 291542 269034 0 0 0 0 0 1 0 0 0 0 0 30
31 284678 264176 0 0 0 0 0 0 1 0 0 0 0 31
32 276475 255198 0 0 0 0 0 0 0 1 0 0 0 32
33 272566 253353 0 0 0 0 0 0 0 0 1 0 0 33
34 264981 246057 0 0 0 0 0 0 0 0 0 1 0 34
35 263290 235372 0 0 0 0 0 0 0 0 0 0 1 35
36 296806 258556 0 0 0 0 0 0 0 0 0 0 0 36
37 303598 260993 1 0 0 0 0 0 0 0 0 0 0 37
38 286994 254663 0 1 0 0 0 0 0 0 0 0 0 38
39 276427 250643 0 0 1 0 0 0 0 0 0 0 0 39
40 266424 243422 0 0 0 1 0 0 0 0 0 0 0 40
41 267153 247105 0 0 0 0 1 0 0 0 0 0 0 41
42 268381 248541 0 0 0 0 0 1 0 0 0 0 0 42
43 262522 245039 0 0 0 0 0 0 1 0 0 0 0 43
44 255542 237080 0 0 0 0 0 0 0 1 0 0 0 44
45 253158 237085 0 0 0 0 0 0 0 0 1 0 0 45
46 243803 225554 0 0 0 0 0 0 0 0 0 1 0 46
47 250741 226839 0 0 0 0 0 0 0 0 0 0 1 47
48 280445 247934 0 0 0 0 0 0 0 0 0 0 0 48
49 285257 248333 1 0 0 0 0 0 0 0 0 0 0 49
50 270976 246969 0 1 0 0 0 0 0 0 0 0 0 50
51 261076 245098 0 0 1 0 0 0 0 0 0 0 0 51
52 255603 246263 0 0 0 1 0 0 0 0 0 0 0 52
53 260376 255765 0 0 0 0 1 0 0 0 0 0 0 53
54 263903 264319 0 0 0 0 0 1 0 0 0 0 0 54
55 264291 268347 0 0 0 0 0 0 1 0 0 0 0 55
56 263276 273046 0 0 0 0 0 0 0 1 0 0 0 56
57 262572 273963 0 0 0 0 0 0 0 0 1 0 0 57
58 256167 267430 0 0 0 0 0 0 0 0 0 1 0 58
59 264221 271993 0 0 0 0 0 0 0 0 0 0 1 59
60 293860 292710 0 0 0 0 0 0 0 0 0 0 0 60
61 300713 295881 1 0 0 0 0 0 0 0 0 0 0 61
62 287224 293299 0 1 0 0 0 0 0 0 0 0 0 62
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Werkl_mannen M1 M2 M3
1.609e+05 6.429e-01 1.206e+03 -6.088e+03 -1.089e+04
M4 M5 M6 M7 M8
-1.635e+04 -1.667e+04 -1.712e+04 -1.919e+04 -2.119e+04
M9 M10 M11 t
-2.286e+04 -2.354e+04 -1.555e+04 -8.609e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-19288.29 -2102.78 32.24 1890.78 10652.38
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.609e+05 1.248e+04 12.892 < 2e-16 ***
Werkl_mannen 6.429e-01 4.204e-02 15.292 < 2e-16 ***
M1 1.206e+03 2.890e+03 0.417 0.678222
M2 -6.088e+03 2.888e+03 -2.108 0.040257 *
M3 -1.089e+04 3.058e+03 -3.562 0.000843 ***
M4 -1.635e+04 3.092e+03 -5.289 2.99e-06 ***
M5 -1.667e+04 3.067e+03 -5.433 1.82e-06 ***
M6 -1.712e+04 3.049e+03 -5.613 9.71e-07 ***
M7 -1.919e+04 3.049e+03 -6.294 8.90e-08 ***
M8 -2.119e+04 3.069e+03 -6.904 1.03e-08 ***
M9 -2.286e+04 3.076e+03 -7.429 1.62e-09 ***
M10 -2.354e+04 3.129e+03 -7.525 1.16e-09 ***
M11 -1.555e+04 3.147e+03 -4.941 9.84e-06 ***
t -8.609e+02 3.881e+01 -22.183 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4764 on 48 degrees of freedom
Multiple R-squared: 0.974, Adjusted R-squared: 0.9669
F-statistic: 138.2 on 13 and 48 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.7161608 5.676785e-01 2.838392e-01
[2,] 0.9976573 4.685410e-03 2.342705e-03
[3,] 0.9943683 1.126345e-02 5.631723e-03
[4,] 0.9888326 2.233487e-02 1.116743e-02
[5,] 0.9915181 1.696373e-02 8.481866e-03
[6,] 0.9946111 1.077771e-02 5.388855e-03
[7,] 0.9999443 1.114663e-04 5.573314e-05
[8,] 0.9998467 3.065426e-04 1.532713e-04
[9,] 0.9999486 1.028165e-04 5.140824e-05
[10,] 0.9999929 1.427176e-05 7.135879e-06
[11,] 0.9999994 1.239129e-06 6.195645e-07
[12,] 0.9999998 3.766060e-07 1.883030e-07
[13,] 0.9999999 1.920694e-07 9.603472e-08
[14,] 0.9999999 2.253091e-07 1.126546e-07
[15,] 0.9999996 7.825755e-07 3.912877e-07
[16,] 0.9999990 1.998213e-06 9.991065e-07
[17,] 0.9999965 6.944764e-06 3.472382e-06
[18,] 0.9999921 1.577935e-05 7.889675e-06
[19,] 0.9999998 4.454675e-07 2.227337e-07
[20,] 0.9999997 5.164768e-07 2.582384e-07
[21,] 0.9999998 4.129769e-07 2.064885e-07
[22,] 0.9999998 4.503836e-07 2.251918e-07
[23,] 0.9999986 2.754559e-06 1.377280e-06
[24,] 0.9999909 1.823464e-05 9.117321e-06
[25,] 0.9999736 5.273150e-05 2.636575e-05
[26,] 0.9999025 1.949543e-04 9.747715e-05
[27,] 0.9996207 7.585321e-04 3.792661e-04
[28,] 0.9982763 3.447485e-03 1.723742e-03
[29,] 0.9895647 2.087061e-02 1.043530e-02
> postscript(file="/var/www/html/rcomp/tmp/12anl1258482727.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/2ardb1258482727.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/3jwxy1258482727.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/4gcp61258482727.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/56kde1258482727.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 = 62
Frequency = 1
1 2 3 4 5 6
-19288.29083 -11337.98261 -1112.99846 -1506.86505 -1065.50807 -3050.22076
7 8 9 10 11 12
-1557.36906 -255.88672 -1027.86443 -176.74274 -2212.04595 -1774.97577
13 14 15 16 17 18
-423.32588 5250.06445 1559.23465 2190.96708 3162.82788 4792.60323
19 20 21 22 23 24
4036.22170 2489.96430 2084.51761 3623.90057 10652.38134 3209.53154
25 26 27 28 29 30
7587.91554 7264.55573 3301.33774 2001.29526 1219.41424 619.67853
31 32 33 34 35 36
-183.16441 244.28537 48.49907 -1295.75938 -3251.68267 670.17309
37 38 39 40 41 42
5549.79628 1171.19553 -1146.07018 -184.15497 -651.20367 964.83890
43 44 45 46 47 48
295.17595 1290.47486 1430.25862 1038.83034 15.97385 1468.92104
49 50 51 52 53 54
5678.84597 430.42950 -2601.50375 -2501.24232 -2665.53039 -3326.89989
55 56 57 58 59 60
-2590.86417 -3768.83781 -2535.41087 -3190.22880 -5204.62656 -3573.64990
61 62
895.05892 -2778.26260
> postscript(file="/var/www/html/rcomp/tmp/6q7ck1258482727.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 = 62
Frequency = 1
lag(myerror, k = 1) myerror
0 -19288.29083 NA
1 -11337.98261 -19288.29083
2 -1112.99846 -11337.98261
3 -1506.86505 -1112.99846
4 -1065.50807 -1506.86505
5 -3050.22076 -1065.50807
6 -1557.36906 -3050.22076
7 -255.88672 -1557.36906
8 -1027.86443 -255.88672
9 -176.74274 -1027.86443
10 -2212.04595 -176.74274
11 -1774.97577 -2212.04595
12 -423.32588 -1774.97577
13 5250.06445 -423.32588
14 1559.23465 5250.06445
15 2190.96708 1559.23465
16 3162.82788 2190.96708
17 4792.60323 3162.82788
18 4036.22170 4792.60323
19 2489.96430 4036.22170
20 2084.51761 2489.96430
21 3623.90057 2084.51761
22 10652.38134 3623.90057
23 3209.53154 10652.38134
24 7587.91554 3209.53154
25 7264.55573 7587.91554
26 3301.33774 7264.55573
27 2001.29526 3301.33774
28 1219.41424 2001.29526
29 619.67853 1219.41424
30 -183.16441 619.67853
31 244.28537 -183.16441
32 48.49907 244.28537
33 -1295.75938 48.49907
34 -3251.68267 -1295.75938
35 670.17309 -3251.68267
36 5549.79628 670.17309
37 1171.19553 5549.79628
38 -1146.07018 1171.19553
39 -184.15497 -1146.07018
40 -651.20367 -184.15497
41 964.83890 -651.20367
42 295.17595 964.83890
43 1290.47486 295.17595
44 1430.25862 1290.47486
45 1038.83034 1430.25862
46 15.97385 1038.83034
47 1468.92104 15.97385
48 5678.84597 1468.92104
49 430.42950 5678.84597
50 -2601.50375 430.42950
51 -2501.24232 -2601.50375
52 -2665.53039 -2501.24232
53 -3326.89989 -2665.53039
54 -2590.86417 -3326.89989
55 -3768.83781 -2590.86417
56 -2535.41087 -3768.83781
57 -3190.22880 -2535.41087
58 -5204.62656 -3190.22880
59 -3573.64990 -5204.62656
60 895.05892 -3573.64990
61 -2778.26260 895.05892
62 NA -2778.26260
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -11337.98261 -19288.29083
[2,] -1112.99846 -11337.98261
[3,] -1506.86505 -1112.99846
[4,] -1065.50807 -1506.86505
[5,] -3050.22076 -1065.50807
[6,] -1557.36906 -3050.22076
[7,] -255.88672 -1557.36906
[8,] -1027.86443 -255.88672
[9,] -176.74274 -1027.86443
[10,] -2212.04595 -176.74274
[11,] -1774.97577 -2212.04595
[12,] -423.32588 -1774.97577
[13,] 5250.06445 -423.32588
[14,] 1559.23465 5250.06445
[15,] 2190.96708 1559.23465
[16,] 3162.82788 2190.96708
[17,] 4792.60323 3162.82788
[18,] 4036.22170 4792.60323
[19,] 2489.96430 4036.22170
[20,] 2084.51761 2489.96430
[21,] 3623.90057 2084.51761
[22,] 10652.38134 3623.90057
[23,] 3209.53154 10652.38134
[24,] 7587.91554 3209.53154
[25,] 7264.55573 7587.91554
[26,] 3301.33774 7264.55573
[27,] 2001.29526 3301.33774
[28,] 1219.41424 2001.29526
[29,] 619.67853 1219.41424
[30,] -183.16441 619.67853
[31,] 244.28537 -183.16441
[32,] 48.49907 244.28537
[33,] -1295.75938 48.49907
[34,] -3251.68267 -1295.75938
[35,] 670.17309 -3251.68267
[36,] 5549.79628 670.17309
[37,] 1171.19553 5549.79628
[38,] -1146.07018 1171.19553
[39,] -184.15497 -1146.07018
[40,] -651.20367 -184.15497
[41,] 964.83890 -651.20367
[42,] 295.17595 964.83890
[43,] 1290.47486 295.17595
[44,] 1430.25862 1290.47486
[45,] 1038.83034 1430.25862
[46,] 15.97385 1038.83034
[47,] 1468.92104 15.97385
[48,] 5678.84597 1468.92104
[49,] 430.42950 5678.84597
[50,] -2601.50375 430.42950
[51,] -2501.24232 -2601.50375
[52,] -2665.53039 -2501.24232
[53,] -3326.89989 -2665.53039
[54,] -2590.86417 -3326.89989
[55,] -3768.83781 -2590.86417
[56,] -2535.41087 -3768.83781
[57,] -3190.22880 -2535.41087
[58,] -5204.62656 -3190.22880
[59,] -3573.64990 -5204.62656
[60,] 895.05892 -3573.64990
[61,] -2778.26260 895.05892
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -11337.98261 -19288.29083
2 -1112.99846 -11337.98261
3 -1506.86505 -1112.99846
4 -1065.50807 -1506.86505
5 -3050.22076 -1065.50807
6 -1557.36906 -3050.22076
7 -255.88672 -1557.36906
8 -1027.86443 -255.88672
9 -176.74274 -1027.86443
10 -2212.04595 -176.74274
11 -1774.97577 -2212.04595
12 -423.32588 -1774.97577
13 5250.06445 -423.32588
14 1559.23465 5250.06445
15 2190.96708 1559.23465
16 3162.82788 2190.96708
17 4792.60323 3162.82788
18 4036.22170 4792.60323
19 2489.96430 4036.22170
20 2084.51761 2489.96430
21 3623.90057 2084.51761
22 10652.38134 3623.90057
23 3209.53154 10652.38134
24 7587.91554 3209.53154
25 7264.55573 7587.91554
26 3301.33774 7264.55573
27 2001.29526 3301.33774
28 1219.41424 2001.29526
29 619.67853 1219.41424
30 -183.16441 619.67853
31 244.28537 -183.16441
32 48.49907 244.28537
33 -1295.75938 48.49907
34 -3251.68267 -1295.75938
35 670.17309 -3251.68267
36 5549.79628 670.17309
37 1171.19553 5549.79628
38 -1146.07018 1171.19553
39 -184.15497 -1146.07018
40 -651.20367 -184.15497
41 964.83890 -651.20367
42 295.17595 964.83890
43 1290.47486 295.17595
44 1430.25862 1290.47486
45 1038.83034 1430.25862
46 15.97385 1038.83034
47 1468.92104 15.97385
48 5678.84597 1468.92104
49 430.42950 5678.84597
50 -2601.50375 430.42950
51 -2501.24232 -2601.50375
52 -2665.53039 -2501.24232
53 -3326.89989 -2665.53039
54 -2590.86417 -3326.89989
55 -3768.83781 -2590.86417
56 -2535.41087 -3768.83781
57 -3190.22880 -2535.41087
58 -5204.62656 -3190.22880
59 -3573.64990 -5204.62656
60 895.05892 -3573.64990
61 -2778.26260 895.05892
> 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/7v6pa1258482727.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/89g791258482727.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/9budw1258482727.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/rcomp/tmp/10kxqy1258482727.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> 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/11r7sm1258482727.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/12xmft1258482728.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/13vuxf1258482728.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/14jh5q1258482728.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/rcomp/tmp/15uysc1258482728.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/rcomp/tmp/16w3nm1258482728.tab")
+ }
>
> system("convert tmp/12anl1258482727.ps tmp/12anl1258482727.png")
> system("convert tmp/2ardb1258482727.ps tmp/2ardb1258482727.png")
> system("convert tmp/3jwxy1258482727.ps tmp/3jwxy1258482727.png")
> system("convert tmp/4gcp61258482727.ps tmp/4gcp61258482727.png")
> system("convert tmp/56kde1258482727.ps tmp/56kde1258482727.png")
> system("convert tmp/6q7ck1258482727.ps tmp/6q7ck1258482727.png")
> system("convert tmp/7v6pa1258482727.ps tmp/7v6pa1258482727.png")
> system("convert tmp/89g791258482727.ps tmp/89g791258482727.png")
> system("convert tmp/9budw1258482727.ps tmp/9budw1258482727.png")
> system("convert tmp/10kxqy1258482727.ps tmp/10kxqy1258482727.png")
>
>
> proc.time()
user system elapsed
2.445 1.581 4.379