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(547344,0,554788,0,562325,0,560854,0,555332,0,543599,0,536662,0,542722,0,593530,0,610763,0,612613,0,611324,0,594167,0,595454,0,590865,0,589379,0,584428,0,573100,0,567456,0,569028,0,620735,0,628884,0,628232,0,612117,0,595404,0,597141,0,593408,0,590072,0,579799,0,574205,0,572775,0,572942,0,619567,0,625809,0,619916,0,587625,0,565742,0,557274,0,560576,0,548854,0,531673,0,525919,0,511038,0,498662,0,555362,0,564591,1,541657,1,527070,1,509846,1,514258,1,516922,1,507561,1,492622,1,490243,1,469357,1,477580,1,528379,1,533590,1,517945,1,506174,1,501866,1),dim=c(2,61),dimnames=list(c('Werkloosheid','Kredietcrisis'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Werkloosheid','Kredietcrisis'),1:61))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Werkloosheid Kredietcrisis M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 547344 0 1 0 0 0 0 0 0 0 0 0 0 1
2 554788 0 0 1 0 0 0 0 0 0 0 0 0 2
3 562325 0 0 0 1 0 0 0 0 0 0 0 0 3
4 560854 0 0 0 0 1 0 0 0 0 0 0 0 4
5 555332 0 0 0 0 0 1 0 0 0 0 0 0 5
6 543599 0 0 0 0 0 0 1 0 0 0 0 0 6
7 536662 0 0 0 0 0 0 0 1 0 0 0 0 7
8 542722 0 0 0 0 0 0 0 0 1 0 0 0 8
9 593530 0 0 0 0 0 0 0 0 0 1 0 0 9
10 610763 0 0 0 0 0 0 0 0 0 0 1 0 10
11 612613 0 0 0 0 0 0 0 0 0 0 0 1 11
12 611324 0 0 0 0 0 0 0 0 0 0 0 0 12
13 594167 0 1 0 0 0 0 0 0 0 0 0 0 13
14 595454 0 0 1 0 0 0 0 0 0 0 0 0 14
15 590865 0 0 0 1 0 0 0 0 0 0 0 0 15
16 589379 0 0 0 0 1 0 0 0 0 0 0 0 16
17 584428 0 0 0 0 0 1 0 0 0 0 0 0 17
18 573100 0 0 0 0 0 0 1 0 0 0 0 0 18
19 567456 0 0 0 0 0 0 0 1 0 0 0 0 19
20 569028 0 0 0 0 0 0 0 0 1 0 0 0 20
21 620735 0 0 0 0 0 0 0 0 0 1 0 0 21
22 628884 0 0 0 0 0 0 0 0 0 0 1 0 22
23 628232 0 0 0 0 0 0 0 0 0 0 0 1 23
24 612117 0 0 0 0 0 0 0 0 0 0 0 0 24
25 595404 0 1 0 0 0 0 0 0 0 0 0 0 25
26 597141 0 0 1 0 0 0 0 0 0 0 0 0 26
27 593408 0 0 0 1 0 0 0 0 0 0 0 0 27
28 590072 0 0 0 0 1 0 0 0 0 0 0 0 28
29 579799 0 0 0 0 0 1 0 0 0 0 0 0 29
30 574205 0 0 0 0 0 0 1 0 0 0 0 0 30
31 572775 0 0 0 0 0 0 0 1 0 0 0 0 31
32 572942 0 0 0 0 0 0 0 0 1 0 0 0 32
33 619567 0 0 0 0 0 0 0 0 0 1 0 0 33
34 625809 0 0 0 0 0 0 0 0 0 0 1 0 34
35 619916 0 0 0 0 0 0 0 0 0 0 0 1 35
36 587625 0 0 0 0 0 0 0 0 0 0 0 0 36
37 565742 0 1 0 0 0 0 0 0 0 0 0 0 37
38 557274 0 0 1 0 0 0 0 0 0 0 0 0 38
39 560576 0 0 0 1 0 0 0 0 0 0 0 0 39
40 548854 0 0 0 0 1 0 0 0 0 0 0 0 40
41 531673 0 0 0 0 0 1 0 0 0 0 0 0 41
42 525919 0 0 0 0 0 0 1 0 0 0 0 0 42
43 511038 0 0 0 0 0 0 0 1 0 0 0 0 43
44 498662 0 0 0 0 0 0 0 0 1 0 0 0 44
45 555362 0 0 0 0 0 0 0 0 0 1 0 0 45
46 564591 1 0 0 0 0 0 0 0 0 0 1 0 46
47 541657 1 0 0 0 0 0 0 0 0 0 0 1 47
48 527070 1 0 0 0 0 0 0 0 0 0 0 0 48
49 509846 1 1 0 0 0 0 0 0 0 0 0 0 49
50 514258 1 0 1 0 0 0 0 0 0 0 0 0 50
51 516922 1 0 0 1 0 0 0 0 0 0 0 0 51
52 507561 1 0 0 0 1 0 0 0 0 0 0 0 52
53 492622 1 0 0 0 0 1 0 0 0 0 0 0 53
54 490243 1 0 0 0 0 0 1 0 0 0 0 0 54
55 469357 1 0 0 0 0 0 0 1 0 0 0 0 55
56 477580 1 0 0 0 0 0 0 0 1 0 0 0 56
57 528379 1 0 0 0 0 0 0 0 0 1 0 0 57
58 533590 1 0 0 0 0 0 0 0 0 0 1 0 58
59 517945 1 0 0 0 0 0 0 0 0 0 0 1 59
60 506174 1 0 0 0 0 0 0 0 0 0 0 0 60
61 501866 1 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Kredietcrisis M1 M2 M3
607102.0 -60547.5 -22451.0 -21083.2 -19657.6
M4 M5 M6 M7 M8
-24743.3 -34927.0 -41895.1 -51461.3 -50342.6
M9 M10 M11 t
1374.7 23086.5 14821.1 -389.5
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-40961 -14896 3179 14888 29208
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 607102.0 10870.4 55.849 < 2e-16 ***
Kredietcrisis -60547.5 9315.9 -6.499 4.71e-08 ***
M1 -22451.0 12303.5 -1.825 0.074395 .
M2 -21083.2 12911.2 -1.633 0.109165
M3 -19657.6 12897.3 -1.524 0.134170
M4 -24743.3 12887.5 -1.920 0.060946 .
M5 -34927.0 12881.9 -2.711 0.009329 **
M6 -41895.1 12880.4 -3.253 0.002120 **
M7 -51461.3 12883.0 -3.995 0.000227 ***
M8 -50342.6 12889.8 -3.906 0.000299 ***
M9 1374.7 12900.8 0.107 0.915594
M10 23086.5 12832.2 1.799 0.078425 .
M11 14821.1 12826.0 1.156 0.253703
t -389.5 231.5 -1.683 0.099064 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20280 on 47 degrees of freedom
Multiple R-squared: 0.8111, Adjusted R-squared: 0.7589
F-statistic: 15.53 on 13 and 47 DF, p-value: 7.932e-13
> 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.18323638 0.366472757 0.816763622
[2,] 0.10140112 0.202802235 0.898598883
[3,] 0.04953803 0.099076068 0.950461966
[4,] 0.03064319 0.061286388 0.969356806
[5,] 0.01819359 0.036387176 0.981806412
[6,] 0.03434645 0.068692909 0.965653546
[7,] 0.04834008 0.096680156 0.951659922
[8,] 0.19835399 0.396707983 0.801646008
[9,] 0.22540379 0.450807571 0.774596214
[10,] 0.21735428 0.434708556 0.782645722
[11,] 0.24388949 0.487778990 0.756110505
[12,] 0.24423266 0.488465322 0.755767339
[13,] 0.26269936 0.525398720 0.737300640
[14,] 0.21525341 0.430506811 0.784746595
[15,] 0.20162637 0.403252738 0.798373631
[16,] 0.25153378 0.503067553 0.748466223
[17,] 0.32346848 0.646936960 0.676531520
[18,] 0.39660415 0.793208308 0.603395846
[19,] 0.87431561 0.251368774 0.125684387
[20,] 0.98900480 0.021990394 0.010995197
[21,] 0.99557401 0.008851987 0.004425994
[22,] 0.99673980 0.006520406 0.003260203
[23,] 0.99634869 0.007302622 0.003651311
[24,] 0.99542691 0.009146186 0.004573093
[25,] 0.99346737 0.013065268 0.006532634
[26,] 0.98594732 0.028105365 0.014052682
[27,] 0.98748273 0.025034534 0.012517267
[28,] 0.96622669 0.067546620 0.033773310
> postscript(file="/var/www/html/rcomp/tmp/1gvsn1229988379.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/2dsq31229988379.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/31k9u1229988380.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/4k2kt1229988380.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/5sact1229988380.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 61
Frequency = 1
1 2 3 4 5 6
-36917.5199 -30451.8465 -23951.0465 -19946.8465 -14895.6465 -19271.0465
7 8 9 10 11 12
-16252.4465 -10921.6465 -11441.4465 -15530.7546 -5025.9546 8895.6454
13 14 15 16 17 18
14579.1493 14887.8227 9262.6227 13251.8227 18874.0227 14903.6227
19 20 21 22 23 24
19215.2227 20058.0227 20437.2227 7263.9146 15266.7146 14362.3146
25 26 27 28 29 30
20489.8185 21248.4919 16479.2919 18618.4919 18918.6919 20682.2919
31 32 33 34 35 36
29207.8919 28645.6919 23942.8919 8862.5838 11624.3838 -5456.0162
37 38 39 40 41 42
-4498.5123 -13944.8389 -11679.0389 -17925.8389 -24533.6389 -22930.0389
43 44 45 46 47 48
-27855.4389 -40960.6389 -35588.4389 12865.7935 -1413.4065 -789.8065
49 50 51 52 53 54
4826.6975 8260.3709 9888.1709 6002.3709 1636.5709 6615.1709
55 56 57 58 59 60
-4315.2291 3178.5709 2649.7709 -13461.5373 -20451.7373 -17012.1373
61
1520.3667
> postscript(file="/var/www/html/rcomp/tmp/64rtb1229988380.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -36917.5199 NA
1 -30451.8465 -36917.5199
2 -23951.0465 -30451.8465
3 -19946.8465 -23951.0465
4 -14895.6465 -19946.8465
5 -19271.0465 -14895.6465
6 -16252.4465 -19271.0465
7 -10921.6465 -16252.4465
8 -11441.4465 -10921.6465
9 -15530.7546 -11441.4465
10 -5025.9546 -15530.7546
11 8895.6454 -5025.9546
12 14579.1493 8895.6454
13 14887.8227 14579.1493
14 9262.6227 14887.8227
15 13251.8227 9262.6227
16 18874.0227 13251.8227
17 14903.6227 18874.0227
18 19215.2227 14903.6227
19 20058.0227 19215.2227
20 20437.2227 20058.0227
21 7263.9146 20437.2227
22 15266.7146 7263.9146
23 14362.3146 15266.7146
24 20489.8185 14362.3146
25 21248.4919 20489.8185
26 16479.2919 21248.4919
27 18618.4919 16479.2919
28 18918.6919 18618.4919
29 20682.2919 18918.6919
30 29207.8919 20682.2919
31 28645.6919 29207.8919
32 23942.8919 28645.6919
33 8862.5838 23942.8919
34 11624.3838 8862.5838
35 -5456.0162 11624.3838
36 -4498.5123 -5456.0162
37 -13944.8389 -4498.5123
38 -11679.0389 -13944.8389
39 -17925.8389 -11679.0389
40 -24533.6389 -17925.8389
41 -22930.0389 -24533.6389
42 -27855.4389 -22930.0389
43 -40960.6389 -27855.4389
44 -35588.4389 -40960.6389
45 12865.7935 -35588.4389
46 -1413.4065 12865.7935
47 -789.8065 -1413.4065
48 4826.6975 -789.8065
49 8260.3709 4826.6975
50 9888.1709 8260.3709
51 6002.3709 9888.1709
52 1636.5709 6002.3709
53 6615.1709 1636.5709
54 -4315.2291 6615.1709
55 3178.5709 -4315.2291
56 2649.7709 3178.5709
57 -13461.5373 2649.7709
58 -20451.7373 -13461.5373
59 -17012.1373 -20451.7373
60 1520.3667 -17012.1373
61 NA 1520.3667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -30451.8465 -36917.5199
[2,] -23951.0465 -30451.8465
[3,] -19946.8465 -23951.0465
[4,] -14895.6465 -19946.8465
[5,] -19271.0465 -14895.6465
[6,] -16252.4465 -19271.0465
[7,] -10921.6465 -16252.4465
[8,] -11441.4465 -10921.6465
[9,] -15530.7546 -11441.4465
[10,] -5025.9546 -15530.7546
[11,] 8895.6454 -5025.9546
[12,] 14579.1493 8895.6454
[13,] 14887.8227 14579.1493
[14,] 9262.6227 14887.8227
[15,] 13251.8227 9262.6227
[16,] 18874.0227 13251.8227
[17,] 14903.6227 18874.0227
[18,] 19215.2227 14903.6227
[19,] 20058.0227 19215.2227
[20,] 20437.2227 20058.0227
[21,] 7263.9146 20437.2227
[22,] 15266.7146 7263.9146
[23,] 14362.3146 15266.7146
[24,] 20489.8185 14362.3146
[25,] 21248.4919 20489.8185
[26,] 16479.2919 21248.4919
[27,] 18618.4919 16479.2919
[28,] 18918.6919 18618.4919
[29,] 20682.2919 18918.6919
[30,] 29207.8919 20682.2919
[31,] 28645.6919 29207.8919
[32,] 23942.8919 28645.6919
[33,] 8862.5838 23942.8919
[34,] 11624.3838 8862.5838
[35,] -5456.0162 11624.3838
[36,] -4498.5123 -5456.0162
[37,] -13944.8389 -4498.5123
[38,] -11679.0389 -13944.8389
[39,] -17925.8389 -11679.0389
[40,] -24533.6389 -17925.8389
[41,] -22930.0389 -24533.6389
[42,] -27855.4389 -22930.0389
[43,] -40960.6389 -27855.4389
[44,] -35588.4389 -40960.6389
[45,] 12865.7935 -35588.4389
[46,] -1413.4065 12865.7935
[47,] -789.8065 -1413.4065
[48,] 4826.6975 -789.8065
[49,] 8260.3709 4826.6975
[50,] 9888.1709 8260.3709
[51,] 6002.3709 9888.1709
[52,] 1636.5709 6002.3709
[53,] 6615.1709 1636.5709
[54,] -4315.2291 6615.1709
[55,] 3178.5709 -4315.2291
[56,] 2649.7709 3178.5709
[57,] -13461.5373 2649.7709
[58,] -20451.7373 -13461.5373
[59,] -17012.1373 -20451.7373
[60,] 1520.3667 -17012.1373
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -30451.8465 -36917.5199
2 -23951.0465 -30451.8465
3 -19946.8465 -23951.0465
4 -14895.6465 -19946.8465
5 -19271.0465 -14895.6465
6 -16252.4465 -19271.0465
7 -10921.6465 -16252.4465
8 -11441.4465 -10921.6465
9 -15530.7546 -11441.4465
10 -5025.9546 -15530.7546
11 8895.6454 -5025.9546
12 14579.1493 8895.6454
13 14887.8227 14579.1493
14 9262.6227 14887.8227
15 13251.8227 9262.6227
16 18874.0227 13251.8227
17 14903.6227 18874.0227
18 19215.2227 14903.6227
19 20058.0227 19215.2227
20 20437.2227 20058.0227
21 7263.9146 20437.2227
22 15266.7146 7263.9146
23 14362.3146 15266.7146
24 20489.8185 14362.3146
25 21248.4919 20489.8185
26 16479.2919 21248.4919
27 18618.4919 16479.2919
28 18918.6919 18618.4919
29 20682.2919 18918.6919
30 29207.8919 20682.2919
31 28645.6919 29207.8919
32 23942.8919 28645.6919
33 8862.5838 23942.8919
34 11624.3838 8862.5838
35 -5456.0162 11624.3838
36 -4498.5123 -5456.0162
37 -13944.8389 -4498.5123
38 -11679.0389 -13944.8389
39 -17925.8389 -11679.0389
40 -24533.6389 -17925.8389
41 -22930.0389 -24533.6389
42 -27855.4389 -22930.0389
43 -40960.6389 -27855.4389
44 -35588.4389 -40960.6389
45 12865.7935 -35588.4389
46 -1413.4065 12865.7935
47 -789.8065 -1413.4065
48 4826.6975 -789.8065
49 8260.3709 4826.6975
50 9888.1709 8260.3709
51 6002.3709 9888.1709
52 1636.5709 6002.3709
53 6615.1709 1636.5709
54 -4315.2291 6615.1709
55 3178.5709 -4315.2291
56 2649.7709 3178.5709
57 -13461.5373 2649.7709
58 -20451.7373 -13461.5373
59 -17012.1373 -20451.7373
60 1520.3667 -17012.1373
> 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/76uwy1229988380.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/8w2fb1229988380.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/9r8dn1229988380.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/10tded1229988380.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/115c2h1229988380.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/12kvg21229988380.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/13e3831229988380.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/14jvsy1229988380.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/15vgri1229988380.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/16ziik1229988380.tab")
+ }
>
> system("convert tmp/1gvsn1229988379.ps tmp/1gvsn1229988379.png")
> system("convert tmp/2dsq31229988379.ps tmp/2dsq31229988379.png")
> system("convert tmp/31k9u1229988380.ps tmp/31k9u1229988380.png")
> system("convert tmp/4k2kt1229988380.ps tmp/4k2kt1229988380.png")
> system("convert tmp/5sact1229988380.ps tmp/5sact1229988380.png")
> system("convert tmp/64rtb1229988380.ps tmp/64rtb1229988380.png")
> system("convert tmp/76uwy1229988380.ps tmp/76uwy1229988380.png")
> system("convert tmp/8w2fb1229988380.ps tmp/8w2fb1229988380.png")
> system("convert tmp/9r8dn1229988380.ps tmp/9r8dn1229988380.png")
> system("convert tmp/10tded1229988380.ps tmp/10tded1229988380.png")
>
>
> proc.time()
user system elapsed
2.432 1.608 3.365