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(19,595,19,18,19,22,591,19,19,18,23,589,22,19,19,20,584,23,22,19,14,573,20,23,22,14,567,14,20,23,14,569,14,14,20,15,621,14,14,14,11,629,15,14,14,17,628,11,15,14,16,612,17,11,15,20,595,16,17,11,24,597,20,16,17,23,593,24,20,16,20,590,23,24,20,21,580,20,23,24,19,574,21,20,23,23,573,19,21,20,23,573,23,19,21,23,620,23,23,19,23,626,23,23,23,27,620,23,23,23,26,588,27,23,23,17,566,26,27,23,24,557,17,26,27,26,561,24,17,26,24,549,26,24,17,27,532,24,26,24,27,526,27,24,26,26,511,27,27,24,24,499,26,27,27,23,555,24,26,27,23,565,23,24,26,24,542,23,23,24,17,527,24,23,23,21,510,17,24,23,19,514,21,17,24,22,517,19,21,17,22,508,22,19,21,18,493,22,22,19,16,490,18,22,22,14,469,16,18,22,12,478,14,16,18,14,528,12,14,16,16,534,14,12,14,8,518,16,14,12,3,506,8,16,14,0,502,3,8,16,5,516,0,3,8,1,528,5,0,3,1,533,1,5,0,3,536,1,1,5,6,537,3,1,1,7,524,6,3,1,8,536,7,6,3,14,587,8,7,6,14,597,14,8,7,13,581,14,14,8),dim=c(5,58),dimnames=list(c('y','x','y1','y2','y3'),1:58))
> y <- array(NA,dim=c(5,58),dimnames=list(c('y','x','y1','y2','y3'),1:58))
> 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 y1 y2 y3 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 19 595 19 18 19 1 0 0 0 0 0 0 0 0 0 0 1
2 22 591 19 19 18 0 1 0 0 0 0 0 0 0 0 0 2
3 23 589 22 19 19 0 0 1 0 0 0 0 0 0 0 0 3
4 20 584 23 22 19 0 0 0 1 0 0 0 0 0 0 0 4
5 14 573 20 23 22 0 0 0 0 1 0 0 0 0 0 0 5
6 14 567 14 20 23 0 0 0 0 0 1 0 0 0 0 0 6
7 14 569 14 14 20 0 0 0 0 0 0 1 0 0 0 0 7
8 15 621 14 14 14 0 0 0 0 0 0 0 1 0 0 0 8
9 11 629 15 14 14 0 0 0 0 0 0 0 0 1 0 0 9
10 17 628 11 15 14 0 0 0 0 0 0 0 0 0 1 0 10
11 16 612 17 11 15 0 0 0 0 0 0 0 0 0 0 1 11
12 20 595 16 17 11 0 0 0 0 0 0 0 0 0 0 0 12
13 24 597 20 16 17 1 0 0 0 0 0 0 0 0 0 0 13
14 23 593 24 20 16 0 1 0 0 0 0 0 0 0 0 0 14
15 20 590 23 24 20 0 0 1 0 0 0 0 0 0 0 0 15
16 21 580 20 23 24 0 0 0 1 0 0 0 0 0 0 0 16
17 19 574 21 20 23 0 0 0 0 1 0 0 0 0 0 0 17
18 23 573 19 21 20 0 0 0 0 0 1 0 0 0 0 0 18
19 23 573 23 19 21 0 0 0 0 0 0 1 0 0 0 0 19
20 23 620 23 23 19 0 0 0 0 0 0 0 1 0 0 0 20
21 23 626 23 23 23 0 0 0 0 0 0 0 0 1 0 0 21
22 27 620 23 23 23 0 0 0 0 0 0 0 0 0 1 0 22
23 26 588 27 23 23 0 0 0 0 0 0 0 0 0 0 1 23
24 17 566 26 27 23 0 0 0 0 0 0 0 0 0 0 0 24
25 24 557 17 26 27 1 0 0 0 0 0 0 0 0 0 0 25
26 26 561 24 17 26 0 1 0 0 0 0 0 0 0 0 0 26
27 24 549 26 24 17 0 0 1 0 0 0 0 0 0 0 0 27
28 27 532 24 26 24 0 0 0 1 0 0 0 0 0 0 0 28
29 27 526 27 24 26 0 0 0 0 1 0 0 0 0 0 0 29
30 26 511 27 27 24 0 0 0 0 0 1 0 0 0 0 0 30
31 24 499 26 27 27 0 0 0 0 0 0 1 0 0 0 0 31
32 23 555 24 26 27 0 0 0 0 0 0 0 1 0 0 0 32
33 23 565 23 24 26 0 0 0 0 0 0 0 0 1 0 0 33
34 24 542 23 23 24 0 0 0 0 0 0 0 0 0 1 0 34
35 17 527 24 23 23 0 0 0 0 0 0 0 0 0 0 1 35
36 21 510 17 24 23 0 0 0 0 0 0 0 0 0 0 0 36
37 19 514 21 17 24 1 0 0 0 0 0 0 0 0 0 0 37
38 22 517 19 21 17 0 1 0 0 0 0 0 0 0 0 0 38
39 22 508 22 19 21 0 0 1 0 0 0 0 0 0 0 0 39
40 18 493 22 22 19 0 0 0 1 0 0 0 0 0 0 0 40
41 16 490 18 22 22 0 0 0 0 1 0 0 0 0 0 0 41
42 14 469 16 18 22 0 0 0 0 0 1 0 0 0 0 0 42
43 12 478 14 16 18 0 0 0 0 0 0 1 0 0 0 0 43
44 14 528 12 14 16 0 0 0 0 0 0 0 1 0 0 0 44
45 16 534 14 12 14 0 0 0 0 0 0 0 0 1 0 0 45
46 8 518 16 14 12 0 0 0 0 0 0 0 0 0 1 0 46
47 3 506 8 16 14 0 0 0 0 0 0 0 0 0 0 1 47
48 0 502 3 8 16 0 0 0 0 0 0 0 0 0 0 0 48
49 5 516 0 3 8 1 0 0 0 0 0 0 0 0 0 0 49
50 1 528 5 0 3 0 1 0 0 0 0 0 0 0 0 0 50
51 1 533 1 5 0 0 0 1 0 0 0 0 0 0 0 0 51
52 3 536 1 1 5 0 0 0 1 0 0 0 0 0 0 0 52
53 6 537 3 1 1 0 0 0 0 1 0 0 0 0 0 0 53
54 7 524 6 3 1 0 0 0 0 0 1 0 0 0 0 0 54
55 8 536 7 6 3 0 0 0 0 0 0 1 0 0 0 0 55
56 14 587 8 7 6 0 0 0 0 0 0 0 1 0 0 0 56
57 14 597 14 8 7 0 0 0 0 0 0 0 0 1 0 0 57
58 13 581 14 14 8 0 0 0 0 0 0 0 0 0 1 0 58
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x y1 y2 y3 M1
-40.29448 0.06385 0.68300 0.12246 0.21821 3.70568
M2 M3 M4 M5 M6 M7
2.87392 1.61334 1.73017 0.64607 2.80949 1.90389
M8 M9 M10 M11 t
0.79388 -1.33057 -0.03759 -2.25856 0.10701
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-7.4983 -1.7986 0.1300 2.0640 5.6072
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -40.29448 15.29372 -2.635 0.01183 *
x 0.06385 0.02315 2.759 0.00863 **
y1 0.68300 0.14992 4.556 4.62e-05 ***
y2 0.12246 0.18517 0.661 0.51210
y3 0.21821 0.15998 1.364 0.18000
M1 3.70568 2.21962 1.670 0.10263
M2 2.87392 2.31071 1.244 0.22066
M3 1.61334 2.24578 0.718 0.47659
M4 1.73017 2.16577 0.799 0.42897
M5 0.64607 2.18984 0.295 0.76946
M6 2.80949 2.18680 1.285 0.20609
M7 1.90389 2.23637 0.851 0.39953
M8 0.79388 2.30301 0.345 0.73207
M9 -1.33057 2.44462 -0.544 0.58919
M10 -0.03759 2.30022 -0.016 0.98704
M11 -2.25856 2.32634 -0.971 0.33731
t 0.10701 0.05772 1.854 0.07095 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.168 on 41 degrees of freedom
Multiple R-squared: 0.8668, Adjusted R-squared: 0.8148
F-statistic: 16.67 on 16 and 41 DF, p-value: 4.493e-13
> 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.3126220 0.6252439 0.6873780
[2,] 0.2056306 0.4112611 0.7943694
[3,] 0.1046795 0.2093591 0.8953205
[4,] 0.2757124 0.5514247 0.7242876
[5,] 0.7654836 0.4690329 0.2345164
[6,] 0.7860686 0.4278628 0.2139314
[7,] 0.6968432 0.6063136 0.3031568
[8,] 0.6232606 0.7534788 0.3767394
[9,] 0.5905773 0.8188454 0.4094227
[10,] 0.5425391 0.9149217 0.4574609
[11,] 0.4415129 0.8830258 0.5584871
[12,] 0.3323132 0.6646265 0.6676868
[13,] 0.3921644 0.7843288 0.6078356
[14,] 0.6045167 0.7909666 0.3954833
[15,] 0.5006383 0.9987235 0.4993617
[16,] 0.5258566 0.9482868 0.4741434
[17,] 0.4211544 0.8423088 0.5788456
[18,] 0.8142260 0.3715480 0.1857740
[19,] 0.7355726 0.5288548 0.2644274
> postscript(file="/var/www/html/rcomp/tmp/1k3sh1258660208.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/2xv3h1258660208.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/370qe1258660208.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/4ik461258660208.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/5xbxk1258660208.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 = 58
Frequency = 1
1 2 3 4 5 6 7
-1.8386265 2.2373026 2.2513661 -1.7035762 -4.7521954 -2.3923434 -0.3320633
8 9 10 11 12 13 14
-0.3401749 -3.5165647 3.7568430 2.0660975 5.6071528 1.7479316 -1.2755106
15 16 17 18 19 20 21
-3.6100745 -0.8967761 -1.6339691 2.0576422 0.1509417 -1.9005924 -1.1391339
22 23 24 25 26 27 28
1.8440057 2.2693030 -7.4983139 1.6602755 0.6689632 0.3295038 3.7847617
29 30 31 32 33 34 35
2.9044681 0.6609054 0.2541016 -1.8302576 0.6947743 2.3223190 -2.0706892
36 37 38 39 40 41 42
5.3077935 -2.8533120 3.0835506 2.1348644 -1.0621098 0.1838967 -0.8897649
43 44 45 46 47 48 49
-0.1820833 1.6755690 4.6252343 -4.9275719 -2.2647113 -3.4166323 1.2837314
50 51 52 53 54 55 56
-4.7143058 -1.1056598 -0.1222996 3.2977996 0.5635607 0.1091033 2.3954559
57 58
-0.6643101 -2.9955958
> postscript(file="/var/www/html/rcomp/tmp/6vlon1258660208.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 = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 -1.8386265 NA
1 2.2373026 -1.8386265
2 2.2513661 2.2373026
3 -1.7035762 2.2513661
4 -4.7521954 -1.7035762
5 -2.3923434 -4.7521954
6 -0.3320633 -2.3923434
7 -0.3401749 -0.3320633
8 -3.5165647 -0.3401749
9 3.7568430 -3.5165647
10 2.0660975 3.7568430
11 5.6071528 2.0660975
12 1.7479316 5.6071528
13 -1.2755106 1.7479316
14 -3.6100745 -1.2755106
15 -0.8967761 -3.6100745
16 -1.6339691 -0.8967761
17 2.0576422 -1.6339691
18 0.1509417 2.0576422
19 -1.9005924 0.1509417
20 -1.1391339 -1.9005924
21 1.8440057 -1.1391339
22 2.2693030 1.8440057
23 -7.4983139 2.2693030
24 1.6602755 -7.4983139
25 0.6689632 1.6602755
26 0.3295038 0.6689632
27 3.7847617 0.3295038
28 2.9044681 3.7847617
29 0.6609054 2.9044681
30 0.2541016 0.6609054
31 -1.8302576 0.2541016
32 0.6947743 -1.8302576
33 2.3223190 0.6947743
34 -2.0706892 2.3223190
35 5.3077935 -2.0706892
36 -2.8533120 5.3077935
37 3.0835506 -2.8533120
38 2.1348644 3.0835506
39 -1.0621098 2.1348644
40 0.1838967 -1.0621098
41 -0.8897649 0.1838967
42 -0.1820833 -0.8897649
43 1.6755690 -0.1820833
44 4.6252343 1.6755690
45 -4.9275719 4.6252343
46 -2.2647113 -4.9275719
47 -3.4166323 -2.2647113
48 1.2837314 -3.4166323
49 -4.7143058 1.2837314
50 -1.1056598 -4.7143058
51 -0.1222996 -1.1056598
52 3.2977996 -0.1222996
53 0.5635607 3.2977996
54 0.1091033 0.5635607
55 2.3954559 0.1091033
56 -0.6643101 2.3954559
57 -2.9955958 -0.6643101
58 NA -2.9955958
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2.2373026 -1.8386265
[2,] 2.2513661 2.2373026
[3,] -1.7035762 2.2513661
[4,] -4.7521954 -1.7035762
[5,] -2.3923434 -4.7521954
[6,] -0.3320633 -2.3923434
[7,] -0.3401749 -0.3320633
[8,] -3.5165647 -0.3401749
[9,] 3.7568430 -3.5165647
[10,] 2.0660975 3.7568430
[11,] 5.6071528 2.0660975
[12,] 1.7479316 5.6071528
[13,] -1.2755106 1.7479316
[14,] -3.6100745 -1.2755106
[15,] -0.8967761 -3.6100745
[16,] -1.6339691 -0.8967761
[17,] 2.0576422 -1.6339691
[18,] 0.1509417 2.0576422
[19,] -1.9005924 0.1509417
[20,] -1.1391339 -1.9005924
[21,] 1.8440057 -1.1391339
[22,] 2.2693030 1.8440057
[23,] -7.4983139 2.2693030
[24,] 1.6602755 -7.4983139
[25,] 0.6689632 1.6602755
[26,] 0.3295038 0.6689632
[27,] 3.7847617 0.3295038
[28,] 2.9044681 3.7847617
[29,] 0.6609054 2.9044681
[30,] 0.2541016 0.6609054
[31,] -1.8302576 0.2541016
[32,] 0.6947743 -1.8302576
[33,] 2.3223190 0.6947743
[34,] -2.0706892 2.3223190
[35,] 5.3077935 -2.0706892
[36,] -2.8533120 5.3077935
[37,] 3.0835506 -2.8533120
[38,] 2.1348644 3.0835506
[39,] -1.0621098 2.1348644
[40,] 0.1838967 -1.0621098
[41,] -0.8897649 0.1838967
[42,] -0.1820833 -0.8897649
[43,] 1.6755690 -0.1820833
[44,] 4.6252343 1.6755690
[45,] -4.9275719 4.6252343
[46,] -2.2647113 -4.9275719
[47,] -3.4166323 -2.2647113
[48,] 1.2837314 -3.4166323
[49,] -4.7143058 1.2837314
[50,] -1.1056598 -4.7143058
[51,] -0.1222996 -1.1056598
[52,] 3.2977996 -0.1222996
[53,] 0.5635607 3.2977996
[54,] 0.1091033 0.5635607
[55,] 2.3954559 0.1091033
[56,] -0.6643101 2.3954559
[57,] -2.9955958 -0.6643101
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2.2373026 -1.8386265
2 2.2513661 2.2373026
3 -1.7035762 2.2513661
4 -4.7521954 -1.7035762
5 -2.3923434 -4.7521954
6 -0.3320633 -2.3923434
7 -0.3401749 -0.3320633
8 -3.5165647 -0.3401749
9 3.7568430 -3.5165647
10 2.0660975 3.7568430
11 5.6071528 2.0660975
12 1.7479316 5.6071528
13 -1.2755106 1.7479316
14 -3.6100745 -1.2755106
15 -0.8967761 -3.6100745
16 -1.6339691 -0.8967761
17 2.0576422 -1.6339691
18 0.1509417 2.0576422
19 -1.9005924 0.1509417
20 -1.1391339 -1.9005924
21 1.8440057 -1.1391339
22 2.2693030 1.8440057
23 -7.4983139 2.2693030
24 1.6602755 -7.4983139
25 0.6689632 1.6602755
26 0.3295038 0.6689632
27 3.7847617 0.3295038
28 2.9044681 3.7847617
29 0.6609054 2.9044681
30 0.2541016 0.6609054
31 -1.8302576 0.2541016
32 0.6947743 -1.8302576
33 2.3223190 0.6947743
34 -2.0706892 2.3223190
35 5.3077935 -2.0706892
36 -2.8533120 5.3077935
37 3.0835506 -2.8533120
38 2.1348644 3.0835506
39 -1.0621098 2.1348644
40 0.1838967 -1.0621098
41 -0.8897649 0.1838967
42 -0.1820833 -0.8897649
43 1.6755690 -0.1820833
44 4.6252343 1.6755690
45 -4.9275719 4.6252343
46 -2.2647113 -4.9275719
47 -3.4166323 -2.2647113
48 1.2837314 -3.4166323
49 -4.7143058 1.2837314
50 -1.1056598 -4.7143058
51 -0.1222996 -1.1056598
52 3.2977996 -0.1222996
53 0.5635607 3.2977996
54 0.1091033 0.5635607
55 2.3954559 0.1091033
56 -0.6643101 2.3954559
57 -2.9955958 -0.6643101
> 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/7llku1258660208.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/8a8tw1258660208.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/9lv0m1258660208.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/10pz941258660208.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/11io811258660208.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/12uioe1258660208.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/139qf41258660209.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/14qm931258660209.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/151y7t1258660209.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/162kth1258660209.tab")
+ }
>
> system("convert tmp/1k3sh1258660208.ps tmp/1k3sh1258660208.png")
> system("convert tmp/2xv3h1258660208.ps tmp/2xv3h1258660208.png")
> system("convert tmp/370qe1258660208.ps tmp/370qe1258660208.png")
> system("convert tmp/4ik461258660208.ps tmp/4ik461258660208.png")
> system("convert tmp/5xbxk1258660208.ps tmp/5xbxk1258660208.png")
> system("convert tmp/6vlon1258660208.ps tmp/6vlon1258660208.png")
> system("convert tmp/7llku1258660208.ps tmp/7llku1258660208.png")
> system("convert tmp/8a8tw1258660208.ps tmp/8a8tw1258660208.png")
> system("convert tmp/9lv0m1258660208.ps tmp/9lv0m1258660208.png")
> system("convert tmp/10pz941258660208.ps tmp/10pz941258660208.png")
>
>
> proc.time()
user system elapsed
2.408 1.617 5.718