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,13.9,561,15.9,555,18.2,544,19.7,537,20.1,543,19.9,594,20,611,22.6,613,20.6,611,20.1,594,20.2,595,21.8,591,22,589,19.5,584,17.5,573,18.2,567,18.8,569,19.7,621,18.8,629,18.5,628,18.7,612,18.5,595,19.3,597,18.9,593,21.4,590,22.5,580,25,574,22.9,573,22.9,573,21.3,620,22.3,626,20.9,620,19.9,588,20.2,566,19.8,557,17.7,561,18.1,549,17.6,532,18.2,526,16,511,16.3,499,17.3,555,19,565,18.6,542,18,527,17.9,510,17.8,514,18.5,517,17.4,508,19,493,17.4,490,20.6,469,18.5,478,20,528,18.8,534,18.8,518,19.7,506,15.3,502,10.6,516,6.1,528,0.9),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 = '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 562 13.9 1 0 0 0 0 0 0 0 0 0 0
2 561 15.9 0 1 0 0 0 0 0 0 0 0 0
3 555 18.2 0 0 1 0 0 0 0 0 0 0 0
4 544 19.7 0 0 0 1 0 0 0 0 0 0 0
5 537 20.1 0 0 0 0 1 0 0 0 0 0 0
6 543 19.9 0 0 0 0 0 1 0 0 0 0 0
7 594 20.0 0 0 0 0 0 0 1 0 0 0 0
8 611 22.6 0 0 0 0 0 0 0 1 0 0 0
9 613 20.6 0 0 0 0 0 0 0 0 1 0 0
10 611 20.1 0 0 0 0 0 0 0 0 0 1 0
11 594 20.2 0 0 0 0 0 0 0 0 0 0 1
12 595 21.8 0 0 0 0 0 0 0 0 0 0 0
13 591 22.0 1 0 0 0 0 0 0 0 0 0 0
14 589 19.5 0 1 0 0 0 0 0 0 0 0 0
15 584 17.5 0 0 1 0 0 0 0 0 0 0 0
16 573 18.2 0 0 0 1 0 0 0 0 0 0 0
17 567 18.8 0 0 0 0 1 0 0 0 0 0 0
18 569 19.7 0 0 0 0 0 1 0 0 0 0 0
19 621 18.8 0 0 0 0 0 0 1 0 0 0 0
20 629 18.5 0 0 0 0 0 0 0 1 0 0 0
21 628 18.7 0 0 0 0 0 0 0 0 1 0 0
22 612 18.5 0 0 0 0 0 0 0 0 0 1 0
23 595 19.3 0 0 0 0 0 0 0 0 0 0 1
24 597 18.9 0 0 0 0 0 0 0 0 0 0 0
25 593 21.4 1 0 0 0 0 0 0 0 0 0 0
26 590 22.5 0 1 0 0 0 0 0 0 0 0 0
27 580 25.0 0 0 1 0 0 0 0 0 0 0 0
28 574 22.9 0 0 0 1 0 0 0 0 0 0 0
29 573 22.9 0 0 0 0 1 0 0 0 0 0 0
30 573 21.3 0 0 0 0 0 1 0 0 0 0 0
31 620 22.3 0 0 0 0 0 0 1 0 0 0 0
32 626 20.9 0 0 0 0 0 0 0 1 0 0 0
33 620 19.9 0 0 0 0 0 0 0 0 1 0 0
34 588 20.2 0 0 0 0 0 0 0 0 0 1 0
35 566 19.8 0 0 0 0 0 0 0 0 0 0 1
36 557 17.7 0 0 0 0 0 0 0 0 0 0 0
37 561 18.1 1 0 0 0 0 0 0 0 0 0 0
38 549 17.6 0 1 0 0 0 0 0 0 0 0 0
39 532 18.2 0 0 1 0 0 0 0 0 0 0 0
40 526 16.0 0 0 0 1 0 0 0 0 0 0 0
41 511 16.3 0 0 0 0 1 0 0 0 0 0 0
42 499 17.3 0 0 0 0 0 1 0 0 0 0 0
43 555 19.0 0 0 0 0 0 0 1 0 0 0 0
44 565 18.6 0 0 0 0 0 0 0 1 0 0 0
45 542 18.0 0 0 0 0 0 0 0 0 1 0 0
46 527 17.9 0 0 0 0 0 0 0 0 0 1 0
47 510 17.8 0 0 0 0 0 0 0 0 0 0 1
48 514 18.5 0 0 0 0 0 0 0 0 0 0 0
49 517 17.4 1 0 0 0 0 0 0 0 0 0 0
50 508 19.0 0 1 0 0 0 0 0 0 0 0 0
51 493 17.4 0 0 1 0 0 0 0 0 0 0 0
52 490 20.6 0 0 0 1 0 0 0 0 0 0 0
53 469 18.5 0 0 0 0 1 0 0 0 0 0 0
54 478 20.0 0 0 0 0 0 1 0 0 0 0 0
55 528 18.8 0 0 0 0 0 0 1 0 0 0 0
56 534 18.8 0 0 0 0 0 0 0 1 0 0 0
57 518 19.7 0 0 0 0 0 0 0 0 1 0 0
58 506 15.3 0 0 0 0 0 0 0 0 0 1 0
59 502 10.6 0 0 0 0 0 0 0 0 0 0 1
60 516 6.1 0 0 0 0 0 0 0 0 0 0 0
61 528 0.9 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) X M1 M2 M3 M4
471.932 5.052 7.835 -8.020 -20.439 -28.951
M5 M6 M7 M8 M9 M10
-38.142 -38.759 11.734 20.629 14.355 3.906
M11
-7.149
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-67.817 -21.578 5.114 27.161 47.236
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 471.932 28.719 16.433 < 2e-16 ***
X 5.052 1.415 3.569 0.000824 ***
M1 7.835 22.405 0.350 0.728101
M2 -8.020 23.582 -0.340 0.735261
M3 -20.439 23.657 -0.864 0.391907
M4 -28.951 23.709 -1.221 0.228019
M5 -38.142 23.671 -1.611 0.113662
M6 -38.759 23.749 -1.632 0.109218
M7 11.734 23.786 0.493 0.624042
M8 20.629 23.813 0.866 0.390648
M9 14.355 23.685 0.606 0.547328
M10 3.906 23.494 0.166 0.868661
M11 -7.149 23.394 -0.306 0.761230
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 36.93 on 48 degrees of freedom
Multiple R-squared: 0.3761, Adjusted R-squared: 0.2201
F-statistic: 2.411 on 12 and 48 DF, p-value: 0.01551
> 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.1264859148 0.2529718296 0.87351409
[2,] 0.0980772438 0.1961544876 0.90192276
[3,] 0.0632981031 0.1265962062 0.93670190
[4,] 0.0476889998 0.0953779996 0.95231100
[5,] 0.0328298486 0.0656596972 0.96717015
[6,] 0.0215490783 0.0430981566 0.97845092
[7,] 0.0125445940 0.0250891879 0.98745541
[8,] 0.0065335810 0.0130671620 0.99346642
[9,] 0.0032124704 0.0064249408 0.99678753
[10,] 0.0015272222 0.0030544443 0.99847278
[11,] 0.0007452624 0.0014905248 0.99925474
[12,] 0.0003208699 0.0006417398 0.99967913
[13,] 0.0001689310 0.0003378620 0.99983107
[14,] 0.0001423248 0.0002846496 0.99985768
[15,] 0.0001762066 0.0003524132 0.99982379
[16,] 0.0001866880 0.0003733760 0.99981331
[17,] 0.0002592077 0.0005184154 0.99974079
[18,] 0.0011629340 0.0023258681 0.99883707
[19,] 0.0061190828 0.0122381656 0.99388092
[20,] 0.0250412081 0.0500824163 0.97495879
[21,] 0.0611072935 0.1222145871 0.93889271
[22,] 0.1096657265 0.2193314530 0.89033427
[23,] 0.1759090465 0.3518180931 0.82409095
[24,] 0.3127884332 0.6255768664 0.68721157
[25,] 0.3760675174 0.7521350349 0.62393248
[26,] 0.5905882562 0.8188234877 0.40941174
[27,] 0.6546092458 0.6907815085 0.34539075
[28,] 0.7436135039 0.5127729922 0.25638650
[29,] 0.8499643762 0.3000712476 0.15003562
[30,] 0.9003556860 0.1992886281 0.09964431
> postscript(file="/var/www/html/rcomp/tmp/1oodq1258757559.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/27mqb1258757559.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/3rj741258757559.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/4lk7m1258757559.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/5lbk51258757559.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
12.00642257 16.75685498 11.55542209 1.48849730 1.65921771 9.28640590
7 8 9 10 11 12
9.28849730 4.25778482 22.63621231 33.61111551 27.16092192 12.92811804
13 14 15 16 17 18
0.08291413 26.56862900 44.09202159 38.06692479 38.22718820 36.29686290
19 20 21 22 23 24
42.35123929 42.97215329 47.23555379 42.69477150 32.70797841 29.57974452
25 26 27 28 29 30
5.11428513 12.41177403 2.19988414 15.32118533 23.51281973 32.21320691
31 32 33 34 35 36
23.66824182 27.84666931 33.17281180 10.10588701 1.18183592 -4.35751349
37 38 39 40 41 42
-10.21317440 -3.83202951 -11.44457791 2.18195177 -5.14209932 -21.57765312
43 44 45 46 47 48
-24.65921771 -21.53307521 -35.22784671 -39.27385750 -44.71359410 -51.39934149
49 50 51 52 53 54
-50.67657490 -51.90522850 -46.40274991 -57.05855919 -58.25712631 -56.21882260
55 56 57 58 59 60
-50.64876071 -53.54353221 -67.81673120 -47.13791652 -16.33714215 13.24899242
61
43.68612747
> postscript(file="/var/www/html/rcomp/tmp/63nch1258757559.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 12.00642257 NA
1 16.75685498 12.00642257
2 11.55542209 16.75685498
3 1.48849730 11.55542209
4 1.65921771 1.48849730
5 9.28640590 1.65921771
6 9.28849730 9.28640590
7 4.25778482 9.28849730
8 22.63621231 4.25778482
9 33.61111551 22.63621231
10 27.16092192 33.61111551
11 12.92811804 27.16092192
12 0.08291413 12.92811804
13 26.56862900 0.08291413
14 44.09202159 26.56862900
15 38.06692479 44.09202159
16 38.22718820 38.06692479
17 36.29686290 38.22718820
18 42.35123929 36.29686290
19 42.97215329 42.35123929
20 47.23555379 42.97215329
21 42.69477150 47.23555379
22 32.70797841 42.69477150
23 29.57974452 32.70797841
24 5.11428513 29.57974452
25 12.41177403 5.11428513
26 2.19988414 12.41177403
27 15.32118533 2.19988414
28 23.51281973 15.32118533
29 32.21320691 23.51281973
30 23.66824182 32.21320691
31 27.84666931 23.66824182
32 33.17281180 27.84666931
33 10.10588701 33.17281180
34 1.18183592 10.10588701
35 -4.35751349 1.18183592
36 -10.21317440 -4.35751349
37 -3.83202951 -10.21317440
38 -11.44457791 -3.83202951
39 2.18195177 -11.44457791
40 -5.14209932 2.18195177
41 -21.57765312 -5.14209932
42 -24.65921771 -21.57765312
43 -21.53307521 -24.65921771
44 -35.22784671 -21.53307521
45 -39.27385750 -35.22784671
46 -44.71359410 -39.27385750
47 -51.39934149 -44.71359410
48 -50.67657490 -51.39934149
49 -51.90522850 -50.67657490
50 -46.40274991 -51.90522850
51 -57.05855919 -46.40274991
52 -58.25712631 -57.05855919
53 -56.21882260 -58.25712631
54 -50.64876071 -56.21882260
55 -53.54353221 -50.64876071
56 -67.81673120 -53.54353221
57 -47.13791652 -67.81673120
58 -16.33714215 -47.13791652
59 13.24899242 -16.33714215
60 43.68612747 13.24899242
61 NA 43.68612747
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 16.75685498 12.00642257
[2,] 11.55542209 16.75685498
[3,] 1.48849730 11.55542209
[4,] 1.65921771 1.48849730
[5,] 9.28640590 1.65921771
[6,] 9.28849730 9.28640590
[7,] 4.25778482 9.28849730
[8,] 22.63621231 4.25778482
[9,] 33.61111551 22.63621231
[10,] 27.16092192 33.61111551
[11,] 12.92811804 27.16092192
[12,] 0.08291413 12.92811804
[13,] 26.56862900 0.08291413
[14,] 44.09202159 26.56862900
[15,] 38.06692479 44.09202159
[16,] 38.22718820 38.06692479
[17,] 36.29686290 38.22718820
[18,] 42.35123929 36.29686290
[19,] 42.97215329 42.35123929
[20,] 47.23555379 42.97215329
[21,] 42.69477150 47.23555379
[22,] 32.70797841 42.69477150
[23,] 29.57974452 32.70797841
[24,] 5.11428513 29.57974452
[25,] 12.41177403 5.11428513
[26,] 2.19988414 12.41177403
[27,] 15.32118533 2.19988414
[28,] 23.51281973 15.32118533
[29,] 32.21320691 23.51281973
[30,] 23.66824182 32.21320691
[31,] 27.84666931 23.66824182
[32,] 33.17281180 27.84666931
[33,] 10.10588701 33.17281180
[34,] 1.18183592 10.10588701
[35,] -4.35751349 1.18183592
[36,] -10.21317440 -4.35751349
[37,] -3.83202951 -10.21317440
[38,] -11.44457791 -3.83202951
[39,] 2.18195177 -11.44457791
[40,] -5.14209932 2.18195177
[41,] -21.57765312 -5.14209932
[42,] -24.65921771 -21.57765312
[43,] -21.53307521 -24.65921771
[44,] -35.22784671 -21.53307521
[45,] -39.27385750 -35.22784671
[46,] -44.71359410 -39.27385750
[47,] -51.39934149 -44.71359410
[48,] -50.67657490 -51.39934149
[49,] -51.90522850 -50.67657490
[50,] -46.40274991 -51.90522850
[51,] -57.05855919 -46.40274991
[52,] -58.25712631 -57.05855919
[53,] -56.21882260 -58.25712631
[54,] -50.64876071 -56.21882260
[55,] -53.54353221 -50.64876071
[56,] -67.81673120 -53.54353221
[57,] -47.13791652 -67.81673120
[58,] -16.33714215 -47.13791652
[59,] 13.24899242 -16.33714215
[60,] 43.68612747 13.24899242
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 16.75685498 12.00642257
2 11.55542209 16.75685498
3 1.48849730 11.55542209
4 1.65921771 1.48849730
5 9.28640590 1.65921771
6 9.28849730 9.28640590
7 4.25778482 9.28849730
8 22.63621231 4.25778482
9 33.61111551 22.63621231
10 27.16092192 33.61111551
11 12.92811804 27.16092192
12 0.08291413 12.92811804
13 26.56862900 0.08291413
14 44.09202159 26.56862900
15 38.06692479 44.09202159
16 38.22718820 38.06692479
17 36.29686290 38.22718820
18 42.35123929 36.29686290
19 42.97215329 42.35123929
20 47.23555379 42.97215329
21 42.69477150 47.23555379
22 32.70797841 42.69477150
23 29.57974452 32.70797841
24 5.11428513 29.57974452
25 12.41177403 5.11428513
26 2.19988414 12.41177403
27 15.32118533 2.19988414
28 23.51281973 15.32118533
29 32.21320691 23.51281973
30 23.66824182 32.21320691
31 27.84666931 23.66824182
32 33.17281180 27.84666931
33 10.10588701 33.17281180
34 1.18183592 10.10588701
35 -4.35751349 1.18183592
36 -10.21317440 -4.35751349
37 -3.83202951 -10.21317440
38 -11.44457791 -3.83202951
39 2.18195177 -11.44457791
40 -5.14209932 2.18195177
41 -21.57765312 -5.14209932
42 -24.65921771 -21.57765312
43 -21.53307521 -24.65921771
44 -35.22784671 -21.53307521
45 -39.27385750 -35.22784671
46 -44.71359410 -39.27385750
47 -51.39934149 -44.71359410
48 -50.67657490 -51.39934149
49 -51.90522850 -50.67657490
50 -46.40274991 -51.90522850
51 -57.05855919 -46.40274991
52 -58.25712631 -57.05855919
53 -56.21882260 -58.25712631
54 -50.64876071 -56.21882260
55 -53.54353221 -50.64876071
56 -67.81673120 -53.54353221
57 -47.13791652 -67.81673120
58 -16.33714215 -47.13791652
59 13.24899242 -16.33714215
60 43.68612747 13.24899242
> 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/7hxu91258757559.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/883gr1258757559.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/9j8oy1258757559.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/10n29a1258757559.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/114ho21258757559.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/127d6k1258757559.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/13gisc1258757559.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/14wacl1258757559.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/15xywe1258757559.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/16etaf1258757559.tab")
+ }
>
> system("convert tmp/1oodq1258757559.ps tmp/1oodq1258757559.png")
> system("convert tmp/27mqb1258757559.ps tmp/27mqb1258757559.png")
> system("convert tmp/3rj741258757559.ps tmp/3rj741258757559.png")
> system("convert tmp/4lk7m1258757559.ps tmp/4lk7m1258757559.png")
> system("convert tmp/5lbk51258757559.ps tmp/5lbk51258757559.png")
> system("convert tmp/63nch1258757559.ps tmp/63nch1258757559.png")
> system("convert tmp/7hxu91258757559.ps tmp/7hxu91258757559.png")
> system("convert tmp/883gr1258757559.ps tmp/883gr1258757559.png")
> system("convert tmp/9j8oy1258757559.ps tmp/9j8oy1258757559.png")
> system("convert tmp/10n29a1258757559.ps tmp/10n29a1258757559.png")
>
>
> proc.time()
user system elapsed
2.357 1.560 2.983