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(562,0,561,0,555,0,544,0,537,0,543,0,594,0,611,0,613,0,611,0,594,0,595,0,591,0,589,0,584,0,573,0,567,0,569,0,621,0,629,0,628,0,612,0,595,0,597,0,593,0,590,0,580,0,574,0,573,0,573,0,620,0,626,0,620,0,588,0,566,0,557,0,561,0,549,0,532,0,526,0,511,0,499,0,555,0,565,0,542,0,527,1,510,1,514,1,517,1,508,1,493,1,490,1,469,1,478,1,528,1,534,1,518,1,506,1,502,1,516,1,528,1),dim=c(2,61),dimnames=list(c('Y','X'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Y','X'),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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 562 0 1 0 0 0 0 0 0 0 0 0 0 1
2 561 0 0 1 0 0 0 0 0 0 0 0 0 2
3 555 0 0 0 1 0 0 0 0 0 0 0 0 3
4 544 0 0 0 0 1 0 0 0 0 0 0 0 4
5 537 0 0 0 0 0 1 0 0 0 0 0 0 5
6 543 0 0 0 0 0 0 1 0 0 0 0 0 6
7 594 0 0 0 0 0 0 0 1 0 0 0 0 7
8 611 0 0 0 0 0 0 0 0 1 0 0 0 8
9 613 0 0 0 0 0 0 0 0 0 1 0 0 9
10 611 0 0 0 0 0 0 0 0 0 0 1 0 10
11 594 0 0 0 0 0 0 0 0 0 0 0 1 11
12 595 0 0 0 0 0 0 0 0 0 0 0 0 12
13 591 0 1 0 0 0 0 0 0 0 0 0 0 13
14 589 0 0 1 0 0 0 0 0 0 0 0 0 14
15 584 0 0 0 1 0 0 0 0 0 0 0 0 15
16 573 0 0 0 0 1 0 0 0 0 0 0 0 16
17 567 0 0 0 0 0 1 0 0 0 0 0 0 17
18 569 0 0 0 0 0 0 1 0 0 0 0 0 18
19 621 0 0 0 0 0 0 0 1 0 0 0 0 19
20 629 0 0 0 0 0 0 0 0 1 0 0 0 20
21 628 0 0 0 0 0 0 0 0 0 1 0 0 21
22 612 0 0 0 0 0 0 0 0 0 0 1 0 22
23 595 0 0 0 0 0 0 0 0 0 0 0 1 23
24 597 0 0 0 0 0 0 0 0 0 0 0 0 24
25 593 0 1 0 0 0 0 0 0 0 0 0 0 25
26 590 0 0 1 0 0 0 0 0 0 0 0 0 26
27 580 0 0 0 1 0 0 0 0 0 0 0 0 27
28 574 0 0 0 0 1 0 0 0 0 0 0 0 28
29 573 0 0 0 0 0 1 0 0 0 0 0 0 29
30 573 0 0 0 0 0 0 1 0 0 0 0 0 30
31 620 0 0 0 0 0 0 0 1 0 0 0 0 31
32 626 0 0 0 0 0 0 0 0 1 0 0 0 32
33 620 0 0 0 0 0 0 0 0 0 1 0 0 33
34 588 0 0 0 0 0 0 0 0 0 0 1 0 34
35 566 0 0 0 0 0 0 0 0 0 0 0 1 35
36 557 0 0 0 0 0 0 0 0 0 0 0 0 36
37 561 0 1 0 0 0 0 0 0 0 0 0 0 37
38 549 0 0 1 0 0 0 0 0 0 0 0 0 38
39 532 0 0 0 1 0 0 0 0 0 0 0 0 39
40 526 0 0 0 0 1 0 0 0 0 0 0 0 40
41 511 0 0 0 0 0 1 0 0 0 0 0 0 41
42 499 0 0 0 0 0 0 1 0 0 0 0 0 42
43 555 0 0 0 0 0 0 0 1 0 0 0 0 43
44 565 0 0 0 0 0 0 0 0 1 0 0 0 44
45 542 0 0 0 0 0 0 0 0 0 1 0 0 45
46 527 1 0 0 0 0 0 0 0 0 0 1 0 46
47 510 1 0 0 0 0 0 0 0 0 0 0 1 47
48 514 1 0 0 0 0 0 0 0 0 0 0 0 48
49 517 1 1 0 0 0 0 0 0 0 0 0 0 49
50 508 1 0 1 0 0 0 0 0 0 0 0 0 50
51 493 1 0 0 1 0 0 0 0 0 0 0 0 51
52 490 1 0 0 0 1 0 0 0 0 0 0 0 52
53 469 1 0 0 0 0 1 0 0 0 0 0 0 53
54 478 1 0 0 0 0 0 1 0 0 0 0 0 54
55 528 1 0 0 0 0 0 0 1 0 0 0 0 55
56 534 1 0 0 0 0 0 0 0 1 0 0 0 56
57 518 1 0 0 0 0 0 0 0 0 1 0 0 57
58 506 1 0 0 0 0 0 0 0 0 0 1 0 58
59 502 1 0 0 0 0 0 0 0 0 0 0 1 59
60 516 1 0 0 0 0 0 0 0 0 0 0 0 60
61 528 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) X M1 M2 M3 M4
603.158 -47.108 -4.234 -13.742 -23.550 -30.158
M5 M6 M7 M8 M9 M10
-39.366 -37.574 14.418 24.610 16.602 11.416
M11 t
-3.192 -0.792
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-42.117 -15.317 1.096 14.852 32.178
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 603.1576 11.1134 54.273 < 2e-16 ***
X -47.1084 9.5242 -4.946 1.01e-05 ***
M1 -4.2342 12.5785 -0.337 0.73790
M2 -13.7423 13.1999 -1.041 0.30316
M3 -23.5502 13.1856 -1.786 0.08054 .
M4 -30.1582 13.1756 -2.289 0.02662 *
M5 -39.3661 13.1698 -2.989 0.00444 **
M6 -37.5741 13.1683 -2.853 0.00641 **
M7 14.4180 13.1710 1.095 0.27924
M8 24.6101 13.1780 1.868 0.06807 .
M9 16.6021 13.1892 1.259 0.21433
M10 11.4159 13.1191 0.870 0.38863
M11 -3.1921 13.1127 -0.243 0.80873
t -0.7921 0.2366 -3.347 0.00161 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.73 on 47 degrees of freedom
Multiple R-squared: 0.8075, Adjusted R-squared: 0.7543
F-statistic: 15.17 on 13 and 47 DF, p-value: 1.213e-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.919417e-05 1.583883e-04 0.999920806
[2,] 5.682262e-05 1.136452e-04 0.999943177
[3,] 5.789089e-06 1.157818e-05 0.999994211
[4,] 2.677368e-04 5.354736e-04 0.999732263
[5,] 6.636292e-04 1.327258e-03 0.999336371
[6,] 1.129028e-02 2.258056e-02 0.988709719
[7,] 2.209591e-02 4.419183e-02 0.977904086
[8,] 2.489657e-02 4.979315e-02 0.975103425
[9,] 2.225341e-02 4.450681e-02 0.977746594
[10,] 1.528577e-02 3.057154e-02 0.984714231
[11,] 1.249916e-02 2.499832e-02 0.987500840
[12,] 6.961432e-03 1.392286e-02 0.993038568
[13,] 4.765994e-03 9.531987e-03 0.995234006
[14,] 4.472059e-03 8.944119e-03 0.995527941
[15,] 5.080942e-03 1.016188e-02 0.994919058
[16,] 1.100544e-02 2.201088e-02 0.988994562
[17,] 3.022986e-01 6.045971e-01 0.697701429
[18,] 8.916030e-01 2.167939e-01 0.108396974
[19,] 9.911811e-01 1.763780e-02 0.008818901
[20,] 9.963944e-01 7.211246e-03 0.003605623
[21,] 9.949634e-01 1.007330e-02 0.005036648
[22,] 9.948224e-01 1.035520e-02 0.005177600
[23,] 9.942908e-01 1.141837e-02 0.005709187
[24,] 9.904309e-01 1.913830e-02 0.009569149
[25,] 9.916681e-01 1.666371e-02 0.008331857
[26,] 9.851549e-01 2.969026e-02 0.014845132
[27,] 9.628119e-01 7.437618e-02 0.037188092
[28,] 9.102714e-01 1.794571e-01 0.089728560
> postscript(file="/var/www/html/rcomp/tmp/13hhj1258570537.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/2m7kl1258570537.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/3vzkc1258570537.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/4k31b1258570537.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/54wwr1258570537.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
-36.1313466 -26.8311810 -22.2311810 -25.8311810 -22.8311810 -17.8311810
7 8 9 10 11 12
-18.0311810 -10.4311810 0.3688190 4.3471302 2.7471302 1.3471302
13 14 15 16 17 18
2.3733996 10.6735651 16.2735651 12.6735651 16.6735651 17.6735651
19 20 21 22 23 24
18.4735651 17.0735651 24.8735651 14.8518764 13.2518764 12.8518764
25 26 27 28 29 30
13.8781457 21.1783113 21.7783113 23.1783113 32.1783113 31.1783113
31 32 33 34 35 36
26.9783113 23.5783113 26.3783113 0.3566225 -6.2433775 -17.6433775
37 38 39 40 41 42
-8.6171082 -10.3169426 -16.7169426 -15.3169426 -20.3169426 -33.3169426
43 44 45 46 47 48
-28.5169426 -27.9169426 -42.1169426 -4.0301876 -5.6301876 -4.0301876
49 50 51 52 53 54
3.9960817 5.2962472 0.8962472 5.2962472 -5.7037528 2.2962472
55 56 57 58 59 60
1.0962472 -2.3037528 -9.5037528 -15.5254415 -4.1254415 7.4745585
61
24.5008278
> postscript(file="/var/www/html/rcomp/tmp/666as1258570537.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 -36.1313466 NA
1 -26.8311810 -36.1313466
2 -22.2311810 -26.8311810
3 -25.8311810 -22.2311810
4 -22.8311810 -25.8311810
5 -17.8311810 -22.8311810
6 -18.0311810 -17.8311810
7 -10.4311810 -18.0311810
8 0.3688190 -10.4311810
9 4.3471302 0.3688190
10 2.7471302 4.3471302
11 1.3471302 2.7471302
12 2.3733996 1.3471302
13 10.6735651 2.3733996
14 16.2735651 10.6735651
15 12.6735651 16.2735651
16 16.6735651 12.6735651
17 17.6735651 16.6735651
18 18.4735651 17.6735651
19 17.0735651 18.4735651
20 24.8735651 17.0735651
21 14.8518764 24.8735651
22 13.2518764 14.8518764
23 12.8518764 13.2518764
24 13.8781457 12.8518764
25 21.1783113 13.8781457
26 21.7783113 21.1783113
27 23.1783113 21.7783113
28 32.1783113 23.1783113
29 31.1783113 32.1783113
30 26.9783113 31.1783113
31 23.5783113 26.9783113
32 26.3783113 23.5783113
33 0.3566225 26.3783113
34 -6.2433775 0.3566225
35 -17.6433775 -6.2433775
36 -8.6171082 -17.6433775
37 -10.3169426 -8.6171082
38 -16.7169426 -10.3169426
39 -15.3169426 -16.7169426
40 -20.3169426 -15.3169426
41 -33.3169426 -20.3169426
42 -28.5169426 -33.3169426
43 -27.9169426 -28.5169426
44 -42.1169426 -27.9169426
45 -4.0301876 -42.1169426
46 -5.6301876 -4.0301876
47 -4.0301876 -5.6301876
48 3.9960817 -4.0301876
49 5.2962472 3.9960817
50 0.8962472 5.2962472
51 5.2962472 0.8962472
52 -5.7037528 5.2962472
53 2.2962472 -5.7037528
54 1.0962472 2.2962472
55 -2.3037528 1.0962472
56 -9.5037528 -2.3037528
57 -15.5254415 -9.5037528
58 -4.1254415 -15.5254415
59 7.4745585 -4.1254415
60 24.5008278 7.4745585
61 NA 24.5008278
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -26.8311810 -36.1313466
[2,] -22.2311810 -26.8311810
[3,] -25.8311810 -22.2311810
[4,] -22.8311810 -25.8311810
[5,] -17.8311810 -22.8311810
[6,] -18.0311810 -17.8311810
[7,] -10.4311810 -18.0311810
[8,] 0.3688190 -10.4311810
[9,] 4.3471302 0.3688190
[10,] 2.7471302 4.3471302
[11,] 1.3471302 2.7471302
[12,] 2.3733996 1.3471302
[13,] 10.6735651 2.3733996
[14,] 16.2735651 10.6735651
[15,] 12.6735651 16.2735651
[16,] 16.6735651 12.6735651
[17,] 17.6735651 16.6735651
[18,] 18.4735651 17.6735651
[19,] 17.0735651 18.4735651
[20,] 24.8735651 17.0735651
[21,] 14.8518764 24.8735651
[22,] 13.2518764 14.8518764
[23,] 12.8518764 13.2518764
[24,] 13.8781457 12.8518764
[25,] 21.1783113 13.8781457
[26,] 21.7783113 21.1783113
[27,] 23.1783113 21.7783113
[28,] 32.1783113 23.1783113
[29,] 31.1783113 32.1783113
[30,] 26.9783113 31.1783113
[31,] 23.5783113 26.9783113
[32,] 26.3783113 23.5783113
[33,] 0.3566225 26.3783113
[34,] -6.2433775 0.3566225
[35,] -17.6433775 -6.2433775
[36,] -8.6171082 -17.6433775
[37,] -10.3169426 -8.6171082
[38,] -16.7169426 -10.3169426
[39,] -15.3169426 -16.7169426
[40,] -20.3169426 -15.3169426
[41,] -33.3169426 -20.3169426
[42,] -28.5169426 -33.3169426
[43,] -27.9169426 -28.5169426
[44,] -42.1169426 -27.9169426
[45,] -4.0301876 -42.1169426
[46,] -5.6301876 -4.0301876
[47,] -4.0301876 -5.6301876
[48,] 3.9960817 -4.0301876
[49,] 5.2962472 3.9960817
[50,] 0.8962472 5.2962472
[51,] 5.2962472 0.8962472
[52,] -5.7037528 5.2962472
[53,] 2.2962472 -5.7037528
[54,] 1.0962472 2.2962472
[55,] -2.3037528 1.0962472
[56,] -9.5037528 -2.3037528
[57,] -15.5254415 -9.5037528
[58,] -4.1254415 -15.5254415
[59,] 7.4745585 -4.1254415
[60,] 24.5008278 7.4745585
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -26.8311810 -36.1313466
2 -22.2311810 -26.8311810
3 -25.8311810 -22.2311810
4 -22.8311810 -25.8311810
5 -17.8311810 -22.8311810
6 -18.0311810 -17.8311810
7 -10.4311810 -18.0311810
8 0.3688190 -10.4311810
9 4.3471302 0.3688190
10 2.7471302 4.3471302
11 1.3471302 2.7471302
12 2.3733996 1.3471302
13 10.6735651 2.3733996
14 16.2735651 10.6735651
15 12.6735651 16.2735651
16 16.6735651 12.6735651
17 17.6735651 16.6735651
18 18.4735651 17.6735651
19 17.0735651 18.4735651
20 24.8735651 17.0735651
21 14.8518764 24.8735651
22 13.2518764 14.8518764
23 12.8518764 13.2518764
24 13.8781457 12.8518764
25 21.1783113 13.8781457
26 21.7783113 21.1783113
27 23.1783113 21.7783113
28 32.1783113 23.1783113
29 31.1783113 32.1783113
30 26.9783113 31.1783113
31 23.5783113 26.9783113
32 26.3783113 23.5783113
33 0.3566225 26.3783113
34 -6.2433775 0.3566225
35 -17.6433775 -6.2433775
36 -8.6171082 -17.6433775
37 -10.3169426 -8.6171082
38 -16.7169426 -10.3169426
39 -15.3169426 -16.7169426
40 -20.3169426 -15.3169426
41 -33.3169426 -20.3169426
42 -28.5169426 -33.3169426
43 -27.9169426 -28.5169426
44 -42.1169426 -27.9169426
45 -4.0301876 -42.1169426
46 -5.6301876 -4.0301876
47 -4.0301876 -5.6301876
48 3.9960817 -4.0301876
49 5.2962472 3.9960817
50 0.8962472 5.2962472
51 5.2962472 0.8962472
52 -5.7037528 5.2962472
53 2.2962472 -5.7037528
54 1.0962472 2.2962472
55 -2.3037528 1.0962472
56 -9.5037528 -2.3037528
57 -15.5254415 -9.5037528
58 -4.1254415 -15.5254415
59 7.4745585 -4.1254415
60 24.5008278 7.4745585
> 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/7pbtr1258570538.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/8vr0s1258570538.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/9c98x1258570538.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/109glc1258570538.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/114bk81258570538.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/12fp6a1258570538.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/13qtyh1258570538.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/14ynwd1258570538.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/15hoy81258570538.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/161vjk1258570538.tab")
+ }
>
> system("convert tmp/13hhj1258570537.ps tmp/13hhj1258570537.png")
> system("convert tmp/2m7kl1258570537.ps tmp/2m7kl1258570537.png")
> system("convert tmp/3vzkc1258570537.ps tmp/3vzkc1258570537.png")
> system("convert tmp/4k31b1258570537.ps tmp/4k31b1258570537.png")
> system("convert tmp/54wwr1258570537.ps tmp/54wwr1258570537.png")
> system("convert tmp/666as1258570537.ps tmp/666as1258570537.png")
> system("convert tmp/7pbtr1258570538.ps tmp/7pbtr1258570538.png")
> system("convert tmp/8vr0s1258570538.ps tmp/8vr0s1258570538.png")
> system("convert tmp/9c98x1258570538.ps tmp/9c98x1258570538.png")
> system("convert tmp/109glc1258570538.ps tmp/109glc1258570538.png")
>
>
> proc.time()
user system elapsed
2.406 1.586 3.201