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(15859.4,0,15258.9,0,15498.6,0,15106.5,0,15023.6,0,12083,0,15761.3,0,16942.6,0,15070.3,0,13659.6,0,14768.9,0,14725.1,0,15998.1,0,15370.6,0,14956.9,0,15469.7,0,15101.8,0,11703.7,0,16283.6,0,16726.5,0,14968.9,0,14861,0,14583.3,0,15305.8,0,17903.9,0,16379.4,0,15420.3,0,17870.5,0,15912.8,0,13866.5,0,17823.2,0,17872,0,17422,0,16704.5,0,15991.2,0,16583.6,0,19123.5,0,17838.7,0,17209.4,0,18586.5,0,16258.1,0,15141.6,1,19202.1,1,17746.5,1,19090.1,1,18040.3,1,17515.5,1,17751.8,1,21072.4,1,17170,1,19439.5,1,19795.4,1,17574.9,1,16165.4,1,19464.6,1,19932.1,1,19961.2,1,17343.4,1,18924.2,1,18574.1,1,21350.6,1,18594.6,1,19823.1,1,20844.4,1,19640.2,1,17735.4,1,19813.6,1,22238.5,1,20682.2,1,17818.6,1,21872.1,1,22117,1,21865.9,1),dim=c(2,73),dimnames=list(c('uitvoer','dummy'),1:73))
> y <- array(NA,dim=c(2,73),dimnames=list(c('uitvoer','dummy'),1:73))
> 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
uitvoer dummy
1 15859.4 0
2 15258.9 0
3 15498.6 0
4 15106.5 0
5 15023.6 0
6 12083.0 0
7 15761.3 0
8 16942.6 0
9 15070.3 0
10 13659.6 0
11 14768.9 0
12 14725.1 0
13 15998.1 0
14 15370.6 0
15 14956.9 0
16 15469.7 0
17 15101.8 0
18 11703.7 0
19 16283.6 0
20 16726.5 0
21 14968.9 0
22 14861.0 0
23 14583.3 0
24 15305.8 0
25 17903.9 0
26 16379.4 0
27 15420.3 0
28 17870.5 0
29 15912.8 0
30 13866.5 0
31 17823.2 0
32 17872.0 0
33 17422.0 0
34 16704.5 0
35 15991.2 0
36 16583.6 0
37 19123.5 0
38 17838.7 0
39 17209.4 0
40 18586.5 0
41 16258.1 0
42 15141.6 1
43 19202.1 1
44 17746.5 1
45 19090.1 1
46 18040.3 1
47 17515.5 1
48 17751.8 1
49 21072.4 1
50 17170.0 1
51 19439.5 1
52 19795.4 1
53 17574.9 1
54 16165.4 1
55 19464.6 1
56 19932.1 1
57 19961.2 1
58 17343.4 1
59 18924.2 1
60 18574.1 1
61 21350.6 1
62 18594.6 1
63 19823.1 1
64 20844.4 1
65 19640.2 1
66 17735.4 1
67 19813.6 1
68 22238.5 1
69 20682.2 1
70 17818.6 1
71 21872.1 1
72 22117.0 1
73 21865.9 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) dummy
15850 3347
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4146.393 -989.093 5.184 876.407 3273.407
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15850.1 258.1 61.418 < 2e-16 ***
dummy 3346.8 389.8 8.586 1.37e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1652 on 71 degrees of freedom
Multiple R-squared: 0.5094, Adjusted R-squared: 0.5025
F-statistic: 73.73 on 1 and 71 DF, p-value: 1.368e-12
> 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.01588171 0.03176343 0.9841183
[2,] 0.50653097 0.98693807 0.4934690
[3,] 0.39901960 0.79803921 0.6009804
[4,] 0.43670539 0.87341078 0.5632946
[5,] 0.32046767 0.64093534 0.6795323
[6,] 0.31721886 0.63443772 0.6827811
[7,] 0.23219078 0.46438155 0.7678092
[8,] 0.16544254 0.33088508 0.8345575
[9,] 0.12989227 0.25978454 0.8701077
[10,] 0.08670895 0.17341790 0.9132911
[11,] 0.05641539 0.11283078 0.9435846
[12,] 0.03591613 0.07183225 0.9640839
[13,] 0.02193054 0.04386108 0.9780695
[14,] 0.19341629 0.38683257 0.8065837
[15,] 0.17578089 0.35156177 0.8242191
[16,] 0.17919145 0.35838290 0.8208086
[17,] 0.14206472 0.28412944 0.8579353
[18,] 0.11426688 0.22853376 0.8857331
[19,] 0.09888683 0.19777366 0.9011132
[20,] 0.07735448 0.15470896 0.9226455
[21,] 0.14708894 0.29417789 0.8529111
[22,] 0.12742346 0.25484691 0.8725765
[23,] 0.10293890 0.20587781 0.8970611
[24,] 0.15109635 0.30219270 0.8489036
[25,] 0.12248585 0.24497171 0.8775141
[26,] 0.17211496 0.34422991 0.8278850
[27,] 0.21286739 0.42573479 0.7871326
[28,] 0.24837048 0.49674096 0.7516295
[29,] 0.24560804 0.49121608 0.7543920
[30,] 0.21396149 0.42792298 0.7860385
[31,] 0.18498097 0.36996195 0.8150190
[32,] 0.16066508 0.32133016 0.8393349
[33,] 0.26223239 0.52446477 0.7377676
[34,] 0.25879952 0.51759905 0.7412005
[35,] 0.22715468 0.45430936 0.7728453
[36,] 0.27492765 0.54985530 0.7250724
[37,] 0.22189586 0.44379172 0.7781041
[38,] 0.35231046 0.70462091 0.6476895
[39,] 0.37531905 0.75063811 0.6246809
[40,] 0.34753050 0.69506100 0.6524695
[41,] 0.30679260 0.61358520 0.6932074
[42,] 0.27109964 0.54219928 0.7289004
[43,] 0.26256570 0.52513140 0.7374343
[44,] 0.24710092 0.49420183 0.7528991
[45,] 0.29595381 0.59190761 0.7040462
[46,] 0.32844208 0.65688416 0.6715579
[47,] 0.27708213 0.55416427 0.7229179
[48,] 0.23301762 0.46603524 0.7669824
[49,] 0.24147857 0.48295714 0.7585214
[50,] 0.46253721 0.92507442 0.5374628
[51,] 0.40097691 0.80195382 0.5990231
[52,] 0.34270284 0.68540569 0.6572972
[53,] 0.28451642 0.56903285 0.7154836
[54,] 0.37781140 0.75562281 0.6221886
[55,] 0.33568248 0.67136496 0.6643175
[56,] 0.32687300 0.65374600 0.6731270
[57,] 0.31609256 0.63218513 0.6839074
[58,] 0.30890006 0.61780012 0.6910999
[59,] 0.23934086 0.47868172 0.7606591
[60,] 0.18301649 0.36603298 0.8169835
[61,] 0.12997961 0.25995923 0.8700204
[62,] 0.25871599 0.51743198 0.7412840
[63,] 0.20153177 0.40306355 0.7984682
[64,] 0.17011926 0.34023853 0.8298807
> postscript(file="/var/www/html/rcomp/tmp/1hglz1230627035.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/2ceeo1230627035.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/31kie1230627035.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/48ph31230627035.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/5ltrc1230627035.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 = 73
Frequency = 1
1 2 3 4 5 6
9.307317 -591.192683 -351.492683 -743.592683 -826.492683 -3767.092683
7 8 9 10 11 12
-88.792683 1092.507317 -779.792683 -2190.492683 -1081.192683 -1124.992683
13 14 15 16 17 18
148.007317 -479.492683 -893.192683 -380.392683 -748.292683 -4146.392683
19 20 21 22 23 24
433.507317 876.407317 -881.192683 -989.092683 -1266.792683 -544.292683
25 26 27 28 29 30
2053.807317 529.307317 -429.792683 2020.407317 62.707317 -1983.592683
31 32 33 34 35 36
1973.107317 2021.907317 1571.907317 854.407317 141.107317 733.507317
37 38 39 40 41 42
3273.407317 1988.607317 1359.307317 2736.407317 408.007317 -4055.315625
43 44 45 46 47 48
5.184375 -1450.415625 -106.815625 -1156.615625 -1681.415625 -1445.115625
49 50 51 52 53 54
1875.484375 -2026.915625 242.584375 598.484375 -1622.015625 -3031.515625
55 56 57 58 59 60
267.684375 735.184375 764.284375 -1853.515625 -272.715625 -622.815625
61 62 63 64 65 66
2153.684375 -602.315625 626.184375 1647.484375 443.284375 -1461.515625
67 68 69 70 71 72
616.684375 3041.584375 1485.284375 -1378.315625 2675.184375 2920.084375
73
2668.984375
> postscript(file="/var/www/html/rcomp/tmp/6k3lt1230627035.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 9.307317 NA
1 -591.192683 9.307317
2 -351.492683 -591.192683
3 -743.592683 -351.492683
4 -826.492683 -743.592683
5 -3767.092683 -826.492683
6 -88.792683 -3767.092683
7 1092.507317 -88.792683
8 -779.792683 1092.507317
9 -2190.492683 -779.792683
10 -1081.192683 -2190.492683
11 -1124.992683 -1081.192683
12 148.007317 -1124.992683
13 -479.492683 148.007317
14 -893.192683 -479.492683
15 -380.392683 -893.192683
16 -748.292683 -380.392683
17 -4146.392683 -748.292683
18 433.507317 -4146.392683
19 876.407317 433.507317
20 -881.192683 876.407317
21 -989.092683 -881.192683
22 -1266.792683 -989.092683
23 -544.292683 -1266.792683
24 2053.807317 -544.292683
25 529.307317 2053.807317
26 -429.792683 529.307317
27 2020.407317 -429.792683
28 62.707317 2020.407317
29 -1983.592683 62.707317
30 1973.107317 -1983.592683
31 2021.907317 1973.107317
32 1571.907317 2021.907317
33 854.407317 1571.907317
34 141.107317 854.407317
35 733.507317 141.107317
36 3273.407317 733.507317
37 1988.607317 3273.407317
38 1359.307317 1988.607317
39 2736.407317 1359.307317
40 408.007317 2736.407317
41 -4055.315625 408.007317
42 5.184375 -4055.315625
43 -1450.415625 5.184375
44 -106.815625 -1450.415625
45 -1156.615625 -106.815625
46 -1681.415625 -1156.615625
47 -1445.115625 -1681.415625
48 1875.484375 -1445.115625
49 -2026.915625 1875.484375
50 242.584375 -2026.915625
51 598.484375 242.584375
52 -1622.015625 598.484375
53 -3031.515625 -1622.015625
54 267.684375 -3031.515625
55 735.184375 267.684375
56 764.284375 735.184375
57 -1853.515625 764.284375
58 -272.715625 -1853.515625
59 -622.815625 -272.715625
60 2153.684375 -622.815625
61 -602.315625 2153.684375
62 626.184375 -602.315625
63 1647.484375 626.184375
64 443.284375 1647.484375
65 -1461.515625 443.284375
66 616.684375 -1461.515625
67 3041.584375 616.684375
68 1485.284375 3041.584375
69 -1378.315625 1485.284375
70 2675.184375 -1378.315625
71 2920.084375 2675.184375
72 2668.984375 2920.084375
73 NA 2668.984375
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -591.192683 9.307317
[2,] -351.492683 -591.192683
[3,] -743.592683 -351.492683
[4,] -826.492683 -743.592683
[5,] -3767.092683 -826.492683
[6,] -88.792683 -3767.092683
[7,] 1092.507317 -88.792683
[8,] -779.792683 1092.507317
[9,] -2190.492683 -779.792683
[10,] -1081.192683 -2190.492683
[11,] -1124.992683 -1081.192683
[12,] 148.007317 -1124.992683
[13,] -479.492683 148.007317
[14,] -893.192683 -479.492683
[15,] -380.392683 -893.192683
[16,] -748.292683 -380.392683
[17,] -4146.392683 -748.292683
[18,] 433.507317 -4146.392683
[19,] 876.407317 433.507317
[20,] -881.192683 876.407317
[21,] -989.092683 -881.192683
[22,] -1266.792683 -989.092683
[23,] -544.292683 -1266.792683
[24,] 2053.807317 -544.292683
[25,] 529.307317 2053.807317
[26,] -429.792683 529.307317
[27,] 2020.407317 -429.792683
[28,] 62.707317 2020.407317
[29,] -1983.592683 62.707317
[30,] 1973.107317 -1983.592683
[31,] 2021.907317 1973.107317
[32,] 1571.907317 2021.907317
[33,] 854.407317 1571.907317
[34,] 141.107317 854.407317
[35,] 733.507317 141.107317
[36,] 3273.407317 733.507317
[37,] 1988.607317 3273.407317
[38,] 1359.307317 1988.607317
[39,] 2736.407317 1359.307317
[40,] 408.007317 2736.407317
[41,] -4055.315625 408.007317
[42,] 5.184375 -4055.315625
[43,] -1450.415625 5.184375
[44,] -106.815625 -1450.415625
[45,] -1156.615625 -106.815625
[46,] -1681.415625 -1156.615625
[47,] -1445.115625 -1681.415625
[48,] 1875.484375 -1445.115625
[49,] -2026.915625 1875.484375
[50,] 242.584375 -2026.915625
[51,] 598.484375 242.584375
[52,] -1622.015625 598.484375
[53,] -3031.515625 -1622.015625
[54,] 267.684375 -3031.515625
[55,] 735.184375 267.684375
[56,] 764.284375 735.184375
[57,] -1853.515625 764.284375
[58,] -272.715625 -1853.515625
[59,] -622.815625 -272.715625
[60,] 2153.684375 -622.815625
[61,] -602.315625 2153.684375
[62,] 626.184375 -602.315625
[63,] 1647.484375 626.184375
[64,] 443.284375 1647.484375
[65,] -1461.515625 443.284375
[66,] 616.684375 -1461.515625
[67,] 3041.584375 616.684375
[68,] 1485.284375 3041.584375
[69,] -1378.315625 1485.284375
[70,] 2675.184375 -1378.315625
[71,] 2920.084375 2675.184375
[72,] 2668.984375 2920.084375
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -591.192683 9.307317
2 -351.492683 -591.192683
3 -743.592683 -351.492683
4 -826.492683 -743.592683
5 -3767.092683 -826.492683
6 -88.792683 -3767.092683
7 1092.507317 -88.792683
8 -779.792683 1092.507317
9 -2190.492683 -779.792683
10 -1081.192683 -2190.492683
11 -1124.992683 -1081.192683
12 148.007317 -1124.992683
13 -479.492683 148.007317
14 -893.192683 -479.492683
15 -380.392683 -893.192683
16 -748.292683 -380.392683
17 -4146.392683 -748.292683
18 433.507317 -4146.392683
19 876.407317 433.507317
20 -881.192683 876.407317
21 -989.092683 -881.192683
22 -1266.792683 -989.092683
23 -544.292683 -1266.792683
24 2053.807317 -544.292683
25 529.307317 2053.807317
26 -429.792683 529.307317
27 2020.407317 -429.792683
28 62.707317 2020.407317
29 -1983.592683 62.707317
30 1973.107317 -1983.592683
31 2021.907317 1973.107317
32 1571.907317 2021.907317
33 854.407317 1571.907317
34 141.107317 854.407317
35 733.507317 141.107317
36 3273.407317 733.507317
37 1988.607317 3273.407317
38 1359.307317 1988.607317
39 2736.407317 1359.307317
40 408.007317 2736.407317
41 -4055.315625 408.007317
42 5.184375 -4055.315625
43 -1450.415625 5.184375
44 -106.815625 -1450.415625
45 -1156.615625 -106.815625
46 -1681.415625 -1156.615625
47 -1445.115625 -1681.415625
48 1875.484375 -1445.115625
49 -2026.915625 1875.484375
50 242.584375 -2026.915625
51 598.484375 242.584375
52 -1622.015625 598.484375
53 -3031.515625 -1622.015625
54 267.684375 -3031.515625
55 735.184375 267.684375
56 764.284375 735.184375
57 -1853.515625 764.284375
58 -272.715625 -1853.515625
59 -622.815625 -272.715625
60 2153.684375 -622.815625
61 -602.315625 2153.684375
62 626.184375 -602.315625
63 1647.484375 626.184375
64 443.284375 1647.484375
65 -1461.515625 443.284375
66 616.684375 -1461.515625
67 3041.584375 616.684375
68 1485.284375 3041.584375
69 -1378.315625 1485.284375
70 2675.184375 -1378.315625
71 2920.084375 2675.184375
72 2668.984375 2920.084375
> 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/7ckeu1230627035.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/8p65m1230627035.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/97kcj1230627035.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/10htpv1230627035.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/11w5091230627035.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/12eg3f1230627035.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/13b3141230627035.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/14wj981230627035.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/15rlj01230627035.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/16v3bz1230627035.tab")
+ }
>
> system("convert tmp/1hglz1230627035.ps tmp/1hglz1230627035.png")
> system("convert tmp/2ceeo1230627035.ps tmp/2ceeo1230627035.png")
> system("convert tmp/31kie1230627035.ps tmp/31kie1230627035.png")
> system("convert tmp/48ph31230627035.ps tmp/48ph31230627035.png")
> system("convert tmp/5ltrc1230627035.ps tmp/5ltrc1230627035.png")
> system("convert tmp/6k3lt1230627035.ps tmp/6k3lt1230627035.png")
> system("convert tmp/7ckeu1230627035.ps tmp/7ckeu1230627035.png")
> system("convert tmp/8p65m1230627035.ps tmp/8p65m1230627035.png")
> system("convert tmp/97kcj1230627035.ps tmp/97kcj1230627035.png")
> system("convert tmp/10htpv1230627035.ps tmp/10htpv1230627035.png")
>
>
> proc.time()
user system elapsed
2.622 1.595 3.660