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(95.2,0,95.00,0,94.00,0,92.2,0,91.00,0,91.2,0,103.4,1,105.00,0,104.6,0,103.8,0,101.8,0,102.4,0,103.8,0,103.4,0,102.00,0,101.8,0,100.2,0,101.4,0,113.8,0,116.00,0,115.6,0,113.00,0,109.4,0,111.00,0,112.4,0,112.2,0,111.00,0,108.8,0,107.4,0,108.6,0,118.8,0,122.2,1,122.6,0,122.2,0,118.8,0,119.00,0,118.2,0,117.8,0,116.8,0,114.6,0,113.4,0,113.8,0,124.2,0,125.8,0,125.6,0,122.4,0,119.00,0,119.4,0,118.6,0,118.00,0,116.00,0,114.8,0,114.6,0,114.6,0,124.00,0,125.2,0,124.00,0,117.6,1,113.2,0,111.4,0,112.2,0,109.8,0,106.4,0,105.2,0,102.2,0,99.8,0,111.00,0,113.00,0,108.4,0,105.4,0,102.00,0,102.8,0),dim=c(2,72),dimnames=list(c('Werkloosheid','Dumivariabele'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Werkloosheid','Dumivariabele'),1:72))
> 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 = 'Do not include Seasonal 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 Dumivariabele
1 95.2 0
2 95.0 0
3 94.0 0
4 92.2 0
5 91.0 0
6 91.2 0
7 103.4 1
8 105.0 0
9 104.6 0
10 103.8 0
11 101.8 0
12 102.4 0
13 103.8 0
14 103.4 0
15 102.0 0
16 101.8 0
17 100.2 0
18 101.4 0
19 113.8 0
20 116.0 0
21 115.6 0
22 113.0 0
23 109.4 0
24 111.0 0
25 112.4 0
26 112.2 0
27 111.0 0
28 108.8 0
29 107.4 0
30 108.6 0
31 118.8 0
32 122.2 1
33 122.6 0
34 122.2 0
35 118.8 0
36 119.0 0
37 118.2 0
38 117.8 0
39 116.8 0
40 114.6 0
41 113.4 0
42 113.8 0
43 124.2 0
44 125.8 0
45 125.6 0
46 122.4 0
47 119.0 0
48 119.4 0
49 118.6 0
50 118.0 0
51 116.0 0
52 114.8 0
53 114.6 0
54 114.6 0
55 124.0 0
56 125.2 0
57 124.0 0
58 117.6 1
59 113.2 0
60 111.4 0
61 112.2 0
62 109.8 0
63 106.4 0
64 105.2 0
65 102.2 0
66 99.8 0
67 111.0 0
68 113.0 0
69 108.4 0
70 105.4 0
71 102.0 0
72 102.8 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dumivariabele
110.623 3.777
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-19.623 -6.923 1.577 7.227 15.177
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 110.623 1.079 102.518 <2e-16 ***
Dumivariabele 3.777 5.286 0.714 0.477
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.963 on 70 degrees of freedom
Multiple R-squared: 0.007239, Adjusted R-squared: -0.006943
F-statistic: 0.5104 on 1 and 70 DF, p-value: 0.4773
> 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.02510841 0.050216825 0.974891588
[2,] 0.01144932 0.022898630 0.988550685
[3,] 0.00303265 0.006065300 0.996967350
[4,] 0.14280555 0.285611102 0.857194449
[5,] 0.23904469 0.478089375 0.760955313
[6,] 0.26410151 0.528203019 0.735898490
[7,] 0.23549082 0.470981645 0.764509177
[8,] 0.21386814 0.427736290 0.786131855
[9,] 0.20764227 0.415284550 0.792357725
[10,] 0.19261947 0.385238933 0.807380533
[11,] 0.16932433 0.338648659 0.830675670
[12,] 0.15037947 0.300758937 0.849620531
[13,] 0.13829677 0.276593546 0.861703227
[14,] 0.13080125 0.261602495 0.869198752
[15,] 0.34395291 0.687905816 0.656047092
[16,] 0.59466535 0.810669304 0.405334652
[17,] 0.73194905 0.536101907 0.268050953
[18,] 0.76697269 0.466054614 0.233027307
[19,] 0.75352649 0.492947014 0.246473507
[20,] 0.74757730 0.504845401 0.252422700
[21,] 0.74885345 0.502293107 0.251146553
[22,] 0.74165938 0.516681242 0.258340621
[23,] 0.72091862 0.558162757 0.279081379
[24,] 0.69053980 0.618920397 0.309460199
[25,] 0.66176428 0.676471439 0.338235720
[26,] 0.63096942 0.738061160 0.369030580
[27,] 0.70196291 0.596074187 0.298037093
[28,] 0.74219664 0.515606720 0.257803360
[29,] 0.83877251 0.322454989 0.161227494
[30,] 0.89153915 0.216921692 0.108460846
[31,] 0.89742972 0.205140557 0.102570278
[32,] 0.90115194 0.197696120 0.098848060
[33,] 0.89691211 0.206175785 0.103087893
[34,] 0.88798450 0.224030992 0.112015496
[35,] 0.87095238 0.258095234 0.129047617
[36,] 0.84071316 0.318573680 0.159286840
[37,] 0.80178736 0.396425278 0.198212639
[38,] 0.75786653 0.484266930 0.242133465
[39,] 0.81606343 0.367873146 0.183936573
[40,] 0.88449546 0.231009090 0.115504545
[41,] 0.93195344 0.136093112 0.068046556
[42,] 0.94500046 0.109999089 0.054999545
[43,] 0.93974241 0.120515175 0.060257588
[44,] 0.93673281 0.126534377 0.063267188
[45,] 0.93012100 0.139757994 0.069878997
[46,] 0.92047999 0.159040022 0.079520011
[47,] 0.89931705 0.201365908 0.100682954
[48,] 0.86815724 0.263685527 0.131842764
[49,] 0.82981828 0.340363435 0.170181718
[50,] 0.78508319 0.429833627 0.214916814
[51,] 0.87490462 0.250190762 0.125095381
[52,] 0.96620053 0.067598944 0.033799472
[53,] 0.99767326 0.004653472 0.002326736
[54,] 0.99490070 0.010198599 0.005099300
[55,] 0.99419548 0.011609045 0.005804523
[56,] 0.99141033 0.017179340 0.008589670
[57,] 0.99026792 0.019464155 0.009732077
[58,] 0.98389311 0.032213788 0.016106894
[59,] 0.96519651 0.069606986 0.034803493
[60,] 0.92869325 0.142613499 0.071306750
[61,] 0.88479139 0.230417215 0.115208607
[62,] 0.88259319 0.234813616 0.117406808
[63,] 0.80956368 0.380872637 0.190436319
> postscript(file="/var/www/html/rcomp/tmp/19r9z1228489851.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/2rtez1228489851.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/33lcu1228489851.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/458ux1228489851.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/5pecg1228489851.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 = 72
Frequency = 1
1 2 3 4 5 6
-15.4231884 -15.6231884 -16.6231884 -18.4231884 -19.6231884 -19.4231884
7 8 9 10 11 12
-11.0000000 -5.6231884 -6.0231884 -6.8231884 -8.8231884 -8.2231884
13 14 15 16 17 18
-6.8231884 -7.2231884 -8.6231884 -8.8231884 -10.4231884 -9.2231884
19 20 21 22 23 24
3.1768116 5.3768116 4.9768116 2.3768116 -1.2231884 0.3768116
25 26 27 28 29 30
1.7768116 1.5768116 0.3768116 -1.8231884 -3.2231884 -2.0231884
31 32 33 34 35 36
8.1768116 7.8000000 11.9768116 11.5768116 8.1768116 8.3768116
37 38 39 40 41 42
7.5768116 7.1768116 6.1768116 3.9768116 2.7768116 3.1768116
43 44 45 46 47 48
13.5768116 15.1768116 14.9768116 11.7768116 8.3768116 8.7768116
49 50 51 52 53 54
7.9768116 7.3768116 5.3768116 4.1768116 3.9768116 3.9768116
55 56 57 58 59 60
13.3768116 14.5768116 13.3768116 3.2000000 2.5768116 0.7768116
61 62 63 64 65 66
1.5768116 -0.8231884 -4.2231884 -5.4231884 -8.4231884 -10.8231884
67 68 69 70 71 72
0.3768116 2.3768116 -2.2231884 -5.2231884 -8.6231884 -7.8231884
> postscript(file="/var/www/html/rcomp/tmp/6jypd1228489851.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -15.4231884 NA
1 -15.6231884 -15.4231884
2 -16.6231884 -15.6231884
3 -18.4231884 -16.6231884
4 -19.6231884 -18.4231884
5 -19.4231884 -19.6231884
6 -11.0000000 -19.4231884
7 -5.6231884 -11.0000000
8 -6.0231884 -5.6231884
9 -6.8231884 -6.0231884
10 -8.8231884 -6.8231884
11 -8.2231884 -8.8231884
12 -6.8231884 -8.2231884
13 -7.2231884 -6.8231884
14 -8.6231884 -7.2231884
15 -8.8231884 -8.6231884
16 -10.4231884 -8.8231884
17 -9.2231884 -10.4231884
18 3.1768116 -9.2231884
19 5.3768116 3.1768116
20 4.9768116 5.3768116
21 2.3768116 4.9768116
22 -1.2231884 2.3768116
23 0.3768116 -1.2231884
24 1.7768116 0.3768116
25 1.5768116 1.7768116
26 0.3768116 1.5768116
27 -1.8231884 0.3768116
28 -3.2231884 -1.8231884
29 -2.0231884 -3.2231884
30 8.1768116 -2.0231884
31 7.8000000 8.1768116
32 11.9768116 7.8000000
33 11.5768116 11.9768116
34 8.1768116 11.5768116
35 8.3768116 8.1768116
36 7.5768116 8.3768116
37 7.1768116 7.5768116
38 6.1768116 7.1768116
39 3.9768116 6.1768116
40 2.7768116 3.9768116
41 3.1768116 2.7768116
42 13.5768116 3.1768116
43 15.1768116 13.5768116
44 14.9768116 15.1768116
45 11.7768116 14.9768116
46 8.3768116 11.7768116
47 8.7768116 8.3768116
48 7.9768116 8.7768116
49 7.3768116 7.9768116
50 5.3768116 7.3768116
51 4.1768116 5.3768116
52 3.9768116 4.1768116
53 3.9768116 3.9768116
54 13.3768116 3.9768116
55 14.5768116 13.3768116
56 13.3768116 14.5768116
57 3.2000000 13.3768116
58 2.5768116 3.2000000
59 0.7768116 2.5768116
60 1.5768116 0.7768116
61 -0.8231884 1.5768116
62 -4.2231884 -0.8231884
63 -5.4231884 -4.2231884
64 -8.4231884 -5.4231884
65 -10.8231884 -8.4231884
66 0.3768116 -10.8231884
67 2.3768116 0.3768116
68 -2.2231884 2.3768116
69 -5.2231884 -2.2231884
70 -8.6231884 -5.2231884
71 -7.8231884 -8.6231884
72 NA -7.8231884
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -15.6231884 -15.4231884
[2,] -16.6231884 -15.6231884
[3,] -18.4231884 -16.6231884
[4,] -19.6231884 -18.4231884
[5,] -19.4231884 -19.6231884
[6,] -11.0000000 -19.4231884
[7,] -5.6231884 -11.0000000
[8,] -6.0231884 -5.6231884
[9,] -6.8231884 -6.0231884
[10,] -8.8231884 -6.8231884
[11,] -8.2231884 -8.8231884
[12,] -6.8231884 -8.2231884
[13,] -7.2231884 -6.8231884
[14,] -8.6231884 -7.2231884
[15,] -8.8231884 -8.6231884
[16,] -10.4231884 -8.8231884
[17,] -9.2231884 -10.4231884
[18,] 3.1768116 -9.2231884
[19,] 5.3768116 3.1768116
[20,] 4.9768116 5.3768116
[21,] 2.3768116 4.9768116
[22,] -1.2231884 2.3768116
[23,] 0.3768116 -1.2231884
[24,] 1.7768116 0.3768116
[25,] 1.5768116 1.7768116
[26,] 0.3768116 1.5768116
[27,] -1.8231884 0.3768116
[28,] -3.2231884 -1.8231884
[29,] -2.0231884 -3.2231884
[30,] 8.1768116 -2.0231884
[31,] 7.8000000 8.1768116
[32,] 11.9768116 7.8000000
[33,] 11.5768116 11.9768116
[34,] 8.1768116 11.5768116
[35,] 8.3768116 8.1768116
[36,] 7.5768116 8.3768116
[37,] 7.1768116 7.5768116
[38,] 6.1768116 7.1768116
[39,] 3.9768116 6.1768116
[40,] 2.7768116 3.9768116
[41,] 3.1768116 2.7768116
[42,] 13.5768116 3.1768116
[43,] 15.1768116 13.5768116
[44,] 14.9768116 15.1768116
[45,] 11.7768116 14.9768116
[46,] 8.3768116 11.7768116
[47,] 8.7768116 8.3768116
[48,] 7.9768116 8.7768116
[49,] 7.3768116 7.9768116
[50,] 5.3768116 7.3768116
[51,] 4.1768116 5.3768116
[52,] 3.9768116 4.1768116
[53,] 3.9768116 3.9768116
[54,] 13.3768116 3.9768116
[55,] 14.5768116 13.3768116
[56,] 13.3768116 14.5768116
[57,] 3.2000000 13.3768116
[58,] 2.5768116 3.2000000
[59,] 0.7768116 2.5768116
[60,] 1.5768116 0.7768116
[61,] -0.8231884 1.5768116
[62,] -4.2231884 -0.8231884
[63,] -5.4231884 -4.2231884
[64,] -8.4231884 -5.4231884
[65,] -10.8231884 -8.4231884
[66,] 0.3768116 -10.8231884
[67,] 2.3768116 0.3768116
[68,] -2.2231884 2.3768116
[69,] -5.2231884 -2.2231884
[70,] -8.6231884 -5.2231884
[71,] -7.8231884 -8.6231884
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -15.6231884 -15.4231884
2 -16.6231884 -15.6231884
3 -18.4231884 -16.6231884
4 -19.6231884 -18.4231884
5 -19.4231884 -19.6231884
6 -11.0000000 -19.4231884
7 -5.6231884 -11.0000000
8 -6.0231884 -5.6231884
9 -6.8231884 -6.0231884
10 -8.8231884 -6.8231884
11 -8.2231884 -8.8231884
12 -6.8231884 -8.2231884
13 -7.2231884 -6.8231884
14 -8.6231884 -7.2231884
15 -8.8231884 -8.6231884
16 -10.4231884 -8.8231884
17 -9.2231884 -10.4231884
18 3.1768116 -9.2231884
19 5.3768116 3.1768116
20 4.9768116 5.3768116
21 2.3768116 4.9768116
22 -1.2231884 2.3768116
23 0.3768116 -1.2231884
24 1.7768116 0.3768116
25 1.5768116 1.7768116
26 0.3768116 1.5768116
27 -1.8231884 0.3768116
28 -3.2231884 -1.8231884
29 -2.0231884 -3.2231884
30 8.1768116 -2.0231884
31 7.8000000 8.1768116
32 11.9768116 7.8000000
33 11.5768116 11.9768116
34 8.1768116 11.5768116
35 8.3768116 8.1768116
36 7.5768116 8.3768116
37 7.1768116 7.5768116
38 6.1768116 7.1768116
39 3.9768116 6.1768116
40 2.7768116 3.9768116
41 3.1768116 2.7768116
42 13.5768116 3.1768116
43 15.1768116 13.5768116
44 14.9768116 15.1768116
45 11.7768116 14.9768116
46 8.3768116 11.7768116
47 8.7768116 8.3768116
48 7.9768116 8.7768116
49 7.3768116 7.9768116
50 5.3768116 7.3768116
51 4.1768116 5.3768116
52 3.9768116 4.1768116
53 3.9768116 3.9768116
54 13.3768116 3.9768116
55 14.5768116 13.3768116
56 13.3768116 14.5768116
57 3.2000000 13.3768116
58 2.5768116 3.2000000
59 0.7768116 2.5768116
60 1.5768116 0.7768116
61 -0.8231884 1.5768116
62 -4.2231884 -0.8231884
63 -5.4231884 -4.2231884
64 -8.4231884 -5.4231884
65 -10.8231884 -8.4231884
66 0.3768116 -10.8231884
67 2.3768116 0.3768116
68 -2.2231884 2.3768116
69 -5.2231884 -2.2231884
70 -8.6231884 -5.2231884
71 -7.8231884 -8.6231884
> 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/7kvcd1228489851.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/88dle1228489851.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/9ws0b1228489851.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/10q3811228489851.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/11k06v1228489851.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/12g0td1228489852.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/136chx1228489852.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/1484tw1228489852.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/15cs761228489852.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/16aj1s1228489852.tab")
+ }
>
> system("convert tmp/19r9z1228489851.ps tmp/19r9z1228489851.png")
> system("convert tmp/2rtez1228489851.ps tmp/2rtez1228489851.png")
> system("convert tmp/33lcu1228489851.ps tmp/33lcu1228489851.png")
> system("convert tmp/458ux1228489851.ps tmp/458ux1228489851.png")
> system("convert tmp/5pecg1228489851.ps tmp/5pecg1228489851.png")
> system("convert tmp/6jypd1228489851.ps tmp/6jypd1228489851.png")
> system("convert tmp/7kvcd1228489851.ps tmp/7kvcd1228489851.png")
> system("convert tmp/88dle1228489851.ps tmp/88dle1228489851.png")
> system("convert tmp/9ws0b1228489851.ps tmp/9ws0b1228489851.png")
> system("convert tmp/10q3811228489851.ps tmp/10q3811228489851.png")
>
>
> proc.time()
user system elapsed
2.637 1.599 3.037