R version 2.7.0 (2008-04-22)
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(308347,0,298427,0,289231,0,291975,0,294912,0,293488,0,290555,0,284736,0,281818,0,287854,0,316263,0,325412,0,326011,0,328282,0,317480,0,317539,0,313737,0,312276,0,309391,0,302950,0,300316,0,304035,0,333476,0,337698,0,335932,0,323931,0,313927,0,314485,1,313218,1,309664,1,302963,1,298989,1,298423,1,301631,1,329765,1,335083,1,327616,1,309119,1,295916,1,291413,1,291542,1,284678,1,276475,1,272566,1,264981,1,263290,1,296806,1,303598,1,286994,1,276427,1,266424,1,267153,1,268381,1,262522,1,255542,1,253158,1,243803,1,250741,1,280445,1,285257,1,270976,1,261076,1,255603,1),dim=c(2,63),dimnames=list(c('Vrouwen','Dummy'),1:63))
> y <- array(NA,dim=c(2,63),dimnames=list(c('Vrouwen','Dummy'),1:63))
> 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
Vrouwen Dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 308347 0 1 0 0 0 0 0 0 0 0 0 0 1
2 298427 0 0 1 0 0 0 0 0 0 0 0 0 2
3 289231 0 0 0 1 0 0 0 0 0 0 0 0 3
4 291975 0 0 0 0 1 0 0 0 0 0 0 0 4
5 294912 0 0 0 0 0 1 0 0 0 0 0 0 5
6 293488 0 0 0 0 0 0 1 0 0 0 0 0 6
7 290555 0 0 0 0 0 0 0 1 0 0 0 0 7
8 284736 0 0 0 0 0 0 0 0 1 0 0 0 8
9 281818 0 0 0 0 0 0 0 0 0 1 0 0 9
10 287854 0 0 0 0 0 0 0 0 0 0 1 0 10
11 316263 0 0 0 0 0 0 0 0 0 0 0 1 11
12 325412 0 0 0 0 0 0 0 0 0 0 0 0 12
13 326011 0 1 0 0 0 0 0 0 0 0 0 0 13
14 328282 0 0 1 0 0 0 0 0 0 0 0 0 14
15 317480 0 0 0 1 0 0 0 0 0 0 0 0 15
16 317539 0 0 0 0 1 0 0 0 0 0 0 0 16
17 313737 0 0 0 0 0 1 0 0 0 0 0 0 17
18 312276 0 0 0 0 0 0 1 0 0 0 0 0 18
19 309391 0 0 0 0 0 0 0 1 0 0 0 0 19
20 302950 0 0 0 0 0 0 0 0 1 0 0 0 20
21 300316 0 0 0 0 0 0 0 0 0 1 0 0 21
22 304035 0 0 0 0 0 0 0 0 0 0 1 0 22
23 333476 0 0 0 0 0 0 0 0 0 0 0 1 23
24 337698 0 0 0 0 0 0 0 0 0 0 0 0 24
25 335932 0 1 0 0 0 0 0 0 0 0 0 0 25
26 323931 0 0 1 0 0 0 0 0 0 0 0 0 26
27 313927 0 0 0 1 0 0 0 0 0 0 0 0 27
28 314485 1 0 0 0 1 0 0 0 0 0 0 0 28
29 313218 1 0 0 0 0 1 0 0 0 0 0 0 29
30 309664 1 0 0 0 0 0 1 0 0 0 0 0 30
31 302963 1 0 0 0 0 0 0 1 0 0 0 0 31
32 298989 1 0 0 0 0 0 0 0 1 0 0 0 32
33 298423 1 0 0 0 0 0 0 0 0 1 0 0 33
34 301631 1 0 0 0 0 0 0 0 0 0 1 0 34
35 329765 1 0 0 0 0 0 0 0 0 0 0 1 35
36 335083 1 0 0 0 0 0 0 0 0 0 0 0 36
37 327616 1 1 0 0 0 0 0 0 0 0 0 0 37
38 309119 1 0 1 0 0 0 0 0 0 0 0 0 38
39 295916 1 0 0 1 0 0 0 0 0 0 0 0 39
40 291413 1 0 0 0 1 0 0 0 0 0 0 0 40
41 291542 1 0 0 0 0 1 0 0 0 0 0 0 41
42 284678 1 0 0 0 0 0 1 0 0 0 0 0 42
43 276475 1 0 0 0 0 0 0 1 0 0 0 0 43
44 272566 1 0 0 0 0 0 0 0 1 0 0 0 44
45 264981 1 0 0 0 0 0 0 0 0 1 0 0 45
46 263290 1 0 0 0 0 0 0 0 0 0 1 0 46
47 296806 1 0 0 0 0 0 0 0 0 0 0 1 47
48 303598 1 0 0 0 0 0 0 0 0 0 0 0 48
49 286994 1 1 0 0 0 0 0 0 0 0 0 0 49
50 276427 1 0 1 0 0 0 0 0 0 0 0 0 50
51 266424 1 0 0 1 0 0 0 0 0 0 0 0 51
52 267153 1 0 0 0 1 0 0 0 0 0 0 0 52
53 268381 1 0 0 0 0 1 0 0 0 0 0 0 53
54 262522 1 0 0 0 0 0 1 0 0 0 0 0 54
55 255542 1 0 0 0 0 0 0 1 0 0 0 0 55
56 253158 1 0 0 0 0 0 0 0 1 0 0 0 56
57 243803 1 0 0 0 0 0 0 0 0 1 0 0 57
58 250741 1 0 0 0 0 0 0 0 0 0 1 0 58
59 280445 1 0 0 0 0 0 0 0 0 0 0 1 59
60 285257 1 0 0 0 0 0 0 0 0 0 0 0 60
61 270976 1 1 0 0 0 0 0 0 0 0 0 0 61
62 261076 1 0 1 0 0 0 0 0 0 0 0 0 62
63 255603 1 0 0 1 0 0 0 0 0 0 0 0 63
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dummy M1 M2 M3 M4
349151 9773 -12343 -21067 -29803 -29253
M5 M6 M7 M8 M9 M10
-28364 -31152 -35647 -39108 -42675 -37989
M11 t
-7103 -1045
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-27568 -10637 -1884 13782 25238
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 349150.8 8539.7 40.886 < 2e-16 ***
Dummy 9773.3 8290.3 1.179 0.244138
M1 -12342.5 9789.7 -1.261 0.213364
M2 -21067.0 9782.8 -2.153 0.036229 *
M3 -29802.5 9781.1 -3.047 0.003718 **
M4 -29253.3 10366.3 -2.822 0.006876 **
M5 -28363.7 10329.1 -2.746 0.008413 **
M6 -31151.5 10296.8 -3.025 0.003948 **
M7 -35647.3 10269.4 -3.471 0.001091 **
M8 -39108.2 10246.9 -3.817 0.000380 ***
M9 -42675.2 10229.4 -4.172 0.000123 ***
M10 -37988.6 10216.9 -3.718 0.000516 ***
M11 -7103.2 10209.4 -0.696 0.489871
t -1044.6 226.4 -4.614 2.86e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 16140 on 49 degrees of freedom
Multiple R-squared: 0.6413, Adjusted R-squared: 0.5462
F-statistic: 6.74 on 13 and 49 DF, p-value: 3.613e-07
> 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.5284327 0.943134652 0.4715673260
[2,] 0.4544154 0.908830863 0.5455845685
[3,] 0.3578708 0.715741552 0.6421292241
[4,] 0.2992279 0.598455823 0.7007720884
[5,] 0.2271782 0.454356448 0.7728217760
[6,] 0.2010760 0.402151974 0.7989240131
[7,] 0.1665041 0.333008283 0.8334958587
[8,] 0.2688952 0.537790371 0.7311048144
[9,] 0.4211700 0.842339910 0.5788300451
[10,] 0.7420168 0.515966439 0.2579832194
[11,] 0.8277103 0.344579454 0.1722897271
[12,] 0.7580417 0.483916529 0.2419582645
[13,] 0.6851715 0.629656929 0.3148284645
[14,] 0.5983851 0.803229859 0.4016149296
[15,] 0.5259959 0.948008288 0.4740041440
[16,] 0.4323150 0.864629986 0.5676850071
[17,] 0.3979989 0.795997725 0.6020011373
[18,] 0.3972938 0.794587526 0.6027062368
[19,] 0.3561690 0.712337912 0.6438310439
[20,] 0.3288413 0.657682651 0.6711586747
[21,] 0.7768964 0.446207137 0.2231035685
[22,] 0.9725026 0.054994786 0.0274973929
[23,] 0.9929440 0.014111943 0.0070559714
[24,] 0.9986050 0.002789914 0.0013949568
[25,] 0.9992143 0.001571461 0.0007857305
[26,] 0.9993852 0.001229669 0.0006148347
[27,] 0.9992515 0.001496971 0.0007484854
[28,] 0.9982830 0.003433930 0.0017169651
[29,] 0.9982452 0.003509594 0.0017547972
[30,] 0.9935293 0.012941399 0.0064706993
> postscript(file="/var/www/html/rcomp/tmp/1akr71229463147.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/2bj5b1229463147.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/38lnv1229463147.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/42ua21229463147.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/5slr01229463147.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 = 63
Frequency = 1
1 2 3 4 5 6
-27416.6778 -27567.6778 -26983.5111 -23744.1400 -20652.1400 -18243.7400
7 8 9 10 11 12
-15636.3400 -16949.9400 -15256.3400 -12862.3400 -14294.1400 -11203.7400
13 14 15 16 17 18
2782.3956 14822.3956 13800.5622 14354.9333 10707.9333 13079.3333
19 20 21 22 23 24
15734.7333 13799.1333 15776.7333 15853.7333 15453.9333 13617.3333
25 26 27 28 29 30
25238.4689 23006.4689 22782.6356 14062.6622 12950.6622 13229.0622
31 32 33 34 35 36
12068.4622 12599.8622 16645.4622 16211.4622 14504.6622 13764.0622
37 38 39 40 41 42
19684.1978 10956.1978 7533.3644 3525.7356 3809.7356 778.1356
43 44 45 46 47 48
-1884.4644 -1288.0644 -4261.4644 -9594.4644 -5919.2644 -5185.8644
49 50 51 52 53 54
-8402.7289 -9200.7289 -9423.5622 -8199.1911 -6816.1911 -8842.7911
55 56 57 58 59 60
-10282.3911 -8160.9911 -12904.3911 -9608.3911 -9745.1911 -10991.7911
61 62 63
-11885.6556 -12016.6556 -7709.4889
> postscript(file="/var/www/html/rcomp/tmp/6oyyr1229463147.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 = 63
Frequency = 1
lag(myerror, k = 1) myerror
0 -27416.6778 NA
1 -27567.6778 -27416.6778
2 -26983.5111 -27567.6778
3 -23744.1400 -26983.5111
4 -20652.1400 -23744.1400
5 -18243.7400 -20652.1400
6 -15636.3400 -18243.7400
7 -16949.9400 -15636.3400
8 -15256.3400 -16949.9400
9 -12862.3400 -15256.3400
10 -14294.1400 -12862.3400
11 -11203.7400 -14294.1400
12 2782.3956 -11203.7400
13 14822.3956 2782.3956
14 13800.5622 14822.3956
15 14354.9333 13800.5622
16 10707.9333 14354.9333
17 13079.3333 10707.9333
18 15734.7333 13079.3333
19 13799.1333 15734.7333
20 15776.7333 13799.1333
21 15853.7333 15776.7333
22 15453.9333 15853.7333
23 13617.3333 15453.9333
24 25238.4689 13617.3333
25 23006.4689 25238.4689
26 22782.6356 23006.4689
27 14062.6622 22782.6356
28 12950.6622 14062.6622
29 13229.0622 12950.6622
30 12068.4622 13229.0622
31 12599.8622 12068.4622
32 16645.4622 12599.8622
33 16211.4622 16645.4622
34 14504.6622 16211.4622
35 13764.0622 14504.6622
36 19684.1978 13764.0622
37 10956.1978 19684.1978
38 7533.3644 10956.1978
39 3525.7356 7533.3644
40 3809.7356 3525.7356
41 778.1356 3809.7356
42 -1884.4644 778.1356
43 -1288.0644 -1884.4644
44 -4261.4644 -1288.0644
45 -9594.4644 -4261.4644
46 -5919.2644 -9594.4644
47 -5185.8644 -5919.2644
48 -8402.7289 -5185.8644
49 -9200.7289 -8402.7289
50 -9423.5622 -9200.7289
51 -8199.1911 -9423.5622
52 -6816.1911 -8199.1911
53 -8842.7911 -6816.1911
54 -10282.3911 -8842.7911
55 -8160.9911 -10282.3911
56 -12904.3911 -8160.9911
57 -9608.3911 -12904.3911
58 -9745.1911 -9608.3911
59 -10991.7911 -9745.1911
60 -11885.6556 -10991.7911
61 -12016.6556 -11885.6556
62 -7709.4889 -12016.6556
63 NA -7709.4889
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -27567.6778 -27416.6778
[2,] -26983.5111 -27567.6778
[3,] -23744.1400 -26983.5111
[4,] -20652.1400 -23744.1400
[5,] -18243.7400 -20652.1400
[6,] -15636.3400 -18243.7400
[7,] -16949.9400 -15636.3400
[8,] -15256.3400 -16949.9400
[9,] -12862.3400 -15256.3400
[10,] -14294.1400 -12862.3400
[11,] -11203.7400 -14294.1400
[12,] 2782.3956 -11203.7400
[13,] 14822.3956 2782.3956
[14,] 13800.5622 14822.3956
[15,] 14354.9333 13800.5622
[16,] 10707.9333 14354.9333
[17,] 13079.3333 10707.9333
[18,] 15734.7333 13079.3333
[19,] 13799.1333 15734.7333
[20,] 15776.7333 13799.1333
[21,] 15853.7333 15776.7333
[22,] 15453.9333 15853.7333
[23,] 13617.3333 15453.9333
[24,] 25238.4689 13617.3333
[25,] 23006.4689 25238.4689
[26,] 22782.6356 23006.4689
[27,] 14062.6622 22782.6356
[28,] 12950.6622 14062.6622
[29,] 13229.0622 12950.6622
[30,] 12068.4622 13229.0622
[31,] 12599.8622 12068.4622
[32,] 16645.4622 12599.8622
[33,] 16211.4622 16645.4622
[34,] 14504.6622 16211.4622
[35,] 13764.0622 14504.6622
[36,] 19684.1978 13764.0622
[37,] 10956.1978 19684.1978
[38,] 7533.3644 10956.1978
[39,] 3525.7356 7533.3644
[40,] 3809.7356 3525.7356
[41,] 778.1356 3809.7356
[42,] -1884.4644 778.1356
[43,] -1288.0644 -1884.4644
[44,] -4261.4644 -1288.0644
[45,] -9594.4644 -4261.4644
[46,] -5919.2644 -9594.4644
[47,] -5185.8644 -5919.2644
[48,] -8402.7289 -5185.8644
[49,] -9200.7289 -8402.7289
[50,] -9423.5622 -9200.7289
[51,] -8199.1911 -9423.5622
[52,] -6816.1911 -8199.1911
[53,] -8842.7911 -6816.1911
[54,] -10282.3911 -8842.7911
[55,] -8160.9911 -10282.3911
[56,] -12904.3911 -8160.9911
[57,] -9608.3911 -12904.3911
[58,] -9745.1911 -9608.3911
[59,] -10991.7911 -9745.1911
[60,] -11885.6556 -10991.7911
[61,] -12016.6556 -11885.6556
[62,] -7709.4889 -12016.6556
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -27567.6778 -27416.6778
2 -26983.5111 -27567.6778
3 -23744.1400 -26983.5111
4 -20652.1400 -23744.1400
5 -18243.7400 -20652.1400
6 -15636.3400 -18243.7400
7 -16949.9400 -15636.3400
8 -15256.3400 -16949.9400
9 -12862.3400 -15256.3400
10 -14294.1400 -12862.3400
11 -11203.7400 -14294.1400
12 2782.3956 -11203.7400
13 14822.3956 2782.3956
14 13800.5622 14822.3956
15 14354.9333 13800.5622
16 10707.9333 14354.9333
17 13079.3333 10707.9333
18 15734.7333 13079.3333
19 13799.1333 15734.7333
20 15776.7333 13799.1333
21 15853.7333 15776.7333
22 15453.9333 15853.7333
23 13617.3333 15453.9333
24 25238.4689 13617.3333
25 23006.4689 25238.4689
26 22782.6356 23006.4689
27 14062.6622 22782.6356
28 12950.6622 14062.6622
29 13229.0622 12950.6622
30 12068.4622 13229.0622
31 12599.8622 12068.4622
32 16645.4622 12599.8622
33 16211.4622 16645.4622
34 14504.6622 16211.4622
35 13764.0622 14504.6622
36 19684.1978 13764.0622
37 10956.1978 19684.1978
38 7533.3644 10956.1978
39 3525.7356 7533.3644
40 3809.7356 3525.7356
41 778.1356 3809.7356
42 -1884.4644 778.1356
43 -1288.0644 -1884.4644
44 -4261.4644 -1288.0644
45 -9594.4644 -4261.4644
46 -5919.2644 -9594.4644
47 -5185.8644 -5919.2644
48 -8402.7289 -5185.8644
49 -9200.7289 -8402.7289
50 -9423.5622 -9200.7289
51 -8199.1911 -9423.5622
52 -6816.1911 -8199.1911
53 -8842.7911 -6816.1911
54 -10282.3911 -8842.7911
55 -8160.9911 -10282.3911
56 -12904.3911 -8160.9911
57 -9608.3911 -12904.3911
58 -9745.1911 -9608.3911
59 -10991.7911 -9745.1911
60 -11885.6556 -10991.7911
61 -12016.6556 -11885.6556
62 -7709.4889 -12016.6556
> 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/7021o1229463147.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/8yx9g1229463147.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/9p5lq1229463147.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/10zmzc1229463147.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/11j9cg1229463147.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/12jvr91229463147.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/13ks361229463147.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/14ivjc1229463147.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/15g5pw1229463148.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/16m9vu1229463148.tab")
+ }
>
> system("convert tmp/1akr71229463147.ps tmp/1akr71229463147.png")
> system("convert tmp/2bj5b1229463147.ps tmp/2bj5b1229463147.png")
> system("convert tmp/38lnv1229463147.ps tmp/38lnv1229463147.png")
> system("convert tmp/42ua21229463147.ps tmp/42ua21229463147.png")
> system("convert tmp/5slr01229463147.ps tmp/5slr01229463147.png")
> system("convert tmp/6oyyr1229463147.ps tmp/6oyyr1229463147.png")
> system("convert tmp/7021o1229463147.ps tmp/7021o1229463147.png")
> system("convert tmp/8yx9g1229463147.ps tmp/8yx9g1229463147.png")
> system("convert tmp/9p5lq1229463147.ps tmp/9p5lq1229463147.png")
> system("convert tmp/10zmzc1229463147.ps tmp/10zmzc1229463147.png")
>
>
> proc.time()
user system elapsed
5.708 2.716 6.132