R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(1901,10436,1395,9314,1639,9717,1643,8997,1751,9062,1797,8885,1373,9058,1558,9095,1555,9149,2061,9857,2010,9848,2119,10269,1985,10341,1963,9690,2017,10125,1975,9349,1589,9224,1679,9224,1392,9454,1511,9347,1449,9430,1767,9933,1899,10148,2179,10677,2217,10735,2049,9760,2343,10567,2175,9333,1607,9409,1702,9502,1764,9348,1766,9319,1615,9594,1953,10160,2091,10182,2411,10810,2550,11105,2351,9874,2786,10958,2525,9311,2474,9610,2332,9398,1978,9784,1789,9425,1904,9557,1997,10166,2207,10337,2453,10770,1948,11265,1384,10183,1989,10941,2140,9628,2100,9709,2045,9637,2083,9579,2022,9741,1950,9754,1422,10508,1859,10749,2147,11079),dim=c(2,60),dimnames=list(c('aanbod','invoer'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('aanbod','invoer'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
aanbod invoer M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 1901 10436 1 0 0 0 0 0 0 0 0 0 0
2 1395 9314 0 1 0 0 0 0 0 0 0 0 0
3 1639 9717 0 0 1 0 0 0 0 0 0 0 0
4 1643 8997 0 0 0 1 0 0 0 0 0 0 0
5 1751 9062 0 0 0 0 1 0 0 0 0 0 0
6 1797 8885 0 0 0 0 0 1 0 0 0 0 0
7 1373 9058 0 0 0 0 0 0 1 0 0 0 0
8 1558 9095 0 0 0 0 0 0 0 1 0 0 0
9 1555 9149 0 0 0 0 0 0 0 0 1 0 0
10 2061 9857 0 0 0 0 0 0 0 0 0 1 0
11 2010 9848 0 0 0 0 0 0 0 0 0 0 1
12 2119 10269 0 0 0 0 0 0 0 0 0 0 0
13 1985 10341 1 0 0 0 0 0 0 0 0 0 0
14 1963 9690 0 1 0 0 0 0 0 0 0 0 0
15 2017 10125 0 0 1 0 0 0 0 0 0 0 0
16 1975 9349 0 0 0 1 0 0 0 0 0 0 0
17 1589 9224 0 0 0 0 1 0 0 0 0 0 0
18 1679 9224 0 0 0 0 0 1 0 0 0 0 0
19 1392 9454 0 0 0 0 0 0 1 0 0 0 0
20 1511 9347 0 0 0 0 0 0 0 1 0 0 0
21 1449 9430 0 0 0 0 0 0 0 0 1 0 0
22 1767 9933 0 0 0 0 0 0 0 0 0 1 0
23 1899 10148 0 0 0 0 0 0 0 0 0 0 1
24 2179 10677 0 0 0 0 0 0 0 0 0 0 0
25 2217 10735 1 0 0 0 0 0 0 0 0 0 0
26 2049 9760 0 1 0 0 0 0 0 0 0 0 0
27 2343 10567 0 0 1 0 0 0 0 0 0 0 0
28 2175 9333 0 0 0 1 0 0 0 0 0 0 0
29 1607 9409 0 0 0 0 1 0 0 0 0 0 0
30 1702 9502 0 0 0 0 0 1 0 0 0 0 0
31 1764 9348 0 0 0 0 0 0 1 0 0 0 0
32 1766 9319 0 0 0 0 0 0 0 1 0 0 0
33 1615 9594 0 0 0 0 0 0 0 0 1 0 0
34 1953 10160 0 0 0 0 0 0 0 0 0 1 0
35 2091 10182 0 0 0 0 0 0 0 0 0 0 1
36 2411 10810 0 0 0 0 0 0 0 0 0 0 0
37 2550 11105 1 0 0 0 0 0 0 0 0 0 0
38 2351 9874 0 1 0 0 0 0 0 0 0 0 0
39 2786 10958 0 0 1 0 0 0 0 0 0 0 0
40 2525 9311 0 0 0 1 0 0 0 0 0 0 0
41 2474 9610 0 0 0 0 1 0 0 0 0 0 0
42 2332 9398 0 0 0 0 0 1 0 0 0 0 0
43 1978 9784 0 0 0 0 0 0 1 0 0 0 0
44 1789 9425 0 0 0 0 0 0 0 1 0 0 0
45 1904 9557 0 0 0 0 0 0 0 0 1 0 0
46 1997 10166 0 0 0 0 0 0 0 0 0 1 0
47 2207 10337 0 0 0 0 0 0 0 0 0 0 1
48 2453 10770 0 0 0 0 0 0 0 0 0 0 0
49 1948 11265 1 0 0 0 0 0 0 0 0 0 0
50 1384 10183 0 1 0 0 0 0 0 0 0 0 0
51 1989 10941 0 0 1 0 0 0 0 0 0 0 0
52 2140 9628 0 0 0 1 0 0 0 0 0 0 0
53 2100 9709 0 0 0 0 1 0 0 0 0 0 0
54 2045 9637 0 0 0 0 0 1 0 0 0 0 0
55 2083 9579 0 0 0 0 0 0 1 0 0 0 0
56 2022 9741 0 0 0 0 0 0 0 1 0 0 0
57 1950 9754 0 0 0 0 0 0 0 0 1 0 0
58 1422 10508 0 0 0 0 0 0 0 0 0 1 0
59 1859 10749 0 0 0 0 0 0 0 0 0 0 1
60 2147 11079 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) invoer M1 M2 M3 M4
-1859.3172 0.3844 -162.8956 -65.6092 -7.2875 366.9560
M5 M6 M7 M8 M9 M10
149.1117 184.2033 -53.1560 -19.1998 -96.6215 -192.6227
M11
-68.6255
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-605.39 -198.15 31.56 153.33 490.15
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1859.3172 1360.0245 -1.367 0.17809
invoer 0.3844 0.1263 3.043 0.00383 **
M1 -162.8956 174.7996 -0.932 0.35615
M2 -65.6092 212.4066 -0.309 0.75877
M3 -7.2875 177.7071 -0.041 0.96746
M4 366.9560 248.3365 1.478 0.14617
M5 149.1117 241.3264 0.618 0.53963
M6 184.2033 247.8340 0.743 0.46103
M7 -53.1560 237.7131 -0.224 0.82403
M8 -19.1998 242.8485 -0.079 0.93732
M9 -96.6215 233.2901 -0.414 0.68063
M10 -192.6227 190.2075 -1.013 0.31639
M11 -68.6255 184.4030 -0.372 0.71145
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 276.2 on 47 degrees of freedom
Multiple R-squared: 0.4252, Adjusted R-squared: 0.2784
F-statistic: 2.897 on 12 and 47 DF, p-value: 0.004514
> 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.05351119 0.10702238 0.9464888
[2,] 0.08871982 0.17743963 0.9112802
[3,] 0.13751788 0.27503575 0.8624821
[4,] 0.11711772 0.23423543 0.8828823
[5,] 0.08285753 0.16571506 0.9171425
[6,] 0.07051182 0.14102364 0.9294882
[7,] 0.06778175 0.13556351 0.9322182
[8,] 0.04989776 0.09979553 0.9501022
[9,] 0.02874991 0.05749983 0.9712501
[10,] 0.01862799 0.03725599 0.9813720
[11,] 0.01635493 0.03270986 0.9836451
[12,] 0.01164472 0.02328943 0.9883553
[13,] 0.01080697 0.02161395 0.9891930
[14,] 0.02474444 0.04948888 0.9752556
[15,] 0.03340322 0.06680645 0.9665968
[16,] 0.06215037 0.12430075 0.9378496
[17,] 0.06182713 0.12365425 0.9381729
[18,] 0.05451469 0.10902937 0.9454853
[19,] 0.03356666 0.06713333 0.9664333
[20,] 0.02721715 0.05443430 0.9727828
[21,] 0.01551979 0.03103957 0.9844802
[22,] 0.01753008 0.03506016 0.9824699
[23,] 0.07971711 0.15943422 0.9202829
[24,] 0.49601071 0.99202142 0.5039893
[25,] 0.54800223 0.90399554 0.4519978
[26,] 0.65681357 0.68637286 0.3431864
[27,] 0.60092996 0.79814009 0.3990700
[28,] 0.44929417 0.89858835 0.5507058
[29,] 0.73888107 0.52223787 0.2611189
> postscript(file="/var/www/html/rcomp/tmp/125ke1258560550.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/2tsxg1258560550.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/3w4pb1258560550.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/4xq4v1258560550.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/502hy1258560550.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 = 60
Frequency = 1
1 2 3 4 5 6
-88.3513590 -260.3446000 -229.5782078 -323.0560337 -22.1976004 56.7490197
7 8 9 10 11 12
-196.3922309 -59.5711947 -5.9068234 323.9414396 152.4037893 30.9473141
13 14 15 16 17 18
32.1663287 163.1222361 -8.4120665 -126.3636765 -246.4698678 -191.5614658
19 20 21 22 23 24
-329.6133291 -203.4391662 -219.9222996 0.7272895 -73.9152245 -65.8865446
25 26 27 28 29 30
112.7140239 222.2144662 147.6845865 79.7866709 -299.5832596 -275.4237519
31 32 33 34 35 36
83.1327224 62.3239417 -116.9633605 99.4692357 105.0152873 114.9886926
37 38 39 40 41 42
303.4872402 480.3932410 440.3854719 438.2433986 490.1530012 394.5535062
43 44 45 46 47 48
129.5357557 44.5778902 186.2593179 141.1628554 161.4337968 172.3645611
49 50 51 52 53 54
-360.0162338 -605.3853432 -350.0797840 -68.6103593 78.0977266 15.6826919
55 56 57 58 59 60
313.3370818 156.1085290 156.5331655 -565.3008203 -344.9376488 -252.4140231
> postscript(file="/var/www/html/rcomp/tmp/60ky01258560550.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -88.3513590 NA
1 -260.3446000 -88.3513590
2 -229.5782078 -260.3446000
3 -323.0560337 -229.5782078
4 -22.1976004 -323.0560337
5 56.7490197 -22.1976004
6 -196.3922309 56.7490197
7 -59.5711947 -196.3922309
8 -5.9068234 -59.5711947
9 323.9414396 -5.9068234
10 152.4037893 323.9414396
11 30.9473141 152.4037893
12 32.1663287 30.9473141
13 163.1222361 32.1663287
14 -8.4120665 163.1222361
15 -126.3636765 -8.4120665
16 -246.4698678 -126.3636765
17 -191.5614658 -246.4698678
18 -329.6133291 -191.5614658
19 -203.4391662 -329.6133291
20 -219.9222996 -203.4391662
21 0.7272895 -219.9222996
22 -73.9152245 0.7272895
23 -65.8865446 -73.9152245
24 112.7140239 -65.8865446
25 222.2144662 112.7140239
26 147.6845865 222.2144662
27 79.7866709 147.6845865
28 -299.5832596 79.7866709
29 -275.4237519 -299.5832596
30 83.1327224 -275.4237519
31 62.3239417 83.1327224
32 -116.9633605 62.3239417
33 99.4692357 -116.9633605
34 105.0152873 99.4692357
35 114.9886926 105.0152873
36 303.4872402 114.9886926
37 480.3932410 303.4872402
38 440.3854719 480.3932410
39 438.2433986 440.3854719
40 490.1530012 438.2433986
41 394.5535062 490.1530012
42 129.5357557 394.5535062
43 44.5778902 129.5357557
44 186.2593179 44.5778902
45 141.1628554 186.2593179
46 161.4337968 141.1628554
47 172.3645611 161.4337968
48 -360.0162338 172.3645611
49 -605.3853432 -360.0162338
50 -350.0797840 -605.3853432
51 -68.6103593 -350.0797840
52 78.0977266 -68.6103593
53 15.6826919 78.0977266
54 313.3370818 15.6826919
55 156.1085290 313.3370818
56 156.5331655 156.1085290
57 -565.3008203 156.5331655
58 -344.9376488 -565.3008203
59 -252.4140231 -344.9376488
60 NA -252.4140231
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -260.3446000 -88.3513590
[2,] -229.5782078 -260.3446000
[3,] -323.0560337 -229.5782078
[4,] -22.1976004 -323.0560337
[5,] 56.7490197 -22.1976004
[6,] -196.3922309 56.7490197
[7,] -59.5711947 -196.3922309
[8,] -5.9068234 -59.5711947
[9,] 323.9414396 -5.9068234
[10,] 152.4037893 323.9414396
[11,] 30.9473141 152.4037893
[12,] 32.1663287 30.9473141
[13,] 163.1222361 32.1663287
[14,] -8.4120665 163.1222361
[15,] -126.3636765 -8.4120665
[16,] -246.4698678 -126.3636765
[17,] -191.5614658 -246.4698678
[18,] -329.6133291 -191.5614658
[19,] -203.4391662 -329.6133291
[20,] -219.9222996 -203.4391662
[21,] 0.7272895 -219.9222996
[22,] -73.9152245 0.7272895
[23,] -65.8865446 -73.9152245
[24,] 112.7140239 -65.8865446
[25,] 222.2144662 112.7140239
[26,] 147.6845865 222.2144662
[27,] 79.7866709 147.6845865
[28,] -299.5832596 79.7866709
[29,] -275.4237519 -299.5832596
[30,] 83.1327224 -275.4237519
[31,] 62.3239417 83.1327224
[32,] -116.9633605 62.3239417
[33,] 99.4692357 -116.9633605
[34,] 105.0152873 99.4692357
[35,] 114.9886926 105.0152873
[36,] 303.4872402 114.9886926
[37,] 480.3932410 303.4872402
[38,] 440.3854719 480.3932410
[39,] 438.2433986 440.3854719
[40,] 490.1530012 438.2433986
[41,] 394.5535062 490.1530012
[42,] 129.5357557 394.5535062
[43,] 44.5778902 129.5357557
[44,] 186.2593179 44.5778902
[45,] 141.1628554 186.2593179
[46,] 161.4337968 141.1628554
[47,] 172.3645611 161.4337968
[48,] -360.0162338 172.3645611
[49,] -605.3853432 -360.0162338
[50,] -350.0797840 -605.3853432
[51,] -68.6103593 -350.0797840
[52,] 78.0977266 -68.6103593
[53,] 15.6826919 78.0977266
[54,] 313.3370818 15.6826919
[55,] 156.1085290 313.3370818
[56,] 156.5331655 156.1085290
[57,] -565.3008203 156.5331655
[58,] -344.9376488 -565.3008203
[59,] -252.4140231 -344.9376488
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -260.3446000 -88.3513590
2 -229.5782078 -260.3446000
3 -323.0560337 -229.5782078
4 -22.1976004 -323.0560337
5 56.7490197 -22.1976004
6 -196.3922309 56.7490197
7 -59.5711947 -196.3922309
8 -5.9068234 -59.5711947
9 323.9414396 -5.9068234
10 152.4037893 323.9414396
11 30.9473141 152.4037893
12 32.1663287 30.9473141
13 163.1222361 32.1663287
14 -8.4120665 163.1222361
15 -126.3636765 -8.4120665
16 -246.4698678 -126.3636765
17 -191.5614658 -246.4698678
18 -329.6133291 -191.5614658
19 -203.4391662 -329.6133291
20 -219.9222996 -203.4391662
21 0.7272895 -219.9222996
22 -73.9152245 0.7272895
23 -65.8865446 -73.9152245
24 112.7140239 -65.8865446
25 222.2144662 112.7140239
26 147.6845865 222.2144662
27 79.7866709 147.6845865
28 -299.5832596 79.7866709
29 -275.4237519 -299.5832596
30 83.1327224 -275.4237519
31 62.3239417 83.1327224
32 -116.9633605 62.3239417
33 99.4692357 -116.9633605
34 105.0152873 99.4692357
35 114.9886926 105.0152873
36 303.4872402 114.9886926
37 480.3932410 303.4872402
38 440.3854719 480.3932410
39 438.2433986 440.3854719
40 490.1530012 438.2433986
41 394.5535062 490.1530012
42 129.5357557 394.5535062
43 44.5778902 129.5357557
44 186.2593179 44.5778902
45 141.1628554 186.2593179
46 161.4337968 141.1628554
47 172.3645611 161.4337968
48 -360.0162338 172.3645611
49 -605.3853432 -360.0162338
50 -350.0797840 -605.3853432
51 -68.6103593 -350.0797840
52 78.0977266 -68.6103593
53 15.6826919 78.0977266
54 313.3370818 15.6826919
55 156.1085290 313.3370818
56 156.5331655 156.1085290
57 -565.3008203 156.5331655
58 -344.9376488 -565.3008203
59 -252.4140231 -344.9376488
> 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/71hgy1258560550.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/8ty201258560550.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/9ta7w1258560550.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/10rffv1258560550.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/11do491258560550.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/12bb5b1258560550.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/13nk9n1258560550.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/1475r51258560551.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/15fhu61258560551.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/1672cg1258560551.tab")
+ }
>
> system("convert tmp/125ke1258560550.ps tmp/125ke1258560550.png")
> system("convert tmp/2tsxg1258560550.ps tmp/2tsxg1258560550.png")
> system("convert tmp/3w4pb1258560550.ps tmp/3w4pb1258560550.png")
> system("convert tmp/4xq4v1258560550.ps tmp/4xq4v1258560550.png")
> system("convert tmp/502hy1258560550.ps tmp/502hy1258560550.png")
> system("convert tmp/60ky01258560550.ps tmp/60ky01258560550.png")
> system("convert tmp/71hgy1258560550.ps tmp/71hgy1258560550.png")
> system("convert tmp/8ty201258560550.ps tmp/8ty201258560550.png")
> system("convert tmp/9ta7w1258560550.ps tmp/9ta7w1258560550.png")
> system("convert tmp/10rffv1258560550.ps tmp/10rffv1258560550.png")
>
>
> proc.time()
user system elapsed
2.363 1.528 2.833