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,0,13480.21,0,13673.28,0,13239.71,0,13557.69,0,13901.28,0,13200.58,0,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 0 0 0 0 0 0 0 0 0 1 0 0 45
46 13480.21 0 0 0 0 0 0 0 0 0 0 1 0 46
47 13673.28 0 0 0 0 0 0 0 0 0 0 0 1 47
48 13239.71 0 0 0 0 0 0 0 0 0 0 0 0 48
49 13557.69 0 1 0 0 0 0 0 0 0 0 0 0 49
50 13901.28 0 0 1 0 0 0 0 0 0 0 0 0 50
51 13200.58 0 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
8934.33 -1235.85 74.41 336.56 339.78 744.84
M5 M6 M7 M8 M9 M10
590.15 585.85 365.88 439.43 501.53 300.74
M11 t
96.25 73.22
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1125.053 -359.205 -9.275 368.981 1201.499
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8934.328 324.567 27.527 < 2e-16 ***
Y -1235.855 278.154 -4.443 5.37e-05 ***
M1 74.411 367.392 0.203 0.8404
M2 336.557 385.929 0.872 0.3876
M3 339.778 385.633 0.881 0.3828
M4 744.843 385.713 1.931 0.0595 .
M5 590.154 385.054 1.533 0.1321
M6 585.853 384.482 1.524 0.1343
M7 365.880 383.997 0.953 0.3456
M8 439.433 383.600 1.146 0.2578
M9 501.526 383.291 1.308 0.1971
M10 300.738 383.071 0.785 0.4364
M11 96.245 382.938 0.251 0.8027
t 73.217 5.818 12.584 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 605.4 on 47 degrees of freedom
Multiple R-squared: 0.7967, Adjusted R-squared: 0.7404
F-statistic: 14.16 on 13 and 47 DF, p-value: 4.101e-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,] 7.976835e-02 1.595367e-01 0.9202317
[2,] 4.076976e-02 8.153953e-02 0.9592302
[3,] 1.510746e-02 3.021492e-02 0.9848925
[4,] 1.492142e-02 2.984284e-02 0.9850786
[5,] 5.131553e-03 1.026311e-02 0.9948684
[6,] 2.070076e-03 4.140153e-03 0.9979299
[7,] 8.082217e-04 1.616443e-03 0.9991918
[8,] 4.245335e-04 8.490670e-04 0.9995755
[9,] 2.549898e-04 5.099796e-04 0.9997450
[10,] 7.621956e-05 1.524391e-04 0.9999238
[11,] 2.792194e-05 5.584388e-05 0.9999721
[12,] 1.069089e-05 2.138178e-05 0.9999893
[13,] 4.071150e-06 8.142300e-06 0.9999959
[14,] 1.502996e-06 3.005992e-06 0.9999985
[15,] 4.974766e-07 9.949532e-07 0.9999995
[16,] 4.979253e-07 9.958506e-07 0.9999995
[17,] 1.798518e-06 3.597036e-06 0.9999982
[18,] 6.503863e-07 1.300773e-06 0.9999993
[19,] 2.109085e-07 4.218171e-07 0.9999998
[20,] 1.475444e-07 2.950888e-07 0.9999999
[21,] 5.324368e-07 1.064874e-06 0.9999995
[22,] 2.057236e-05 4.114472e-05 0.9999794
[23,] 5.348663e-05 1.069733e-04 0.9999465
[24,] 8.876373e-04 1.775275e-03 0.9991124
[25,] 1.987848e-03 3.975695e-03 0.9980122
[26,] 3.155695e-03 6.311390e-03 0.9968443
[27,] 6.930464e-03 1.386093e-02 0.9930695
[28,] 6.214466e-02 1.242893e-01 0.9378553
> postscript(file="/var/www/html/freestat/rcomp/tmp/1syqr1227780991.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/freestat/rcomp/tmp/2cju11227780991.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/freestat/rcomp/tmp/3pche1227780991.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/freestat/rcomp/tmp/4zx4f1227780991.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/freestat/rcomp/tmp/5w7g11227780991.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
410.533024 265.030355 268.362355 152.589451 649.481451 642.125451
7 8 9 10 11 12
511.001451 458.901451 -1.848549 397.671451 316.127451 219.865451
13 14 15 16 17 18
244.026847 -294.325822 39.386178 -177.266726 -229.664726 -114.310726
19 20 21 22 23 24
-9.274726 -554.914726 -596.234726 -359.204726 -169.188726 -137.270726
25 26 27 28 29 30
-306.629331 -850.222000 -555.720000 -901.442904 -775.300904 -745.506904
31 32 33 34 35 36
-424.290904 -482.030904 -518.140904 -726.480904 -556.284904 -312.796904
37 38 39 40 41 42
-184.185509 -90.018178 55.573822 -230.239082 -13.497082 36.176918
43 44 45 46 47 48
-180.017082 159.482918 677.122918 877.152918 1201.498918 790.956918
49 50 51 52 53 54
961.308313 969.535645 192.397645 1156.359260 368.981260 181.515260
55 56 57 58 59 60
102.581260 418.561260 439.101260 -189.138740 -792.152740 -560.754740
61
-1125.053344
> postscript(file="/var/www/html/freestat/rcomp/tmp/6qawf1227780991.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 410.533024 NA
1 265.030355 410.533024
2 268.362355 265.030355
3 152.589451 268.362355
4 649.481451 152.589451
5 642.125451 649.481451
6 511.001451 642.125451
7 458.901451 511.001451
8 -1.848549 458.901451
9 397.671451 -1.848549
10 316.127451 397.671451
11 219.865451 316.127451
12 244.026847 219.865451
13 -294.325822 244.026847
14 39.386178 -294.325822
15 -177.266726 39.386178
16 -229.664726 -177.266726
17 -114.310726 -229.664726
18 -9.274726 -114.310726
19 -554.914726 -9.274726
20 -596.234726 -554.914726
21 -359.204726 -596.234726
22 -169.188726 -359.204726
23 -137.270726 -169.188726
24 -306.629331 -137.270726
25 -850.222000 -306.629331
26 -555.720000 -850.222000
27 -901.442904 -555.720000
28 -775.300904 -901.442904
29 -745.506904 -775.300904
30 -424.290904 -745.506904
31 -482.030904 -424.290904
32 -518.140904 -482.030904
33 -726.480904 -518.140904
34 -556.284904 -726.480904
35 -312.796904 -556.284904
36 -184.185509 -312.796904
37 -90.018178 -184.185509
38 55.573822 -90.018178
39 -230.239082 55.573822
40 -13.497082 -230.239082
41 36.176918 -13.497082
42 -180.017082 36.176918
43 159.482918 -180.017082
44 677.122918 159.482918
45 877.152918 677.122918
46 1201.498918 877.152918
47 790.956918 1201.498918
48 961.308313 790.956918
49 969.535645 961.308313
50 192.397645 969.535645
51 1156.359260 192.397645
52 368.981260 1156.359260
53 181.515260 368.981260
54 102.581260 181.515260
55 418.561260 102.581260
56 439.101260 418.561260
57 -189.138740 439.101260
58 -792.152740 -189.138740
59 -560.754740 -792.152740
60 -1125.053344 -560.754740
61 NA -1125.053344
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 265.030355 410.533024
[2,] 268.362355 265.030355
[3,] 152.589451 268.362355
[4,] 649.481451 152.589451
[5,] 642.125451 649.481451
[6,] 511.001451 642.125451
[7,] 458.901451 511.001451
[8,] -1.848549 458.901451
[9,] 397.671451 -1.848549
[10,] 316.127451 397.671451
[11,] 219.865451 316.127451
[12,] 244.026847 219.865451
[13,] -294.325822 244.026847
[14,] 39.386178 -294.325822
[15,] -177.266726 39.386178
[16,] -229.664726 -177.266726
[17,] -114.310726 -229.664726
[18,] -9.274726 -114.310726
[19,] -554.914726 -9.274726
[20,] -596.234726 -554.914726
[21,] -359.204726 -596.234726
[22,] -169.188726 -359.204726
[23,] -137.270726 -169.188726
[24,] -306.629331 -137.270726
[25,] -850.222000 -306.629331
[26,] -555.720000 -850.222000
[27,] -901.442904 -555.720000
[28,] -775.300904 -901.442904
[29,] -745.506904 -775.300904
[30,] -424.290904 -745.506904
[31,] -482.030904 -424.290904
[32,] -518.140904 -482.030904
[33,] -726.480904 -518.140904
[34,] -556.284904 -726.480904
[35,] -312.796904 -556.284904
[36,] -184.185509 -312.796904
[37,] -90.018178 -184.185509
[38,] 55.573822 -90.018178
[39,] -230.239082 55.573822
[40,] -13.497082 -230.239082
[41,] 36.176918 -13.497082
[42,] -180.017082 36.176918
[43,] 159.482918 -180.017082
[44,] 677.122918 159.482918
[45,] 877.152918 677.122918
[46,] 1201.498918 877.152918
[47,] 790.956918 1201.498918
[48,] 961.308313 790.956918
[49,] 969.535645 961.308313
[50,] 192.397645 969.535645
[51,] 1156.359260 192.397645
[52,] 368.981260 1156.359260
[53,] 181.515260 368.981260
[54,] 102.581260 181.515260
[55,] 418.561260 102.581260
[56,] 439.101260 418.561260
[57,] -189.138740 439.101260
[58,] -792.152740 -189.138740
[59,] -560.754740 -792.152740
[60,] -1125.053344 -560.754740
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 265.030355 410.533024
2 268.362355 265.030355
3 152.589451 268.362355
4 649.481451 152.589451
5 642.125451 649.481451
6 511.001451 642.125451
7 458.901451 511.001451
8 -1.848549 458.901451
9 397.671451 -1.848549
10 316.127451 397.671451
11 219.865451 316.127451
12 244.026847 219.865451
13 -294.325822 244.026847
14 39.386178 -294.325822
15 -177.266726 39.386178
16 -229.664726 -177.266726
17 -114.310726 -229.664726
18 -9.274726 -114.310726
19 -554.914726 -9.274726
20 -596.234726 -554.914726
21 -359.204726 -596.234726
22 -169.188726 -359.204726
23 -137.270726 -169.188726
24 -306.629331 -137.270726
25 -850.222000 -306.629331
26 -555.720000 -850.222000
27 -901.442904 -555.720000
28 -775.300904 -901.442904
29 -745.506904 -775.300904
30 -424.290904 -745.506904
31 -482.030904 -424.290904
32 -518.140904 -482.030904
33 -726.480904 -518.140904
34 -556.284904 -726.480904
35 -312.796904 -556.284904
36 -184.185509 -312.796904
37 -90.018178 -184.185509
38 55.573822 -90.018178
39 -230.239082 55.573822
40 -13.497082 -230.239082
41 36.176918 -13.497082
42 -180.017082 36.176918
43 159.482918 -180.017082
44 677.122918 159.482918
45 877.152918 677.122918
46 1201.498918 877.152918
47 790.956918 1201.498918
48 961.308313 790.956918
49 969.535645 961.308313
50 192.397645 969.535645
51 1156.359260 192.397645
52 368.981260 1156.359260
53 181.515260 368.981260
54 102.581260 181.515260
55 418.561260 102.581260
56 439.101260 418.561260
57 -189.138740 439.101260
58 -792.152740 -189.138740
59 -560.754740 -792.152740
60 -1125.053344 -560.754740
> 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/freestat/rcomp/tmp/7hn3l1227780991.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/freestat/rcomp/tmp/8ri3u1227780991.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/freestat/rcomp/tmp/9vezd1227780991.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/freestat/rcomp/tmp/10r3b21227780991.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11ehsl1227780991.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/freestat/rcomp/tmp/12dnrs1227780991.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/freestat/rcomp/tmp/13n64w1227780991.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/freestat/rcomp/tmp/14wjbu1227780991.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/freestat/rcomp/tmp/15kvfc1227780991.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/freestat/rcomp/tmp/16lpcy1227780991.tab")
+ }
>
> system("convert tmp/1syqr1227780991.ps tmp/1syqr1227780991.png")
> system("convert tmp/2cju11227780991.ps tmp/2cju11227780991.png")
> system("convert tmp/3pche1227780991.ps tmp/3pche1227780991.png")
> system("convert tmp/4zx4f1227780991.ps tmp/4zx4f1227780991.png")
> system("convert tmp/5w7g11227780991.ps tmp/5w7g11227780991.png")
> system("convert tmp/6qawf1227780991.ps tmp/6qawf1227780991.png")
> system("convert tmp/7hn3l1227780991.ps tmp/7hn3l1227780991.png")
> system("convert tmp/8ri3u1227780991.ps tmp/8ri3u1227780991.png")
> system("convert tmp/9vezd1227780991.ps tmp/9vezd1227780991.png")
> system("convert tmp/10r3b21227780991.ps tmp/10r3b21227780991.png")
>
>
> proc.time()
user system elapsed
3.596 2.467 4.011