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(10144,112,10751,304,11752,794,13808,901,16203,1232,17432,1240,18014,1032,16956,1145,17982,1588,19435,2264,19990,2209,20154,2917,10327,243,9807,558,10862,1238,13743,1502,16458,2000,18466,2146,18810,2066,17361,2046,17411,1952,18517,2771,18525,3278,17859,4000,9499,410,9490,1107,9255,1622,10758,1986,12375,2036,14617,2400,15427,2736,14136,2901,14308,2883,15293,3747,15679,4075,16319,4996,11196,575,11169,999,12158,1411,14251,1493,16237,1846,19706,2899,18960,2372,18537,2856,19103,3468,19691,4193,19464,4440,17264,4186,8957,655,9703,1453,9166,1989,9519,2209,10535,2667,11526,3005,9630,2195,7061,2236,6021,2489,4728,2651,2657,2636,1264,2819),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 10144 112 1 0 0 0 0 0 0 0 0 0 0
2 10751 304 0 1 0 0 0 0 0 0 0 0 0
3 11752 794 0 0 1 0 0 0 0 0 0 0 0
4 13808 901 0 0 0 1 0 0 0 0 0 0 0
5 16203 1232 0 0 0 0 1 0 0 0 0 0 0
6 17432 1240 0 0 0 0 0 1 0 0 0 0 0
7 18014 1032 0 0 0 0 0 0 1 0 0 0 0
8 16956 1145 0 0 0 0 0 0 0 1 0 0 0
9 17982 1588 0 0 0 0 0 0 0 0 1 0 0
10 19435 2264 0 0 0 0 0 0 0 0 0 1 0
11 19990 2209 0 0 0 0 0 0 0 0 0 0 1
12 20154 2917 0 0 0 0 0 0 0 0 0 0 0
13 10327 243 1 0 0 0 0 0 0 0 0 0 0
14 9807 558 0 1 0 0 0 0 0 0 0 0 0
15 10862 1238 0 0 1 0 0 0 0 0 0 0 0
16 13743 1502 0 0 0 1 0 0 0 0 0 0 0
17 16458 2000 0 0 0 0 1 0 0 0 0 0 0
18 18466 2146 0 0 0 0 0 1 0 0 0 0 0
19 18810 2066 0 0 0 0 0 0 1 0 0 0 0
20 17361 2046 0 0 0 0 0 0 0 1 0 0 0
21 17411 1952 0 0 0 0 0 0 0 0 1 0 0
22 18517 2771 0 0 0 0 0 0 0 0 0 1 0
23 18525 3278 0 0 0 0 0 0 0 0 0 0 1
24 17859 4000 0 0 0 0 0 0 0 0 0 0 0
25 9499 410 1 0 0 0 0 0 0 0 0 0 0
26 9490 1107 0 1 0 0 0 0 0 0 0 0 0
27 9255 1622 0 0 1 0 0 0 0 0 0 0 0
28 10758 1986 0 0 0 1 0 0 0 0 0 0 0
29 12375 2036 0 0 0 0 1 0 0 0 0 0 0
30 14617 2400 0 0 0 0 0 1 0 0 0 0 0
31 15427 2736 0 0 0 0 0 0 1 0 0 0 0
32 14136 2901 0 0 0 0 0 0 0 1 0 0 0
33 14308 2883 0 0 0 0 0 0 0 0 1 0 0
34 15293 3747 0 0 0 0 0 0 0 0 0 1 0
35 15679 4075 0 0 0 0 0 0 0 0 0 0 1
36 16319 4996 0 0 0 0 0 0 0 0 0 0 0
37 11196 575 1 0 0 0 0 0 0 0 0 0 0
38 11169 999 0 1 0 0 0 0 0 0 0 0 0
39 12158 1411 0 0 1 0 0 0 0 0 0 0 0
40 14251 1493 0 0 0 1 0 0 0 0 0 0 0
41 16237 1846 0 0 0 0 1 0 0 0 0 0 0
42 19706 2899 0 0 0 0 0 1 0 0 0 0 0
43 18960 2372 0 0 0 0 0 0 1 0 0 0 0
44 18537 2856 0 0 0 0 0 0 0 1 0 0 0
45 19103 3468 0 0 0 0 0 0 0 0 1 0 0
46 19691 4193 0 0 0 0 0 0 0 0 0 1 0
47 19464 4440 0 0 0 0 0 0 0 0 0 0 1
48 17264 4186 0 0 0 0 0 0 0 0 0 0 0
49 8957 655 1 0 0 0 0 0 0 0 0 0 0
50 9703 1453 0 1 0 0 0 0 0 0 0 0 0
51 9166 1989 0 0 1 0 0 0 0 0 0 0 0
52 9519 2209 0 0 0 1 0 0 0 0 0 0 0
53 10535 2667 0 0 0 0 1 0 0 0 0 0 0
54 11526 3005 0 0 0 0 0 1 0 0 0 0 0
55 9630 2195 0 0 0 0 0 0 1 0 0 0 0
56 7061 2236 0 0 0 0 0 0 0 1 0 0 0
57 6021 2489 0 0 0 0 0 0 0 0 1 0 0
58 4728 2651 0 0 0 0 0 0 0 0 0 1 0
59 2657 2636 0 0 0 0 0 0 0 0 0 0 1
60 1264 2819 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) X M1 M2 M3 M4
13159.5926 0.3733 -3283.9382 -3305.6620 -3047.6403 -1347.8622
M5 M6 M7 M8 M9 M10
471.7634 2317.0385 2232.0745 815.6162 881.1235 1206.5789
M11
861.2235
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12948 -1030 1200 2627 5905
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13159.5926 4215.5893 3.122 0.00307 **
X 0.3733 0.9781 0.382 0.70443
M1 -3283.9382 4371.6769 -0.751 0.45629
M2 -3305.6620 4024.2760 -0.821 0.41555
M3 -3047.6403 3679.5200 -0.828 0.41170
M4 -1347.8622 3555.0620 -0.379 0.70629
M5 471.7634 3368.5947 0.140 0.88922
M6 2317.0385 3186.2196 0.727 0.47071
M7 2232.0745 3305.8432 0.675 0.50286
M8 815.6162 3231.3608 0.252 0.80183
M9 881.1235 3128.6629 0.282 0.77946
M10 1206.5789 2927.0393 0.412 0.68205
M11 861.2235 2889.9458 0.298 0.76701
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4515 on 47 degrees of freedom
Multiple R-squared: 0.2396, Adjusted R-squared: 0.04544
F-statistic: 1.234 on 12 and 47 DF, p-value: 0.2893
> 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,] 1.033663e-03 2.067326e-03 0.9989663
[2,] 1.783977e-04 3.567954e-04 0.9998216
[3,] 5.798838e-05 1.159768e-04 0.9999420
[4,] 7.719033e-06 1.543807e-05 0.9999923
[5,] 8.931420e-07 1.786284e-06 0.9999991
[6,] 2.501240e-07 5.002480e-07 0.9999997
[7,] 1.738460e-07 3.476920e-07 0.9999998
[8,] 2.984665e-07 5.969330e-07 0.9999997
[9,] 6.950495e-07 1.390099e-06 0.9999993
[10,] 1.531686e-07 3.063372e-07 0.9999998
[11,] 2.825173e-08 5.650345e-08 1.0000000
[12,] 2.789950e-08 5.579900e-08 1.0000000
[13,] 7.632081e-08 1.526416e-07 0.9999999
[14,] 7.683902e-07 1.536780e-06 0.9999992
[15,] 7.816446e-07 1.563289e-06 0.9999992
[16,] 2.713944e-07 5.427887e-07 0.9999997
[17,] 7.543285e-08 1.508657e-07 0.9999999
[18,] 2.885719e-08 5.771438e-08 1.0000000
[19,] 9.190551e-09 1.838110e-08 1.0000000
[20,] 2.258821e-09 4.517643e-09 1.0000000
[21,] 7.183316e-10 1.436663e-09 1.0000000
[22,] 4.076237e-10 8.152474e-10 1.0000000
[23,] 3.008741e-10 6.017482e-10 1.0000000
[24,] 4.967484e-10 9.934967e-10 1.0000000
[25,] 4.027975e-09 8.055950e-09 1.0000000
[26,] 2.454846e-06 4.909692e-06 0.9999975
[27,] 1.773773e-03 3.547546e-03 0.9982262
[28,] 5.415563e-02 1.083113e-01 0.9458444
[29,] 3.057390e-01 6.114780e-01 0.6942610
> postscript(file="/var/www/html/rcomp/tmp/1umwh1258721721.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/21x6p1258721721.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/34c8t1258721721.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/4qtrd1258721721.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/57x7w1258721721.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
226.5363 783.5871 1343.6497 1659.9288 2111.7419 1492.4804
7 8 9 10 11 12
2237.0902 2553.3659 3348.4880 4223.6836 5144.5703 5905.4994
13 14 15 16 17 18
360.6344 -255.2304 287.9058 1370.5771 2080.0496 2188.2731
19 20 21 22 23 24
2647.1008 2622.0251 2641.6078 3116.4219 3280.5155 3206.2185
25 26 27 28 29 30
-529.7063 -777.1706 -1462.4404 -1795.0987 -2016.3891 -1755.5444
31 32 33 34 35 36
-986.0083 -922.1440 -808.9320 -471.9162 136.9976 1294.4144
37 38 39 40 41 42
1105.6997 942.1455 1519.3253 1881.9368 1916.5374 3147.1802
43 44 45 46 47 48
2682.8719 3495.6543 3767.6891 3759.5932 3785.7441 2541.7852
49 50 51 52 53 54
-1163.1641 -693.3315 -1688.4405 -3117.3440 -4091.9397 -5072.3893
55 56 57 58 59 60
-6581.0545 -7748.9014 -8948.8529 -10627.7824 -12347.8276 -12947.9175
> postscript(file="/var/www/html/rcomp/tmp/66y9r1258721721.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 226.5363 NA
1 783.5871 226.5363
2 1343.6497 783.5871
3 1659.9288 1343.6497
4 2111.7419 1659.9288
5 1492.4804 2111.7419
6 2237.0902 1492.4804
7 2553.3659 2237.0902
8 3348.4880 2553.3659
9 4223.6836 3348.4880
10 5144.5703 4223.6836
11 5905.4994 5144.5703
12 360.6344 5905.4994
13 -255.2304 360.6344
14 287.9058 -255.2304
15 1370.5771 287.9058
16 2080.0496 1370.5771
17 2188.2731 2080.0496
18 2647.1008 2188.2731
19 2622.0251 2647.1008
20 2641.6078 2622.0251
21 3116.4219 2641.6078
22 3280.5155 3116.4219
23 3206.2185 3280.5155
24 -529.7063 3206.2185
25 -777.1706 -529.7063
26 -1462.4404 -777.1706
27 -1795.0987 -1462.4404
28 -2016.3891 -1795.0987
29 -1755.5444 -2016.3891
30 -986.0083 -1755.5444
31 -922.1440 -986.0083
32 -808.9320 -922.1440
33 -471.9162 -808.9320
34 136.9976 -471.9162
35 1294.4144 136.9976
36 1105.6997 1294.4144
37 942.1455 1105.6997
38 1519.3253 942.1455
39 1881.9368 1519.3253
40 1916.5374 1881.9368
41 3147.1802 1916.5374
42 2682.8719 3147.1802
43 3495.6543 2682.8719
44 3767.6891 3495.6543
45 3759.5932 3767.6891
46 3785.7441 3759.5932
47 2541.7852 3785.7441
48 -1163.1641 2541.7852
49 -693.3315 -1163.1641
50 -1688.4405 -693.3315
51 -3117.3440 -1688.4405
52 -4091.9397 -3117.3440
53 -5072.3893 -4091.9397
54 -6581.0545 -5072.3893
55 -7748.9014 -6581.0545
56 -8948.8529 -7748.9014
57 -10627.7824 -8948.8529
58 -12347.8276 -10627.7824
59 -12947.9175 -12347.8276
60 NA -12947.9175
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 783.5871 226.5363
[2,] 1343.6497 783.5871
[3,] 1659.9288 1343.6497
[4,] 2111.7419 1659.9288
[5,] 1492.4804 2111.7419
[6,] 2237.0902 1492.4804
[7,] 2553.3659 2237.0902
[8,] 3348.4880 2553.3659
[9,] 4223.6836 3348.4880
[10,] 5144.5703 4223.6836
[11,] 5905.4994 5144.5703
[12,] 360.6344 5905.4994
[13,] -255.2304 360.6344
[14,] 287.9058 -255.2304
[15,] 1370.5771 287.9058
[16,] 2080.0496 1370.5771
[17,] 2188.2731 2080.0496
[18,] 2647.1008 2188.2731
[19,] 2622.0251 2647.1008
[20,] 2641.6078 2622.0251
[21,] 3116.4219 2641.6078
[22,] 3280.5155 3116.4219
[23,] 3206.2185 3280.5155
[24,] -529.7063 3206.2185
[25,] -777.1706 -529.7063
[26,] -1462.4404 -777.1706
[27,] -1795.0987 -1462.4404
[28,] -2016.3891 -1795.0987
[29,] -1755.5444 -2016.3891
[30,] -986.0083 -1755.5444
[31,] -922.1440 -986.0083
[32,] -808.9320 -922.1440
[33,] -471.9162 -808.9320
[34,] 136.9976 -471.9162
[35,] 1294.4144 136.9976
[36,] 1105.6997 1294.4144
[37,] 942.1455 1105.6997
[38,] 1519.3253 942.1455
[39,] 1881.9368 1519.3253
[40,] 1916.5374 1881.9368
[41,] 3147.1802 1916.5374
[42,] 2682.8719 3147.1802
[43,] 3495.6543 2682.8719
[44,] 3767.6891 3495.6543
[45,] 3759.5932 3767.6891
[46,] 3785.7441 3759.5932
[47,] 2541.7852 3785.7441
[48,] -1163.1641 2541.7852
[49,] -693.3315 -1163.1641
[50,] -1688.4405 -693.3315
[51,] -3117.3440 -1688.4405
[52,] -4091.9397 -3117.3440
[53,] -5072.3893 -4091.9397
[54,] -6581.0545 -5072.3893
[55,] -7748.9014 -6581.0545
[56,] -8948.8529 -7748.9014
[57,] -10627.7824 -8948.8529
[58,] -12347.8276 -10627.7824
[59,] -12947.9175 -12347.8276
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 783.5871 226.5363
2 1343.6497 783.5871
3 1659.9288 1343.6497
4 2111.7419 1659.9288
5 1492.4804 2111.7419
6 2237.0902 1492.4804
7 2553.3659 2237.0902
8 3348.4880 2553.3659
9 4223.6836 3348.4880
10 5144.5703 4223.6836
11 5905.4994 5144.5703
12 360.6344 5905.4994
13 -255.2304 360.6344
14 287.9058 -255.2304
15 1370.5771 287.9058
16 2080.0496 1370.5771
17 2188.2731 2080.0496
18 2647.1008 2188.2731
19 2622.0251 2647.1008
20 2641.6078 2622.0251
21 3116.4219 2641.6078
22 3280.5155 3116.4219
23 3206.2185 3280.5155
24 -529.7063 3206.2185
25 -777.1706 -529.7063
26 -1462.4404 -777.1706
27 -1795.0987 -1462.4404
28 -2016.3891 -1795.0987
29 -1755.5444 -2016.3891
30 -986.0083 -1755.5444
31 -922.1440 -986.0083
32 -808.9320 -922.1440
33 -471.9162 -808.9320
34 136.9976 -471.9162
35 1294.4144 136.9976
36 1105.6997 1294.4144
37 942.1455 1105.6997
38 1519.3253 942.1455
39 1881.9368 1519.3253
40 1916.5374 1881.9368
41 3147.1802 1916.5374
42 2682.8719 3147.1802
43 3495.6543 2682.8719
44 3767.6891 3495.6543
45 3759.5932 3767.6891
46 3785.7441 3759.5932
47 2541.7852 3785.7441
48 -1163.1641 2541.7852
49 -693.3315 -1163.1641
50 -1688.4405 -693.3315
51 -3117.3440 -1688.4405
52 -4091.9397 -3117.3440
53 -5072.3893 -4091.9397
54 -6581.0545 -5072.3893
55 -7748.9014 -6581.0545
56 -8948.8529 -7748.9014
57 -10627.7824 -8948.8529
58 -12347.8276 -10627.7824
59 -12947.9175 -12347.8276
> 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/7yxe71258721721.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/812bl1258721721.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/9mb391258721721.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/1021nm1258721721.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/11yfah1258721721.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/12d1vo1258721721.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/132vdy1258721721.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/14kd5g1258721721.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/15iux81258721721.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/168agb1258721721.tab")
+ }
>
> system("convert tmp/1umwh1258721721.ps tmp/1umwh1258721721.png")
> system("convert tmp/21x6p1258721721.ps tmp/21x6p1258721721.png")
> system("convert tmp/34c8t1258721721.ps tmp/34c8t1258721721.png")
> system("convert tmp/4qtrd1258721721.ps tmp/4qtrd1258721721.png")
> system("convert tmp/57x7w1258721721.ps tmp/57x7w1258721721.png")
> system("convert tmp/66y9r1258721721.ps tmp/66y9r1258721721.png")
> system("convert tmp/7yxe71258721721.ps tmp/7yxe71258721721.png")
> system("convert tmp/812bl1258721721.ps tmp/812bl1258721721.png")
> system("convert tmp/9mb391258721721.ps tmp/9mb391258721721.png")
> system("convert tmp/1021nm1258721721.ps tmp/1021nm1258721721.png")
>
>
> proc.time()
user system elapsed
2.440 1.584 6.571