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 = '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)
> 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
1 325412 285351 1 0 0 0 0 0 0 0 0 0 0
2 326011 286602 0 1 0 0 0 0 0 0 0 0 0
3 328282 283042 0 0 1 0 0 0 0 0 0 0 0
4 317480 276687 0 0 0 1 0 0 0 0 0 0 0
5 317539 277915 0 0 0 0 1 0 0 0 0 0 0
6 313737 277128 0 0 0 0 0 1 0 0 0 0 0
7 312276 277103 0 0 0 0 0 0 1 0 0 0 0
8 309391 275037 0 0 0 0 0 0 0 1 0 0 0
9 302950 270150 0 0 0 0 0 0 0 0 1 0 0
10 300316 267140 0 0 0 0 0 0 0 0 0 1 0
11 304035 264993 0 0 0 0 0 0 0 0 0 0 1
12 333476 287259 0 0 0 0 0 0 0 0 0 0 0
13 337698 291186 1 0 0 0 0 0 0 0 0 0 0
14 335932 292300 0 1 0 0 0 0 0 0 0 0 0
15 323931 288186 0 0 1 0 0 0 0 0 0 0 0
16 313927 281477 0 0 0 1 0 0 0 0 0 0 0
17 314485 282656 0 0 0 0 1 0 0 0 0 0 0
18 313218 280190 0 0 0 0 0 1 0 0 0 0 0
19 309664 280408 0 0 0 0 0 0 1 0 0 0 0
20 302963 276836 0 0 0 0 0 0 0 1 0 0 0
21 298989 275216 0 0 0 0 0 0 0 0 1 0 0
22 298423 274352 0 0 0 0 0 0 0 0 0 1 0
23 310631 271311 0 0 0 0 0 0 0 0 0 0 1
24 329765 289802 0 0 0 0 0 0 0 0 0 0 0
25 335083 290726 1 0 0 0 0 0 0 0 0 0 0
26 327616 292300 0 1 0 0 0 0 0 0 0 0 0
27 309119 278506 0 0 1 0 0 0 0 0 0 0 0
28 295916 269826 0 0 0 1 0 0 0 0 0 0 0
29 291413 265861 0 0 0 0 1 0 0 0 0 0 0
30 291542 269034 0 0 0 0 0 1 0 0 0 0 0
31 284678 264176 0 0 0 0 0 0 1 0 0 0 0
32 276475 255198 0 0 0 0 0 0 0 1 0 0 0
33 272566 253353 0 0 0 0 0 0 0 0 1 0 0
34 264981 246057 0 0 0 0 0 0 0 0 0 1 0
35 263290 235372 0 0 0 0 0 0 0 0 0 0 1
36 296806 258556 0 0 0 0 0 0 0 0 0 0 0
37 303598 260993 1 0 0 0 0 0 0 0 0 0 0
38 286994 254663 0 1 0 0 0 0 0 0 0 0 0
39 276427 250643 0 0 1 0 0 0 0 0 0 0 0
40 266424 243422 0 0 0 1 0 0 0 0 0 0 0
41 267153 247105 0 0 0 0 1 0 0 0 0 0 0
42 268381 248541 0 0 0 0 0 1 0 0 0 0 0
43 262522 245039 0 0 0 0 0 0 1 0 0 0 0
44 255542 237080 0 0 0 0 0 0 0 1 0 0 0
45 253158 237085 0 0 0 0 0 0 0 0 1 0 0
46 243803 225554 0 0 0 0 0 0 0 0 0 1 0
47 250741 226839 0 0 0 0 0 0 0 0 0 0 1
48 280445 247934 0 0 0 0 0 0 0 0 0 0 0
49 285257 248333 1 0 0 0 0 0 0 0 0 0 0
50 270976 246969 0 1 0 0 0 0 0 0 0 0 0
51 261076 245098 0 0 1 0 0 0 0 0 0 0 0
52 255603 246263 0 0 0 1 0 0 0 0 0 0 0
53 260376 255765 0 0 0 0 1 0 0 0 0 0 0
54 263903 264319 0 0 0 0 0 1 0 0 0 0 0
55 264291 268347 0 0 0 0 0 0 1 0 0 0 0
56 263276 273046 0 0 0 0 0 0 0 1 0 0 0
57 262572 273963 0 0 0 0 0 0 0 0 1 0 0
58 256167 267430 0 0 0 0 0 0 0 0 0 1 0
59 264221 271993 0 0 0 0 0 0 0 0 0 0 1
60 293860 292710 0 0 0 0 0 0 0 0 0 0 0
61 300713 295881 1 0 0 0 0 0 0 0 0 0 0
62 287224 293299 0 1 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) Werkl_mannen M1 M2 M3
8149.593 1.085 3965.830 -3722.619 -421.224
M4 M5 M6 M7 M8
-4284.166 -6484.635 -8672.624 -11244.244 -12521.015
M9 M10 M11
-14390.715 -13354.404 -5332.853
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-35509 -1323 3744 8391 16008
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.150e+03 3.456e+04 0.236 0.815
Werkl_mannen 1.085e+00 1.229e-01 8.831 1.05e-11 ***
M1 3.966e+03 9.586e+03 0.414 0.681
M2 -3.723e+03 9.581e+03 -0.389 0.699
M3 -4.212e+02 1.003e+04 -0.042 0.967
M4 -4.284e+03 1.011e+04 -0.424 0.673
M5 -6.485e+03 1.007e+04 -0.644 0.523
M6 -8.673e+03 1.004e+04 -0.863 0.392
M7 -1.124e+04 1.005e+04 -1.118 0.269
M8 -1.252e+04 1.011e+04 -1.239 0.221
M9 -1.439e+04 1.014e+04 -1.420 0.162
M10 -1.335e+04 1.028e+04 -1.300 0.200
M11 -5.333e+03 1.033e+04 -0.516 0.608
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15820 on 49 degrees of freedom
Multiple R-squared: 0.7072, Adjusted R-squared: 0.6355
F-statistic: 9.862 on 12 and 49 DF, p-value: 2.138e-09
> 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.0666004508 0.1332009015 0.9333995492
[2,] 0.0283276872 0.0566553744 0.9716723128
[3,] 0.0096712585 0.0193425170 0.9903287415
[4,] 0.0036587456 0.0073174911 0.9963412544
[5,] 0.0019775071 0.0039550141 0.9980224929
[6,] 0.0008854286 0.0017708573 0.9991145714
[7,] 0.0003559048 0.0007118097 0.9996440952
[8,] 0.0002944270 0.0005888540 0.9997055730
[9,] 0.0001761266 0.0003522531 0.9998238734
[10,] 0.0001175430 0.0002350861 0.9998824570
[11,] 0.0001704123 0.0003408245 0.9998295877
[12,] 0.0024050102 0.0048100203 0.9975949898
[13,] 0.0116965643 0.0233931285 0.9883034357
[14,] 0.0250731367 0.0501462734 0.9749268633
[15,] 0.0574842246 0.1149684492 0.9425157754
[16,] 0.0900255093 0.1800510186 0.9099744907
[17,] 0.1151237330 0.2302474659 0.8848762670
[18,] 0.1414856937 0.2829713875 0.8585143063
[19,] 0.1928784170 0.3857568341 0.8071215830
[20,] 0.1807951616 0.3615903233 0.8192048384
[21,] 0.2546473669 0.5092947337 0.7453526331
[22,] 0.4282237076 0.8564474153 0.5717762924
[23,] 0.5710236555 0.8579526890 0.4289763445
[24,] 0.7640466666 0.4719066667 0.2359533334
[25,] 0.8937183515 0.2125632970 0.1062816485
[26,] 0.9555808737 0.0888382526 0.0444191263
[27,] 0.9951880588 0.0096238825 0.0048119412
[28,] 0.9986494032 0.0027011936 0.0013505968
[29,] 0.9990800956 0.0018398088 0.0009199044
[30,] 0.9985979883 0.0028040233 0.0014020117
[31,] 0.9931880766 0.0136238468 0.0068119234
> postscript(file="/var/www/html/rcomp/tmp/172ay1258482009.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/2nmnl1258482009.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/346on1258482009.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/4mv901258482009.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/52y491258482009.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
3615.92460 10545.71051 13378.84816 13336.63161 14263.39834 13503.48936
7 8 9 10 11 12
14641.24087 15275.16316 16007.53989 15603.86809 13631.37533 13575.07387
13 14 15 16 17 18
9569.41969 14282.88652 3445.25943 4585.22569 6064.17027 9661.41652
19 20 21 22 23 24
8442.44931 6894.77647 6548.60161 5883.95714 13370.68878 7104.25211
25 26 27 28 29 30
7453.64030 5966.88652 -861.40252 -781.38447 1219.14884 92.60157
31 32 33 34 35 36
1072.42527 3889.67986 3852.68898 3149.45096 5032.92681 8055.35467
37 38 39 40 41 42
8236.74116 6190.89977 -3314.74197 -1618.12146 -2685.67342 -828.12026
43 44 45 46 47 48
-314.91210 2619.46031 2099.73438 4222.58175 1744.46913 3222.00971
49 50 51 52 53 54
3635.16055 -1477.09290 -12647.96310 -15522.35136 -18861.04403 -22429.38718
55 56 57 58 59 60
-23841.20335 -28679.07981 -28508.56486 -28859.85794 -33779.46004 -31956.69035
61 62
-32510.88631 -35509.29041
> postscript(file="/var/www/html/rcomp/tmp/6rycb1258482009.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 3615.92460 NA
1 10545.71051 3615.92460
2 13378.84816 10545.71051
3 13336.63161 13378.84816
4 14263.39834 13336.63161
5 13503.48936 14263.39834
6 14641.24087 13503.48936
7 15275.16316 14641.24087
8 16007.53989 15275.16316
9 15603.86809 16007.53989
10 13631.37533 15603.86809
11 13575.07387 13631.37533
12 9569.41969 13575.07387
13 14282.88652 9569.41969
14 3445.25943 14282.88652
15 4585.22569 3445.25943
16 6064.17027 4585.22569
17 9661.41652 6064.17027
18 8442.44931 9661.41652
19 6894.77647 8442.44931
20 6548.60161 6894.77647
21 5883.95714 6548.60161
22 13370.68878 5883.95714
23 7104.25211 13370.68878
24 7453.64030 7104.25211
25 5966.88652 7453.64030
26 -861.40252 5966.88652
27 -781.38447 -861.40252
28 1219.14884 -781.38447
29 92.60157 1219.14884
30 1072.42527 92.60157
31 3889.67986 1072.42527
32 3852.68898 3889.67986
33 3149.45096 3852.68898
34 5032.92681 3149.45096
35 8055.35467 5032.92681
36 8236.74116 8055.35467
37 6190.89977 8236.74116
38 -3314.74197 6190.89977
39 -1618.12146 -3314.74197
40 -2685.67342 -1618.12146
41 -828.12026 -2685.67342
42 -314.91210 -828.12026
43 2619.46031 -314.91210
44 2099.73438 2619.46031
45 4222.58175 2099.73438
46 1744.46913 4222.58175
47 3222.00971 1744.46913
48 3635.16055 3222.00971
49 -1477.09290 3635.16055
50 -12647.96310 -1477.09290
51 -15522.35136 -12647.96310
52 -18861.04403 -15522.35136
53 -22429.38718 -18861.04403
54 -23841.20335 -22429.38718
55 -28679.07981 -23841.20335
56 -28508.56486 -28679.07981
57 -28859.85794 -28508.56486
58 -33779.46004 -28859.85794
59 -31956.69035 -33779.46004
60 -32510.88631 -31956.69035
61 -35509.29041 -32510.88631
62 NA -35509.29041
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 10545.71051 3615.92460
[2,] 13378.84816 10545.71051
[3,] 13336.63161 13378.84816
[4,] 14263.39834 13336.63161
[5,] 13503.48936 14263.39834
[6,] 14641.24087 13503.48936
[7,] 15275.16316 14641.24087
[8,] 16007.53989 15275.16316
[9,] 15603.86809 16007.53989
[10,] 13631.37533 15603.86809
[11,] 13575.07387 13631.37533
[12,] 9569.41969 13575.07387
[13,] 14282.88652 9569.41969
[14,] 3445.25943 14282.88652
[15,] 4585.22569 3445.25943
[16,] 6064.17027 4585.22569
[17,] 9661.41652 6064.17027
[18,] 8442.44931 9661.41652
[19,] 6894.77647 8442.44931
[20,] 6548.60161 6894.77647
[21,] 5883.95714 6548.60161
[22,] 13370.68878 5883.95714
[23,] 7104.25211 13370.68878
[24,] 7453.64030 7104.25211
[25,] 5966.88652 7453.64030
[26,] -861.40252 5966.88652
[27,] -781.38447 -861.40252
[28,] 1219.14884 -781.38447
[29,] 92.60157 1219.14884
[30,] 1072.42527 92.60157
[31,] 3889.67986 1072.42527
[32,] 3852.68898 3889.67986
[33,] 3149.45096 3852.68898
[34,] 5032.92681 3149.45096
[35,] 8055.35467 5032.92681
[36,] 8236.74116 8055.35467
[37,] 6190.89977 8236.74116
[38,] -3314.74197 6190.89977
[39,] -1618.12146 -3314.74197
[40,] -2685.67342 -1618.12146
[41,] -828.12026 -2685.67342
[42,] -314.91210 -828.12026
[43,] 2619.46031 -314.91210
[44,] 2099.73438 2619.46031
[45,] 4222.58175 2099.73438
[46,] 1744.46913 4222.58175
[47,] 3222.00971 1744.46913
[48,] 3635.16055 3222.00971
[49,] -1477.09290 3635.16055
[50,] -12647.96310 -1477.09290
[51,] -15522.35136 -12647.96310
[52,] -18861.04403 -15522.35136
[53,] -22429.38718 -18861.04403
[54,] -23841.20335 -22429.38718
[55,] -28679.07981 -23841.20335
[56,] -28508.56486 -28679.07981
[57,] -28859.85794 -28508.56486
[58,] -33779.46004 -28859.85794
[59,] -31956.69035 -33779.46004
[60,] -32510.88631 -31956.69035
[61,] -35509.29041 -32510.88631
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 10545.71051 3615.92460
2 13378.84816 10545.71051
3 13336.63161 13378.84816
4 14263.39834 13336.63161
5 13503.48936 14263.39834
6 14641.24087 13503.48936
7 15275.16316 14641.24087
8 16007.53989 15275.16316
9 15603.86809 16007.53989
10 13631.37533 15603.86809
11 13575.07387 13631.37533
12 9569.41969 13575.07387
13 14282.88652 9569.41969
14 3445.25943 14282.88652
15 4585.22569 3445.25943
16 6064.17027 4585.22569
17 9661.41652 6064.17027
18 8442.44931 9661.41652
19 6894.77647 8442.44931
20 6548.60161 6894.77647
21 5883.95714 6548.60161
22 13370.68878 5883.95714
23 7104.25211 13370.68878
24 7453.64030 7104.25211
25 5966.88652 7453.64030
26 -861.40252 5966.88652
27 -781.38447 -861.40252
28 1219.14884 -781.38447
29 92.60157 1219.14884
30 1072.42527 92.60157
31 3889.67986 1072.42527
32 3852.68898 3889.67986
33 3149.45096 3852.68898
34 5032.92681 3149.45096
35 8055.35467 5032.92681
36 8236.74116 8055.35467
37 6190.89977 8236.74116
38 -3314.74197 6190.89977
39 -1618.12146 -3314.74197
40 -2685.67342 -1618.12146
41 -828.12026 -2685.67342
42 -314.91210 -828.12026
43 2619.46031 -314.91210
44 2099.73438 2619.46031
45 4222.58175 2099.73438
46 1744.46913 4222.58175
47 3222.00971 1744.46913
48 3635.16055 3222.00971
49 -1477.09290 3635.16055
50 -12647.96310 -1477.09290
51 -15522.35136 -12647.96310
52 -18861.04403 -15522.35136
53 -22429.38718 -18861.04403
54 -23841.20335 -22429.38718
55 -28679.07981 -23841.20335
56 -28508.56486 -28679.07981
57 -28859.85794 -28508.56486
58 -33779.46004 -28859.85794
59 -31956.69035 -33779.46004
60 -32510.88631 -31956.69035
61 -35509.29041 -32510.88631
> 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/7ha301258482009.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/8ew1c1258482009.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/9nhmj1258482009.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/106jk01258482009.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/11w88b1258482009.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/12yonz1258482009.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/13v2en1258482009.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/14p40l1258482009.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/15eddx1258482009.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/169zlf1258482009.tab")
+ }
>
> system("convert tmp/172ay1258482009.ps tmp/172ay1258482009.png")
> system("convert tmp/2nmnl1258482009.ps tmp/2nmnl1258482009.png")
> system("convert tmp/346on1258482009.ps tmp/346on1258482009.png")
> system("convert tmp/4mv901258482009.ps tmp/4mv901258482009.png")
> system("convert tmp/52y491258482009.ps tmp/52y491258482009.png")
> system("convert tmp/6rycb1258482009.ps tmp/6rycb1258482009.png")
> system("convert tmp/7ha301258482009.ps tmp/7ha301258482009.png")
> system("convert tmp/8ew1c1258482009.ps tmp/8ew1c1258482009.png")
> system("convert tmp/9nhmj1258482009.ps tmp/9nhmj1258482009.png")
> system("convert tmp/106jk01258482009.ps tmp/106jk01258482009.png")
>
>
> proc.time()
user system elapsed
2.444 1.569 3.216