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(121148,0,114624,0,109822,0,112081,0,113534,0,112110,0,109826,0,107423,0,105540,0,108573,0,128591,0,139145,0,129700,0,132828,0,126868,0,128390,0,126830,0,124105,0,122323,0,119296,0,116822,0,119224,0,139357,0,144322,0,133676,0,128283,0,121640,0,122877,1,117284,1,116463,1,112685,1,113235,1,111692,1,113152,1,129889,1,131153,1,123770,1,112516,1,105940,1,104320,1,103582,1,99064,1,94989,1,92241,1,89752,1,90610,1,109456,1,110213,1,97694,1,91844,1,87572,1,89812,1,89050,1,85990,1,85070,1,83277,1,79586,1,84215,1,99708,1,100698,1,90861,1,86700,1),dim=c(2,62),dimnames=list(c('werkl.vrouwen','Wetswijziging'),1:62))
> y <- array(NA,dim=c(2,62),dimnames=list(c('werkl.vrouwen','Wetswijziging'),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 Wetswijziging M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 121148 0 1 0 0 0 0 0 0 0 0 0 0 1
2 114624 0 0 1 0 0 0 0 0 0 0 0 0 2
3 109822 0 0 0 1 0 0 0 0 0 0 0 0 3
4 112081 0 0 0 0 1 0 0 0 0 0 0 0 4
5 113534 0 0 0 0 0 1 0 0 0 0 0 0 5
6 112110 0 0 0 0 0 0 1 0 0 0 0 0 6
7 109826 0 0 0 0 0 0 0 1 0 0 0 0 7
8 107423 0 0 0 0 0 0 0 0 1 0 0 0 8
9 105540 0 0 0 0 0 0 0 0 0 1 0 0 9
10 108573 0 0 0 0 0 0 0 0 0 0 1 0 10
11 128591 0 0 0 0 0 0 0 0 0 0 0 1 11
12 139145 0 0 0 0 0 0 0 0 0 0 0 0 12
13 129700 0 1 0 0 0 0 0 0 0 0 0 0 13
14 132828 0 0 1 0 0 0 0 0 0 0 0 0 14
15 126868 0 0 0 1 0 0 0 0 0 0 0 0 15
16 128390 0 0 0 0 1 0 0 0 0 0 0 0 16
17 126830 0 0 0 0 0 1 0 0 0 0 0 0 17
18 124105 0 0 0 0 0 0 1 0 0 0 0 0 18
19 122323 0 0 0 0 0 0 0 1 0 0 0 0 19
20 119296 0 0 0 0 0 0 0 0 1 0 0 0 20
21 116822 0 0 0 0 0 0 0 0 0 1 0 0 21
22 119224 0 0 0 0 0 0 0 0 0 0 1 0 22
23 139357 0 0 0 0 0 0 0 0 0 0 0 1 23
24 144322 0 0 0 0 0 0 0 0 0 0 0 0 24
25 133676 0 1 0 0 0 0 0 0 0 0 0 0 25
26 128283 0 0 1 0 0 0 0 0 0 0 0 0 26
27 121640 0 0 0 1 0 0 0 0 0 0 0 0 27
28 122877 1 0 0 0 1 0 0 0 0 0 0 0 28
29 117284 1 0 0 0 0 1 0 0 0 0 0 0 29
30 116463 1 0 0 0 0 0 1 0 0 0 0 0 30
31 112685 1 0 0 0 0 0 0 1 0 0 0 0 31
32 113235 1 0 0 0 0 0 0 0 1 0 0 0 32
33 111692 1 0 0 0 0 0 0 0 0 1 0 0 33
34 113152 1 0 0 0 0 0 0 0 0 0 1 0 34
35 129889 1 0 0 0 0 0 0 0 0 0 0 1 35
36 131153 1 0 0 0 0 0 0 0 0 0 0 0 36
37 123770 1 1 0 0 0 0 0 0 0 0 0 0 37
38 112516 1 0 1 0 0 0 0 0 0 0 0 0 38
39 105940 1 0 0 1 0 0 0 0 0 0 0 0 39
40 104320 1 0 0 0 1 0 0 0 0 0 0 0 40
41 103582 1 0 0 0 0 1 0 0 0 0 0 0 41
42 99064 1 0 0 0 0 0 1 0 0 0 0 0 42
43 94989 1 0 0 0 0 0 0 1 0 0 0 0 43
44 92241 1 0 0 0 0 0 0 0 1 0 0 0 44
45 89752 1 0 0 0 0 0 0 0 0 1 0 0 45
46 90610 1 0 0 0 0 0 0 0 0 0 1 0 46
47 109456 1 0 0 0 0 0 0 0 0 0 0 1 47
48 110213 1 0 0 0 0 0 0 0 0 0 0 0 48
49 97694 1 1 0 0 0 0 0 0 0 0 0 0 49
50 91844 1 0 1 0 0 0 0 0 0 0 0 0 50
51 87572 1 0 0 1 0 0 0 0 0 0 0 0 51
52 89812 1 0 0 0 1 0 0 0 0 0 0 0 52
53 89050 1 0 0 0 0 1 0 0 0 0 0 0 53
54 85990 1 0 0 0 0 0 1 0 0 0 0 0 54
55 85070 1 0 0 0 0 0 0 1 0 0 0 0 55
56 83277 1 0 0 0 0 0 0 0 1 0 0 0 56
57 79586 1 0 0 0 0 0 0 0 0 1 0 0 57
58 84215 1 0 0 0 0 0 0 0 0 0 1 0 58
59 99708 1 0 0 0 0 0 0 0 0 0 0 1 59
60 100698 1 0 0 0 0 0 0 0 0 0 0 0 60
61 90861 1 1 0 0 0 0 0 0 0 0 0 0 61
62 86700 1 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) Wetswijziging M1 M2 M3
148463.3 -623.5 -12219.1 -16589.7 -20608.3
M4 M5 M6 M7 M8
-18717.5 -19519.1 -21390.3 -23319.7 -24565.5
M9 M10 M11 t
-26343.1 -23228.2 -4344.4 -638.4
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-16118 -6513 -1868 8701 13392
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 148463.3 5313.9 27.939 < 2e-16 ***
Wetswijziging -623.5 5145.1 -0.121 0.904058
M1 -12219.1 6031.2 -2.026 0.048344 *
M2 -16589.7 6026.5 -2.753 0.008317 **
M3 -20608.3 6320.7 -3.260 0.002049 **
M4 -18717.5 6390.6 -2.929 0.005190 **
M5 -19519.1 6366.5 -3.066 0.003558 **
M6 -21390.3 6345.6 -3.371 0.001487 **
M7 -23319.7 6327.9 -3.685 0.000580 ***
M8 -24565.5 6313.4 -3.891 0.000307 ***
M9 -26343.1 6302.0 -4.180 0.000123 ***
M10 -23228.2 6293.9 -3.691 0.000571 ***
M11 -4344.4 6289.1 -0.691 0.493025
t -638.4 142.9 -4.467 4.82e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9941 on 48 degrees of freedom
Multiple R-squared: 0.7059, Adjusted R-squared: 0.6262
F-statistic: 8.861 on 13 and 48 DF, p-value: 7.68e-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.5152811 0.9694377995 0.4847188998
[2,] 0.3898748 0.7797496196 0.6101251902
[3,] 0.2663990 0.5327980749 0.7336009626
[4,] 0.1922279 0.3844558023 0.8077720989
[5,] 0.1489379 0.2978757118 0.8510621441
[6,] 0.1323566 0.2647132811 0.8676433595
[7,] 0.1014551 0.2029102525 0.8985448737
[8,] 0.2159597 0.4319193660 0.7840403170
[9,] 0.4910315 0.9820630918 0.5089684541
[10,] 0.7396780 0.5206439884 0.2603219942
[11,] 0.8365462 0.3269075892 0.1634537946
[12,] 0.7831354 0.4337292227 0.2168646113
[13,] 0.7312839 0.5374322154 0.2687161077
[14,] 0.6504176 0.6991648161 0.3495824081
[15,] 0.5654673 0.8690653092 0.4345326546
[16,] 0.4908105 0.9816209387 0.5091895306
[17,] 0.4546959 0.9093917485 0.5453041258
[18,] 0.4022687 0.8045373164 0.5977313418
[19,] 0.3720116 0.7440231748 0.6279884126
[20,] 0.4797929 0.9595857187 0.5202071407
[21,] 0.8017131 0.3965738409 0.1982869205
[22,] 0.9403643 0.1192713558 0.0596356779
[23,] 0.9890804 0.0218392779 0.0109196390
[24,] 0.9971317 0.0057365319 0.0028682660
[25,] 0.9990712 0.0018575438 0.0009287719
[26,] 0.9996051 0.0007897887 0.0003948944
[27,] 0.9989914 0.0020171474 0.0010085737
[28,] 0.9962154 0.0075691231 0.0037845616
[29,] 0.9893612 0.0212776625 0.0106388312
> postscript(file="/var/www/html/rcomp/tmp/1uisy1229538047.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/2hqfb1229538047.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/3021m1229538047.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/4o3dx1229538047.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/5y49i1229538047.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
-14457.7857 -15972.7857 -16117.8286 -15111.1214 -12218.1214 -11132.5214
7 8 9 10 11 12
-10848.7214 -11367.5214 -10834.5214 -10277.9214 -8505.3214 -1657.3214
13 14 15 16 17 18
1755.2357 9892.2357 8589.1929 8858.9000 8738.9000 8523.5000
19 20 21 22 23 24
9309.3000 8166.5000 8108.5000 8034.1000 9921.7000 11180.7000
25 26 27 28 29 30
13392.2571 13008.2571 11022.2143 11630.3857 7477.3857 9165.9857
31 32 33 34 35 36
7955.7857 10389.9857 11262.9857 10246.5857 8738.1857 6296.1857
37 38 39 40 41 42
11770.7429 5525.7429 3606.7000 734.4071 1436.4071 -571.9929
43 44 45 46 47 48
-2079.1929 -2942.9929 -3015.9929 -4634.3929 -4033.7929 -6982.7929
49 50 51 52 53 54
-6644.2357 -7485.2357 -7100.2786 -6112.5714 -5434.5714 -5984.9714
55 56 57 58 59 60
-4337.1714 -4245.9714 -5520.9714 -3368.3714 -6120.7714 -8836.7714
61 62
-5816.2143 -4968.2143
> postscript(file="/var/www/html/rcomp/tmp/6el7p1229538047.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 -14457.7857 NA
1 -15972.7857 -14457.7857
2 -16117.8286 -15972.7857
3 -15111.1214 -16117.8286
4 -12218.1214 -15111.1214
5 -11132.5214 -12218.1214
6 -10848.7214 -11132.5214
7 -11367.5214 -10848.7214
8 -10834.5214 -11367.5214
9 -10277.9214 -10834.5214
10 -8505.3214 -10277.9214
11 -1657.3214 -8505.3214
12 1755.2357 -1657.3214
13 9892.2357 1755.2357
14 8589.1929 9892.2357
15 8858.9000 8589.1929
16 8738.9000 8858.9000
17 8523.5000 8738.9000
18 9309.3000 8523.5000
19 8166.5000 9309.3000
20 8108.5000 8166.5000
21 8034.1000 8108.5000
22 9921.7000 8034.1000
23 11180.7000 9921.7000
24 13392.2571 11180.7000
25 13008.2571 13392.2571
26 11022.2143 13008.2571
27 11630.3857 11022.2143
28 7477.3857 11630.3857
29 9165.9857 7477.3857
30 7955.7857 9165.9857
31 10389.9857 7955.7857
32 11262.9857 10389.9857
33 10246.5857 11262.9857
34 8738.1857 10246.5857
35 6296.1857 8738.1857
36 11770.7429 6296.1857
37 5525.7429 11770.7429
38 3606.7000 5525.7429
39 734.4071 3606.7000
40 1436.4071 734.4071
41 -571.9929 1436.4071
42 -2079.1929 -571.9929
43 -2942.9929 -2079.1929
44 -3015.9929 -2942.9929
45 -4634.3929 -3015.9929
46 -4033.7929 -4634.3929
47 -6982.7929 -4033.7929
48 -6644.2357 -6982.7929
49 -7485.2357 -6644.2357
50 -7100.2786 -7485.2357
51 -6112.5714 -7100.2786
52 -5434.5714 -6112.5714
53 -5984.9714 -5434.5714
54 -4337.1714 -5984.9714
55 -4245.9714 -4337.1714
56 -5520.9714 -4245.9714
57 -3368.3714 -5520.9714
58 -6120.7714 -3368.3714
59 -8836.7714 -6120.7714
60 -5816.2143 -8836.7714
61 -4968.2143 -5816.2143
62 NA -4968.2143
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -15972.7857 -14457.7857
[2,] -16117.8286 -15972.7857
[3,] -15111.1214 -16117.8286
[4,] -12218.1214 -15111.1214
[5,] -11132.5214 -12218.1214
[6,] -10848.7214 -11132.5214
[7,] -11367.5214 -10848.7214
[8,] -10834.5214 -11367.5214
[9,] -10277.9214 -10834.5214
[10,] -8505.3214 -10277.9214
[11,] -1657.3214 -8505.3214
[12,] 1755.2357 -1657.3214
[13,] 9892.2357 1755.2357
[14,] 8589.1929 9892.2357
[15,] 8858.9000 8589.1929
[16,] 8738.9000 8858.9000
[17,] 8523.5000 8738.9000
[18,] 9309.3000 8523.5000
[19,] 8166.5000 9309.3000
[20,] 8108.5000 8166.5000
[21,] 8034.1000 8108.5000
[22,] 9921.7000 8034.1000
[23,] 11180.7000 9921.7000
[24,] 13392.2571 11180.7000
[25,] 13008.2571 13392.2571
[26,] 11022.2143 13008.2571
[27,] 11630.3857 11022.2143
[28,] 7477.3857 11630.3857
[29,] 9165.9857 7477.3857
[30,] 7955.7857 9165.9857
[31,] 10389.9857 7955.7857
[32,] 11262.9857 10389.9857
[33,] 10246.5857 11262.9857
[34,] 8738.1857 10246.5857
[35,] 6296.1857 8738.1857
[36,] 11770.7429 6296.1857
[37,] 5525.7429 11770.7429
[38,] 3606.7000 5525.7429
[39,] 734.4071 3606.7000
[40,] 1436.4071 734.4071
[41,] -571.9929 1436.4071
[42,] -2079.1929 -571.9929
[43,] -2942.9929 -2079.1929
[44,] -3015.9929 -2942.9929
[45,] -4634.3929 -3015.9929
[46,] -4033.7929 -4634.3929
[47,] -6982.7929 -4033.7929
[48,] -6644.2357 -6982.7929
[49,] -7485.2357 -6644.2357
[50,] -7100.2786 -7485.2357
[51,] -6112.5714 -7100.2786
[52,] -5434.5714 -6112.5714
[53,] -5984.9714 -5434.5714
[54,] -4337.1714 -5984.9714
[55,] -4245.9714 -4337.1714
[56,] -5520.9714 -4245.9714
[57,] -3368.3714 -5520.9714
[58,] -6120.7714 -3368.3714
[59,] -8836.7714 -6120.7714
[60,] -5816.2143 -8836.7714
[61,] -4968.2143 -5816.2143
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -15972.7857 -14457.7857
2 -16117.8286 -15972.7857
3 -15111.1214 -16117.8286
4 -12218.1214 -15111.1214
5 -11132.5214 -12218.1214
6 -10848.7214 -11132.5214
7 -11367.5214 -10848.7214
8 -10834.5214 -11367.5214
9 -10277.9214 -10834.5214
10 -8505.3214 -10277.9214
11 -1657.3214 -8505.3214
12 1755.2357 -1657.3214
13 9892.2357 1755.2357
14 8589.1929 9892.2357
15 8858.9000 8589.1929
16 8738.9000 8858.9000
17 8523.5000 8738.9000
18 9309.3000 8523.5000
19 8166.5000 9309.3000
20 8108.5000 8166.5000
21 8034.1000 8108.5000
22 9921.7000 8034.1000
23 11180.7000 9921.7000
24 13392.2571 11180.7000
25 13008.2571 13392.2571
26 11022.2143 13008.2571
27 11630.3857 11022.2143
28 7477.3857 11630.3857
29 9165.9857 7477.3857
30 7955.7857 9165.9857
31 10389.9857 7955.7857
32 11262.9857 10389.9857
33 10246.5857 11262.9857
34 8738.1857 10246.5857
35 6296.1857 8738.1857
36 11770.7429 6296.1857
37 5525.7429 11770.7429
38 3606.7000 5525.7429
39 734.4071 3606.7000
40 1436.4071 734.4071
41 -571.9929 1436.4071
42 -2079.1929 -571.9929
43 -2942.9929 -2079.1929
44 -3015.9929 -2942.9929
45 -4634.3929 -3015.9929
46 -4033.7929 -4634.3929
47 -6982.7929 -4033.7929
48 -6644.2357 -6982.7929
49 -7485.2357 -6644.2357
50 -7100.2786 -7485.2357
51 -6112.5714 -7100.2786
52 -5434.5714 -6112.5714
53 -5984.9714 -5434.5714
54 -4337.1714 -5984.9714
55 -4245.9714 -4337.1714
56 -5520.9714 -4245.9714
57 -3368.3714 -5520.9714
58 -6120.7714 -3368.3714
59 -8836.7714 -6120.7714
60 -5816.2143 -8836.7714
61 -4968.2143 -5816.2143
> 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/7qfdg1229538047.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/8gi401229538047.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/99cqf1229538047.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/10ow5h1229538048.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/1151x31229538048.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/12k1dv1229538048.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/13senl1229538048.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/14rw9c1229538048.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/15e01r1229538048.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/16exey1229538048.tab")
+ }
>
> system("convert tmp/1uisy1229538047.ps tmp/1uisy1229538047.png")
> system("convert tmp/2hqfb1229538047.ps tmp/2hqfb1229538047.png")
> system("convert tmp/3021m1229538047.ps tmp/3021m1229538047.png")
> system("convert tmp/4o3dx1229538047.ps tmp/4o3dx1229538047.png")
> system("convert tmp/5y49i1229538047.ps tmp/5y49i1229538047.png")
> system("convert tmp/6el7p1229538047.ps tmp/6el7p1229538047.png")
> system("convert tmp/7qfdg1229538047.ps tmp/7qfdg1229538047.png")
> system("convert tmp/8gi401229538047.ps tmp/8gi401229538047.png")
> system("convert tmp/99cqf1229538047.ps tmp/99cqf1229538047.png")
> system("convert tmp/10ow5h1229538048.ps tmp/10ow5h1229538048.png")
>
>
> proc.time()
user system elapsed
2.741 1.764 4.707