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(512927,0,502831,0,470984,0,471067,0,476049,0,474605,0,470439,0,461251,0,454724,0,455626,0,516847,0,525192,0,522975,0,518585,0,509239,0,512238,0,519164,0,517009,0,509933,0,509127,0,500857,0,506971,0,569323,0,579714,0,577992,0,565464,0,547344,0,554788,0,562325,0,560854,0,555332,0,543599,0,536662,0,542722,0,593530,0,610763,1,612613,1,611324,1,594167,1,595454,1,590865,1,589379,1,584428,1,573100,1,567456,1,569028,1,620735,1,628884,1,628232,1,612117,1,595404,1,597141,1,593408,1,590072,1,579799,1,574205,1,572775,1,572942,1,619567,1,625809,1,619916,1),dim=c(2,61),dimnames=list(c('y','d'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('y','d'),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 = '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
y d M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 512927 0 1 0 0 0 0 0 0 0 0 0 0
2 502831 0 0 1 0 0 0 0 0 0 0 0 0
3 470984 0 0 0 1 0 0 0 0 0 0 0 0
4 471067 0 0 0 0 1 0 0 0 0 0 0 0
5 476049 0 0 0 0 0 1 0 0 0 0 0 0
6 474605 0 0 0 0 0 0 1 0 0 0 0 0
7 470439 0 0 0 0 0 0 0 1 0 0 0 0
8 461251 0 0 0 0 0 0 0 0 1 0 0 0
9 454724 0 0 0 0 0 0 0 0 0 1 0 0
10 455626 0 0 0 0 0 0 0 0 0 0 1 0
11 516847 0 0 0 0 0 0 0 0 0 0 0 1
12 525192 0 0 0 0 0 0 0 0 0 0 0 0
13 522975 0 1 0 0 0 0 0 0 0 0 0 0
14 518585 0 0 1 0 0 0 0 0 0 0 0 0
15 509239 0 0 0 1 0 0 0 0 0 0 0 0
16 512238 0 0 0 0 1 0 0 0 0 0 0 0
17 519164 0 0 0 0 0 1 0 0 0 0 0 0
18 517009 0 0 0 0 0 0 1 0 0 0 0 0
19 509933 0 0 0 0 0 0 0 1 0 0 0 0
20 509127 0 0 0 0 0 0 0 0 1 0 0 0
21 500857 0 0 0 0 0 0 0 0 0 1 0 0
22 506971 0 0 0 0 0 0 0 0 0 0 1 0
23 569323 0 0 0 0 0 0 0 0 0 0 0 1
24 579714 0 0 0 0 0 0 0 0 0 0 0 0
25 577992 0 1 0 0 0 0 0 0 0 0 0 0
26 565464 0 0 1 0 0 0 0 0 0 0 0 0
27 547344 0 0 0 1 0 0 0 0 0 0 0 0
28 554788 0 0 0 0 1 0 0 0 0 0 0 0
29 562325 0 0 0 0 0 1 0 0 0 0 0 0
30 560854 0 0 0 0 0 0 1 0 0 0 0 0
31 555332 0 0 0 0 0 0 0 1 0 0 0 0
32 543599 0 0 0 0 0 0 0 0 1 0 0 0
33 536662 0 0 0 0 0 0 0 0 0 1 0 0
34 542722 0 0 0 0 0 0 0 0 0 0 1 0
35 593530 0 0 0 0 0 0 0 0 0 0 0 1
36 610763 1 0 0 0 0 0 0 0 0 0 0 0
37 612613 1 1 0 0 0 0 0 0 0 0 0 0
38 611324 1 0 1 0 0 0 0 0 0 0 0 0
39 594167 1 0 0 1 0 0 0 0 0 0 0 0
40 595454 1 0 0 0 1 0 0 0 0 0 0 0
41 590865 1 0 0 0 0 1 0 0 0 0 0 0
42 589379 1 0 0 0 0 0 1 0 0 0 0 0
43 584428 1 0 0 0 0 0 0 1 0 0 0 0
44 573100 1 0 0 0 0 0 0 0 1 0 0 0
45 567456 1 0 0 0 0 0 0 0 0 1 0 0
46 569028 1 0 0 0 0 0 0 0 0 0 1 0
47 620735 1 0 0 0 0 0 0 0 0 0 0 1
48 628884 1 0 0 0 0 0 0 0 0 0 0 0
49 628232 1 1 0 0 0 0 0 0 0 0 0 0
50 612117 1 0 1 0 0 0 0 0 0 0 0 0
51 595404 1 0 0 1 0 0 0 0 0 0 0 0
52 597141 1 0 0 0 1 0 0 0 0 0 0 0
53 593408 1 0 0 0 0 1 0 0 0 0 0 0
54 590072 1 0 0 0 0 0 1 0 0 0 0 0
55 579799 1 0 0 0 0 0 0 1 0 0 0 0
56 574205 1 0 0 0 0 0 0 0 1 0 0 0
57 572775 1 0 0 0 0 0 0 0 0 1 0 0
58 572942 1 0 0 0 0 0 0 0 0 0 1 0
59 619567 1 0 0 0 0 0 0 0 0 0 0 1
60 625809 1 0 0 0 0 0 0 0 0 0 0 0
61 619916 1 1 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) d M1 M2 M3 M4
549464 74347 -7529 -17139 -35775 -33065
M5 M6 M7 M8 M9 M10
-30841 -32819 -39217 -46947 -52708 -49745
M11
4797
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-45331.9 -9041.4 -166.2 7251.9 45084.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 549464 13399 41.009 < 2e-16 ***
d 74347 7376 10.080 1.95e-13 ***
M1 -7529 17140 -0.439 0.66245
M2 -17139 17946 -0.955 0.34435
M3 -35776 17946 -1.994 0.05191 .
M4 -33066 17946 -1.842 0.07158 .
M5 -30841 17946 -1.719 0.09214 .
M6 -32819 17946 -1.829 0.07365 .
M7 -39217 17946 -2.185 0.03378 *
M8 -46947 17946 -2.616 0.01186 *
M9 -52708 17946 -2.937 0.00508 **
M10 -49745 17946 -2.772 0.00791 **
M11 4797 17946 0.267 0.79037
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 28280 on 48 degrees of freedom
Multiple R-squared: 0.7401, Adjusted R-squared: 0.6751
F-statistic: 11.39 on 12 and 48 DF, p-value: 2.561e-10
> 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.9682194 6.356116e-02 3.178058e-02
[2,] 0.9930800 1.383995e-02 6.919974e-03
[3,] 0.9988345 2.330978e-03 1.165489e-03
[4,] 0.9998734 2.532780e-04 1.266390e-04
[5,] 0.9999897 2.068639e-05 1.034320e-05
[6,] 0.9999998 4.075336e-07 2.037668e-07
[7,] 1.0000000 1.527247e-09 7.636234e-10
[8,] 1.0000000 7.734184e-11 3.867092e-11
[9,] 1.0000000 5.228694e-11 2.614347e-11
[10,] 1.0000000 1.852061e-11 9.260304e-12
[11,] 1.0000000 4.802389e-12 2.401195e-12
[12,] 1.0000000 3.805457e-13 1.902729e-13
[13,] 1.0000000 9.281420e-14 4.640710e-14
[14,] 1.0000000 2.403306e-13 1.201653e-13
[15,] 1.0000000 7.602795e-13 3.801397e-13
[16,] 1.0000000 2.553991e-12 1.276995e-12
[17,] 1.0000000 1.229375e-11 6.146873e-12
[18,] 1.0000000 4.503362e-11 2.251681e-11
[19,] 1.0000000 2.319081e-10 1.159541e-10
[20,] 1.0000000 1.463607e-09 7.318033e-10
[21,] 1.0000000 1.503146e-10 7.515732e-11
[22,] 1.0000000 2.838859e-11 1.419429e-11
[23,] 1.0000000 4.168378e-10 2.084189e-10
[24,] 1.0000000 5.640692e-09 2.820346e-09
[25,] 1.0000000 7.005472e-08 3.502736e-08
[26,] 0.9999996 7.563933e-07 3.781967e-07
[27,] 0.9999956 8.863818e-06 4.431909e-06
[28,] 0.9999675 6.496728e-05 3.248364e-05
[29,] 0.9996617 6.765170e-04 3.382585e-04
[30,] 0.9978109 4.378119e-03 2.189059e-03
> postscript(file="/var/www/html/rcomp/tmp/1nbm61227104436.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/27zpn1227104436.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/33rzq1227104436.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/4c27j1227104436.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/5vy8x1227104436.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
-29008.8231 -29494.5252 -42704.9252 -45331.9252 -42574.5252 -42040.1252
7 8 9 10 11 12
-39808.5252 -41266.7252 -42032.1252 -44093.1252 -37414.7252 -24272.3878
13 14 15 16 17 18
-18960.8231 -13740.5252 -4449.9252 -4160.9252 540.4748 363.8748
19 20 21 22 23 24
-314.5252 6609.2748 4100.8748 7251.8748 15061.2748 30249.6122
25 26 27 28 29 30
36056.1769 33138.4748 33655.0748 38389.0748 43701.4748 44208.8748
31 32 33 34 35 36
45084.4748 41081.2748 39905.8748 43002.8748 39268.2748 -13048.0748
37 38 39 40 41 42
-3669.5102 4651.7878 6131.3878 4708.3878 -2105.2122 -1612.8122
43 44 45 46 47 48
-166.2122 -3764.4122 -3646.8122 -5037.8122 -7873.4122 5072.9252
49 50 51 52 53 54
11949.4898 5444.7878 7368.3878 6395.3878 437.7878 -919.8122
55 56 57 58 59 60
-4795.2122 -2659.4122 1672.1878 -1123.8122 -9041.4122 1997.9252
61
3633.4898
> postscript(file="/var/www/html/rcomp/tmp/6ma0s1227104436.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 -29008.8231 NA
1 -29494.5252 -29008.8231
2 -42704.9252 -29494.5252
3 -45331.9252 -42704.9252
4 -42574.5252 -45331.9252
5 -42040.1252 -42574.5252
6 -39808.5252 -42040.1252
7 -41266.7252 -39808.5252
8 -42032.1252 -41266.7252
9 -44093.1252 -42032.1252
10 -37414.7252 -44093.1252
11 -24272.3878 -37414.7252
12 -18960.8231 -24272.3878
13 -13740.5252 -18960.8231
14 -4449.9252 -13740.5252
15 -4160.9252 -4449.9252
16 540.4748 -4160.9252
17 363.8748 540.4748
18 -314.5252 363.8748
19 6609.2748 -314.5252
20 4100.8748 6609.2748
21 7251.8748 4100.8748
22 15061.2748 7251.8748
23 30249.6122 15061.2748
24 36056.1769 30249.6122
25 33138.4748 36056.1769
26 33655.0748 33138.4748
27 38389.0748 33655.0748
28 43701.4748 38389.0748
29 44208.8748 43701.4748
30 45084.4748 44208.8748
31 41081.2748 45084.4748
32 39905.8748 41081.2748
33 43002.8748 39905.8748
34 39268.2748 43002.8748
35 -13048.0748 39268.2748
36 -3669.5102 -13048.0748
37 4651.7878 -3669.5102
38 6131.3878 4651.7878
39 4708.3878 6131.3878
40 -2105.2122 4708.3878
41 -1612.8122 -2105.2122
42 -166.2122 -1612.8122
43 -3764.4122 -166.2122
44 -3646.8122 -3764.4122
45 -5037.8122 -3646.8122
46 -7873.4122 -5037.8122
47 5072.9252 -7873.4122
48 11949.4898 5072.9252
49 5444.7878 11949.4898
50 7368.3878 5444.7878
51 6395.3878 7368.3878
52 437.7878 6395.3878
53 -919.8122 437.7878
54 -4795.2122 -919.8122
55 -2659.4122 -4795.2122
56 1672.1878 -2659.4122
57 -1123.8122 1672.1878
58 -9041.4122 -1123.8122
59 1997.9252 -9041.4122
60 3633.4898 1997.9252
61 NA 3633.4898
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -29494.5252 -29008.8231
[2,] -42704.9252 -29494.5252
[3,] -45331.9252 -42704.9252
[4,] -42574.5252 -45331.9252
[5,] -42040.1252 -42574.5252
[6,] -39808.5252 -42040.1252
[7,] -41266.7252 -39808.5252
[8,] -42032.1252 -41266.7252
[9,] -44093.1252 -42032.1252
[10,] -37414.7252 -44093.1252
[11,] -24272.3878 -37414.7252
[12,] -18960.8231 -24272.3878
[13,] -13740.5252 -18960.8231
[14,] -4449.9252 -13740.5252
[15,] -4160.9252 -4449.9252
[16,] 540.4748 -4160.9252
[17,] 363.8748 540.4748
[18,] -314.5252 363.8748
[19,] 6609.2748 -314.5252
[20,] 4100.8748 6609.2748
[21,] 7251.8748 4100.8748
[22,] 15061.2748 7251.8748
[23,] 30249.6122 15061.2748
[24,] 36056.1769 30249.6122
[25,] 33138.4748 36056.1769
[26,] 33655.0748 33138.4748
[27,] 38389.0748 33655.0748
[28,] 43701.4748 38389.0748
[29,] 44208.8748 43701.4748
[30,] 45084.4748 44208.8748
[31,] 41081.2748 45084.4748
[32,] 39905.8748 41081.2748
[33,] 43002.8748 39905.8748
[34,] 39268.2748 43002.8748
[35,] -13048.0748 39268.2748
[36,] -3669.5102 -13048.0748
[37,] 4651.7878 -3669.5102
[38,] 6131.3878 4651.7878
[39,] 4708.3878 6131.3878
[40,] -2105.2122 4708.3878
[41,] -1612.8122 -2105.2122
[42,] -166.2122 -1612.8122
[43,] -3764.4122 -166.2122
[44,] -3646.8122 -3764.4122
[45,] -5037.8122 -3646.8122
[46,] -7873.4122 -5037.8122
[47,] 5072.9252 -7873.4122
[48,] 11949.4898 5072.9252
[49,] 5444.7878 11949.4898
[50,] 7368.3878 5444.7878
[51,] 6395.3878 7368.3878
[52,] 437.7878 6395.3878
[53,] -919.8122 437.7878
[54,] -4795.2122 -919.8122
[55,] -2659.4122 -4795.2122
[56,] 1672.1878 -2659.4122
[57,] -1123.8122 1672.1878
[58,] -9041.4122 -1123.8122
[59,] 1997.9252 -9041.4122
[60,] 3633.4898 1997.9252
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -29494.5252 -29008.8231
2 -42704.9252 -29494.5252
3 -45331.9252 -42704.9252
4 -42574.5252 -45331.9252
5 -42040.1252 -42574.5252
6 -39808.5252 -42040.1252
7 -41266.7252 -39808.5252
8 -42032.1252 -41266.7252
9 -44093.1252 -42032.1252
10 -37414.7252 -44093.1252
11 -24272.3878 -37414.7252
12 -18960.8231 -24272.3878
13 -13740.5252 -18960.8231
14 -4449.9252 -13740.5252
15 -4160.9252 -4449.9252
16 540.4748 -4160.9252
17 363.8748 540.4748
18 -314.5252 363.8748
19 6609.2748 -314.5252
20 4100.8748 6609.2748
21 7251.8748 4100.8748
22 15061.2748 7251.8748
23 30249.6122 15061.2748
24 36056.1769 30249.6122
25 33138.4748 36056.1769
26 33655.0748 33138.4748
27 38389.0748 33655.0748
28 43701.4748 38389.0748
29 44208.8748 43701.4748
30 45084.4748 44208.8748
31 41081.2748 45084.4748
32 39905.8748 41081.2748
33 43002.8748 39905.8748
34 39268.2748 43002.8748
35 -13048.0748 39268.2748
36 -3669.5102 -13048.0748
37 4651.7878 -3669.5102
38 6131.3878 4651.7878
39 4708.3878 6131.3878
40 -2105.2122 4708.3878
41 -1612.8122 -2105.2122
42 -166.2122 -1612.8122
43 -3764.4122 -166.2122
44 -3646.8122 -3764.4122
45 -5037.8122 -3646.8122
46 -7873.4122 -5037.8122
47 5072.9252 -7873.4122
48 11949.4898 5072.9252
49 5444.7878 11949.4898
50 7368.3878 5444.7878
51 6395.3878 7368.3878
52 437.7878 6395.3878
53 -919.8122 437.7878
54 -4795.2122 -919.8122
55 -2659.4122 -4795.2122
56 1672.1878 -2659.4122
57 -1123.8122 1672.1878
58 -9041.4122 -1123.8122
59 1997.9252 -9041.4122
60 3633.4898 1997.9252
> 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/76yqy1227104436.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/8iu4t1227104436.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/9osu21227104436.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/10hd1e1227104436.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/11dyhq1227104436.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/120q8s1227104436.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/13phel1227104436.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/141kpt1227104436.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/15gho01227104436.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/16712s1227104436.tab")
+ }
>
> system("convert tmp/1nbm61227104436.ps tmp/1nbm61227104436.png")
> system("convert tmp/27zpn1227104436.ps tmp/27zpn1227104436.png")
> system("convert tmp/33rzq1227104436.ps tmp/33rzq1227104436.png")
> system("convert tmp/4c27j1227104436.ps tmp/4c27j1227104436.png")
> system("convert tmp/5vy8x1227104436.ps tmp/5vy8x1227104436.png")
> system("convert tmp/6ma0s1227104436.ps tmp/6ma0s1227104436.png")
> system("convert tmp/76yqy1227104436.ps tmp/76yqy1227104436.png")
> system("convert tmp/8iu4t1227104436.ps tmp/8iu4t1227104436.png")
> system("convert tmp/9osu21227104436.ps tmp/9osu21227104436.png")
> system("convert tmp/10hd1e1227104436.ps tmp/10hd1e1227104436.png")
>
>
> proc.time()
user system elapsed
2.401 1.564 2.983