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(9492.49,0,9682.35,0,9762.12,0,10124.63,0,10540.05,0,10601.61,0,10323.73,0,10418.4,0,10092.96,0,10364.91,0,10152.09,0,10032.8,0,10204.59,0,10001.6,0,10411.75,0,10673.38,0,10539.51,0,10723.78,0,10682.06,0,10283.19,0,10377.18,0,10486.64,0,10545.38,0,10554.27,0,10532.54,0,10324.31,0,10695.25,0,10827.81,0,10872.48,0,10971.19,0,11145.65,0,11234.68,0,11333.88,0,10997.97,0,11036.89,0,11257.35,0,11533.59,0,11963.12,0,12185.15,0,12377.62,0,12512.89,0,12631.48,0,12268.53,0,12754.8,0,13407.75,1,13480.21,1,13673.28,1,13239.71,1,13557.69,1,13901.28,1,13200.58,1,13406.97,1,12538.12,1,12419.57,1,12193.88,1,12656.63,1,12812.48,1,12056.67,1,11322.38,1,11530.75,1,11114.08,1),dim=c(2,61),dimnames=list(c('X','Y'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('X','Y'),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
X Y M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 9492.49 0 1 0 0 0 0 0 0 0 0 0 0 1
2 9682.35 0 0 1 0 0 0 0 0 0 0 0 0 2
3 9762.12 0 0 0 1 0 0 0 0 0 0 0 0 3
4 10124.63 0 0 0 0 1 0 0 0 0 0 0 0 4
5 10540.05 0 0 0 0 0 1 0 0 0 0 0 0 5
6 10601.61 0 0 0 0 0 0 1 0 0 0 0 0 6
7 10323.73 0 0 0 0 0 0 0 1 0 0 0 0 7
8 10418.40 0 0 0 0 0 0 0 0 1 0 0 0 8
9 10092.96 0 0 0 0 0 0 0 0 0 1 0 0 9
10 10364.91 0 0 0 0 0 0 0 0 0 0 1 0 10
11 10152.09 0 0 0 0 0 0 0 0 0 0 0 1 11
12 10032.80 0 0 0 0 0 0 0 0 0 0 0 0 12
13 10204.59 0 1 0 0 0 0 0 0 0 0 0 0 13
14 10001.60 0 0 1 0 0 0 0 0 0 0 0 0 14
15 10411.75 0 0 0 1 0 0 0 0 0 0 0 0 15
16 10673.38 0 0 0 0 1 0 0 0 0 0 0 0 16
17 10539.51 0 0 0 0 0 1 0 0 0 0 0 0 17
18 10723.78 0 0 0 0 0 0 1 0 0 0 0 0 18
19 10682.06 0 0 0 0 0 0 0 1 0 0 0 0 19
20 10283.19 0 0 0 0 0 0 0 0 1 0 0 0 20
21 10377.18 0 0 0 0 0 0 0 0 0 1 0 0 21
22 10486.64 0 0 0 0 0 0 0 0 0 0 1 0 22
23 10545.38 0 0 0 0 0 0 0 0 0 0 0 1 23
24 10554.27 0 0 0 0 0 0 0 0 0 0 0 0 24
25 10532.54 0 1 0 0 0 0 0 0 0 0 0 0 25
26 10324.31 0 0 1 0 0 0 0 0 0 0 0 0 26
27 10695.25 0 0 0 1 0 0 0 0 0 0 0 0 27
28 10827.81 0 0 0 0 1 0 0 0 0 0 0 0 28
29 10872.48 0 0 0 0 0 1 0 0 0 0 0 0 29
30 10971.19 0 0 0 0 0 0 1 0 0 0 0 0 30
31 11145.65 0 0 0 0 0 0 0 1 0 0 0 0 31
32 11234.68 0 0 0 0 0 0 0 0 1 0 0 0 32
33 11333.88 0 0 0 0 0 0 0 0 0 1 0 0 33
34 10997.97 0 0 0 0 0 0 0 0 0 0 1 0 34
35 11036.89 0 0 0 0 0 0 0 0 0 0 0 1 35
36 11257.35 0 0 0 0 0 0 0 0 0 0 0 0 36
37 11533.59 0 1 0 0 0 0 0 0 0 0 0 0 37
38 11963.12 0 0 1 0 0 0 0 0 0 0 0 0 38
39 12185.15 0 0 0 1 0 0 0 0 0 0 0 0 39
40 12377.62 0 0 0 0 1 0 0 0 0 0 0 0 40
41 12512.89 0 0 0 0 0 1 0 0 0 0 0 0 41
42 12631.48 0 0 0 0 0 0 1 0 0 0 0 0 42
43 12268.53 0 0 0 0 0 0 0 1 0 0 0 0 43
44 12754.80 0 0 0 0 0 0 0 0 1 0 0 0 44
45 13407.75 1 0 0 0 0 0 0 0 0 1 0 0 45
46 13480.21 1 0 0 0 0 0 0 0 0 0 1 0 46
47 13673.28 1 0 0 0 0 0 0 0 0 0 0 1 47
48 13239.71 1 0 0 0 0 0 0 0 0 0 0 0 48
49 13557.69 1 1 0 0 0 0 0 0 0 0 0 0 49
50 13901.28 1 0 1 0 0 0 0 0 0 0 0 0 50
51 13200.58 1 0 0 1 0 0 0 0 0 0 0 0 51
52 13406.97 1 0 0 0 1 0 0 0 0 0 0 0 52
53 12538.12 1 0 0 0 0 1 0 0 0 0 0 0 53
54 12419.57 1 0 0 0 0 0 1 0 0 0 0 0 54
55 12193.88 1 0 0 0 0 0 0 1 0 0 0 0 55
56 12656.63 1 0 0 0 0 0 0 0 1 0 0 0 56
57 12812.48 1 0 0 0 0 0 0 0 0 1 0 0 57
58 12056.67 1 0 0 0 0 0 0 0 0 0 1 0 58
59 11322.38 1 0 0 0 0 0 0 0 0 0 0 1 59
60 11530.75 1 0 0 0 0 0 0 0 0 0 0 0 60
61 11114.08 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) Y M1 M2 M3 M4
9472.90 619.91 13.36 420.57 452.51 639.11
M5 M6 M7 M8 M9 M10
513.14 537.55 346.29 448.56 415.38 243.31
M11 t
67.53 44.50
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1706.79 -374.35 -53.14 360.88 1421.29
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9472.90 374.36 25.305 < 2e-16 ***
Y 619.91 321.19 1.930 0.0597 .
M1 13.36 421.40 0.032 0.9748
M2 420.57 442.18 0.951 0.3464
M3 452.50 441.68 1.025 0.3108
M4 639.11 441.33 1.448 0.1542
M5 513.14 441.12 1.163 0.2506
M6 537.55 441.07 1.219 0.2290
M7 346.29 441.17 0.785 0.4364
M8 448.56 441.41 1.016 0.3147
M9 415.38 439.88 0.944 0.3498
M10 243.31 439.50 0.554 0.5825
M11 67.53 439.28 0.154 0.8785
t 44.50 8.12 5.481 1.63e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 694.4 on 47 degrees of freedom
Multiple R-squared: 0.7324, Adjusted R-squared: 0.6584
F-statistic: 9.897 on 13 and 47 DF, p-value: 1.674e-09
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 2.940152e-02 5.880304e-02 0.9705985
[2,] 1.010467e-02 2.020934e-02 0.9898953
[3,] 2.250879e-03 4.501758e-03 0.9977491
[4,] 1.713304e-03 3.426607e-03 0.9982867
[5,] 4.355288e-04 8.710576e-04 0.9995645
[6,] 1.177188e-04 2.354376e-04 0.9998823
[7,] 2.591447e-05 5.182894e-05 0.9999741
[8,] 6.735986e-06 1.347197e-05 0.9999933
[9,] 1.734663e-06 3.469327e-06 0.9999983
[10,] 6.653892e-07 1.330778e-06 0.9999993
[11,] 2.151738e-07 4.303476e-07 0.9999998
[12,] 9.524211e-08 1.904842e-07 0.9999999
[13,] 4.892747e-08 9.785494e-08 1.0000000
[14,] 3.283050e-08 6.566099e-08 1.0000000
[15,] 1.696934e-08 3.393868e-08 1.0000000
[16,] 1.177935e-07 2.355870e-07 0.9999999
[17,] 9.870022e-07 1.974004e-06 0.9999990
[18,] 1.331901e-06 2.663803e-06 0.9999987
[19,] 2.182574e-06 4.365148e-06 0.9999978
[20,] 8.615768e-06 1.723154e-05 0.9999914
[21,] 1.294448e-04 2.588895e-04 0.9998706
[22,] 2.449101e-02 4.898202e-02 0.9755090
[23,] 8.004454e-02 1.600891e-01 0.9199555
[24,] 2.221596e-01 4.443192e-01 0.7778404
[25,] 1.954082e-01 3.908164e-01 0.8045918
[26,] 1.522313e-01 3.044627e-01 0.8477687
[27,] 8.869005e-02 1.773801e-01 0.9113100
[28,] 6.214466e-02 1.242893e-01 0.9378553
> postscript(file="/var/www/html/rcomp/tmp/1aqbc1227779405.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/2gu9q1227779405.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/3jhsj1227779405.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/44kkg1227779405.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/5nx3u1227779405.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
-38.27291 -300.12242 -296.79042 -165.39242 331.49958 324.14358
7 8 9 10 11 12
193.01958 140.91958 -195.84814 203.67186 122.12786 25.86586
13 14 15 16 17 18
139.78844 -514.91107 -181.19907 -150.68107 -203.07907 -87.72507
19 20 21 22 23 24
17.31093 -528.32907 -445.66679 -208.63679 -18.62079 13.29721
25 26 27 28 29 30
-66.30021 -726.23972 -431.73772 -530.28972 -404.14772 -374.35372
31 32 33 34 35 36
-53.13772 -110.87772 -23.00544 -231.34544 -61.14944 182.33856
37 38 39 40 41 42
400.71115 378.53163 524.12363 485.48163 702.22363 751.89763
43 44 45 46 47 48
535.70363 875.20363 896.91451 1096.94451 1421.29051 1010.74851
49 50 51 52 53 54
1270.86109 1162.74157 385.60357 360.88157 -426.49643 -613.96243
55 56 57 58 59 60
-692.89643 -376.91643 -232.39414 -860.63414 -1463.64814 -1232.25014
61
-1706.78756
> postscript(file="/var/www/html/rcomp/tmp/616j51227779405.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 -38.27291 NA
1 -300.12242 -38.27291
2 -296.79042 -300.12242
3 -165.39242 -296.79042
4 331.49958 -165.39242
5 324.14358 331.49958
6 193.01958 324.14358
7 140.91958 193.01958
8 -195.84814 140.91958
9 203.67186 -195.84814
10 122.12786 203.67186
11 25.86586 122.12786
12 139.78844 25.86586
13 -514.91107 139.78844
14 -181.19907 -514.91107
15 -150.68107 -181.19907
16 -203.07907 -150.68107
17 -87.72507 -203.07907
18 17.31093 -87.72507
19 -528.32907 17.31093
20 -445.66679 -528.32907
21 -208.63679 -445.66679
22 -18.62079 -208.63679
23 13.29721 -18.62079
24 -66.30021 13.29721
25 -726.23972 -66.30021
26 -431.73772 -726.23972
27 -530.28972 -431.73772
28 -404.14772 -530.28972
29 -374.35372 -404.14772
30 -53.13772 -374.35372
31 -110.87772 -53.13772
32 -23.00544 -110.87772
33 -231.34544 -23.00544
34 -61.14944 -231.34544
35 182.33856 -61.14944
36 400.71115 182.33856
37 378.53163 400.71115
38 524.12363 378.53163
39 485.48163 524.12363
40 702.22363 485.48163
41 751.89763 702.22363
42 535.70363 751.89763
43 875.20363 535.70363
44 896.91451 875.20363
45 1096.94451 896.91451
46 1421.29051 1096.94451
47 1010.74851 1421.29051
48 1270.86109 1010.74851
49 1162.74157 1270.86109
50 385.60357 1162.74157
51 360.88157 385.60357
52 -426.49643 360.88157
53 -613.96243 -426.49643
54 -692.89643 -613.96243
55 -376.91643 -692.89643
56 -232.39414 -376.91643
57 -860.63414 -232.39414
58 -1463.64814 -860.63414
59 -1232.25014 -1463.64814
60 -1706.78756 -1232.25014
61 NA -1706.78756
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -300.12242 -38.27291
[2,] -296.79042 -300.12242
[3,] -165.39242 -296.79042
[4,] 331.49958 -165.39242
[5,] 324.14358 331.49958
[6,] 193.01958 324.14358
[7,] 140.91958 193.01958
[8,] -195.84814 140.91958
[9,] 203.67186 -195.84814
[10,] 122.12786 203.67186
[11,] 25.86586 122.12786
[12,] 139.78844 25.86586
[13,] -514.91107 139.78844
[14,] -181.19907 -514.91107
[15,] -150.68107 -181.19907
[16,] -203.07907 -150.68107
[17,] -87.72507 -203.07907
[18,] 17.31093 -87.72507
[19,] -528.32907 17.31093
[20,] -445.66679 -528.32907
[21,] -208.63679 -445.66679
[22,] -18.62079 -208.63679
[23,] 13.29721 -18.62079
[24,] -66.30021 13.29721
[25,] -726.23972 -66.30021
[26,] -431.73772 -726.23972
[27,] -530.28972 -431.73772
[28,] -404.14772 -530.28972
[29,] -374.35372 -404.14772
[30,] -53.13772 -374.35372
[31,] -110.87772 -53.13772
[32,] -23.00544 -110.87772
[33,] -231.34544 -23.00544
[34,] -61.14944 -231.34544
[35,] 182.33856 -61.14944
[36,] 400.71115 182.33856
[37,] 378.53163 400.71115
[38,] 524.12363 378.53163
[39,] 485.48163 524.12363
[40,] 702.22363 485.48163
[41,] 751.89763 702.22363
[42,] 535.70363 751.89763
[43,] 875.20363 535.70363
[44,] 896.91451 875.20363
[45,] 1096.94451 896.91451
[46,] 1421.29051 1096.94451
[47,] 1010.74851 1421.29051
[48,] 1270.86109 1010.74851
[49,] 1162.74157 1270.86109
[50,] 385.60357 1162.74157
[51,] 360.88157 385.60357
[52,] -426.49643 360.88157
[53,] -613.96243 -426.49643
[54,] -692.89643 -613.96243
[55,] -376.91643 -692.89643
[56,] -232.39414 -376.91643
[57,] -860.63414 -232.39414
[58,] -1463.64814 -860.63414
[59,] -1232.25014 -1463.64814
[60,] -1706.78756 -1232.25014
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -300.12242 -38.27291
2 -296.79042 -300.12242
3 -165.39242 -296.79042
4 331.49958 -165.39242
5 324.14358 331.49958
6 193.01958 324.14358
7 140.91958 193.01958
8 -195.84814 140.91958
9 203.67186 -195.84814
10 122.12786 203.67186
11 25.86586 122.12786
12 139.78844 25.86586
13 -514.91107 139.78844
14 -181.19907 -514.91107
15 -150.68107 -181.19907
16 -203.07907 -150.68107
17 -87.72507 -203.07907
18 17.31093 -87.72507
19 -528.32907 17.31093
20 -445.66679 -528.32907
21 -208.63679 -445.66679
22 -18.62079 -208.63679
23 13.29721 -18.62079
24 -66.30021 13.29721
25 -726.23972 -66.30021
26 -431.73772 -726.23972
27 -530.28972 -431.73772
28 -404.14772 -530.28972
29 -374.35372 -404.14772
30 -53.13772 -374.35372
31 -110.87772 -53.13772
32 -23.00544 -110.87772
33 -231.34544 -23.00544
34 -61.14944 -231.34544
35 182.33856 -61.14944
36 400.71115 182.33856
37 378.53163 400.71115
38 524.12363 378.53163
39 485.48163 524.12363
40 702.22363 485.48163
41 751.89763 702.22363
42 535.70363 751.89763
43 875.20363 535.70363
44 896.91451 875.20363
45 1096.94451 896.91451
46 1421.29051 1096.94451
47 1010.74851 1421.29051
48 1270.86109 1010.74851
49 1162.74157 1270.86109
50 385.60357 1162.74157
51 360.88157 385.60357
52 -426.49643 360.88157
53 -613.96243 -426.49643
54 -692.89643 -613.96243
55 -376.91643 -692.89643
56 -232.39414 -376.91643
57 -860.63414 -232.39414
58 -1463.64814 -860.63414
59 -1232.25014 -1463.64814
60 -1706.78756 -1232.25014
> 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/78lkp1227779405.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/8t4o61227779405.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/98snu1227779405.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/1021eu1227779405.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/114tc41227779405.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/12eit91227779405.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/13bfyr1227779406.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/14g6lt1227779406.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/155ho61227779406.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/16lvm81227779406.tab")
+ }
>
> system("convert tmp/1aqbc1227779405.ps tmp/1aqbc1227779405.png")
> system("convert tmp/2gu9q1227779405.ps tmp/2gu9q1227779405.png")
> system("convert tmp/3jhsj1227779405.ps tmp/3jhsj1227779405.png")
> system("convert tmp/44kkg1227779405.ps tmp/44kkg1227779405.png")
> system("convert tmp/5nx3u1227779405.ps tmp/5nx3u1227779405.png")
> system("convert tmp/616j51227779405.ps tmp/616j51227779405.png")
> system("convert tmp/78lkp1227779405.ps tmp/78lkp1227779405.png")
> system("convert tmp/8t4o61227779405.ps tmp/8t4o61227779405.png")
> system("convert tmp/98snu1227779405.ps tmp/98snu1227779405.png")
> system("convert tmp/1021eu1227779405.ps tmp/1021eu1227779405.png")
>
>
> proc.time()
user system elapsed
2.645 1.754 3.255