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(449,0,452,0,462,0,455,0,461,0,461,0,463,0,462,0,456,0,455,0,456,0,472,0,472,0,471,0,465,0,459,0,465,0,468,0,467,0,463,0,460,0,462,0,461,0,476,0,476,0,471,0,453,0,443,0,442,0,444,0,438,0,427,0,424,0,416,0,406,0,431,0,434,0,418,0,412,0,404,0,409,0,412,1,406,1,398,1,397,1,385,1,390,1,413,1,413,1,401,1,397,1,397,1,409,1,419,1,424,1,428,1,430,1,424,1,433,1,456,1,459,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 449 0 1 0 0 0 0 0 0 0 0 0 0 1
2 452 0 0 1 0 0 0 0 0 0 0 0 0 2
3 462 0 0 0 1 0 0 0 0 0 0 0 0 3
4 455 0 0 0 0 1 0 0 0 0 0 0 0 4
5 461 0 0 0 0 0 1 0 0 0 0 0 0 5
6 461 0 0 0 0 0 0 1 0 0 0 0 0 6
7 463 0 0 0 0 0 0 0 1 0 0 0 0 7
8 462 0 0 0 0 0 0 0 0 1 0 0 0 8
9 456 0 0 0 0 0 0 0 0 0 1 0 0 9
10 455 0 0 0 0 0 0 0 0 0 0 1 0 10
11 456 0 0 0 0 0 0 0 0 0 0 0 1 11
12 472 0 0 0 0 0 0 0 0 0 0 0 0 12
13 472 0 1 0 0 0 0 0 0 0 0 0 0 13
14 471 0 0 1 0 0 0 0 0 0 0 0 0 14
15 465 0 0 0 1 0 0 0 0 0 0 0 0 15
16 459 0 0 0 0 1 0 0 0 0 0 0 0 16
17 465 0 0 0 0 0 1 0 0 0 0 0 0 17
18 468 0 0 0 0 0 0 1 0 0 0 0 0 18
19 467 0 0 0 0 0 0 0 1 0 0 0 0 19
20 463 0 0 0 0 0 0 0 0 1 0 0 0 20
21 460 0 0 0 0 0 0 0 0 0 1 0 0 21
22 462 0 0 0 0 0 0 0 0 0 0 1 0 22
23 461 0 0 0 0 0 0 0 0 0 0 0 1 23
24 476 0 0 0 0 0 0 0 0 0 0 0 0 24
25 476 0 1 0 0 0 0 0 0 0 0 0 0 25
26 471 0 0 1 0 0 0 0 0 0 0 0 0 26
27 453 0 0 0 1 0 0 0 0 0 0 0 0 27
28 443 0 0 0 0 1 0 0 0 0 0 0 0 28
29 442 0 0 0 0 0 1 0 0 0 0 0 0 29
30 444 0 0 0 0 0 0 1 0 0 0 0 0 30
31 438 0 0 0 0 0 0 0 1 0 0 0 0 31
32 427 0 0 0 0 0 0 0 0 1 0 0 0 32
33 424 0 0 0 0 0 0 0 0 0 1 0 0 33
34 416 0 0 0 0 0 0 0 0 0 0 1 0 34
35 406 0 0 0 0 0 0 0 0 0 0 0 1 35
36 431 0 0 0 0 0 0 0 0 0 0 0 0 36
37 434 0 1 0 0 0 0 0 0 0 0 0 0 37
38 418 0 0 1 0 0 0 0 0 0 0 0 0 38
39 412 0 0 0 1 0 0 0 0 0 0 0 0 39
40 404 0 0 0 0 1 0 0 0 0 0 0 0 40
41 409 0 0 0 0 0 1 0 0 0 0 0 0 41
42 412 1 0 0 0 0 0 1 0 0 0 0 0 42
43 406 1 0 0 0 0 0 0 1 0 0 0 0 43
44 398 1 0 0 0 0 0 0 0 1 0 0 0 44
45 397 1 0 0 0 0 0 0 0 0 1 0 0 45
46 385 1 0 0 0 0 0 0 0 0 0 1 0 46
47 390 1 0 0 0 0 0 0 0 0 0 0 1 47
48 413 1 0 0 0 0 0 0 0 0 0 0 0 48
49 413 1 1 0 0 0 0 0 0 0 0 0 0 49
50 401 1 0 1 0 0 0 0 0 0 0 0 0 50
51 397 1 0 0 1 0 0 0 0 0 0 0 0 51
52 397 1 0 0 0 1 0 0 0 0 0 0 0 52
53 409 1 0 0 0 0 1 0 0 0 0 0 0 53
54 419 1 0 0 0 0 0 1 0 0 0 0 0 54
55 424 1 0 0 0 0 0 0 1 0 0 0 0 55
56 428 1 0 0 0 0 0 0 0 1 0 0 0 56
57 430 1 0 0 0 0 0 0 0 0 1 0 0 57
58 424 1 0 0 0 0 0 0 0 0 0 1 0 58
59 433 1 0 0 0 0 0 0 0 0 0 0 1 59
60 456 1 0 0 0 0 0 0 0 0 0 0 0 60
61 459 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
483.6914 -11.4134 -3.9617 -17.4844 -21.4642 -26.8440
M5 M6 M7 M8 M9 M10
-20.4239 -13.7210 -14.1008 -17.2807 -18.6605 -22.8403
M11 t
-21.2202 -0.8202
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-29.910 -15.128 1.433 12.993 40.714
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 483.6914 10.4940 46.092 < 2e-16 ***
X -11.4134 9.0988 -1.254 0.21591
M1 -3.9617 11.5742 -0.342 0.73366
M2 -17.4844 12.1429 -1.440 0.15653
M3 -21.4642 12.1266 -1.770 0.08321 .
M4 -26.8440 12.1152 -2.216 0.03159 *
M5 -20.4239 12.1086 -1.687 0.09828 .
M6 -13.7210 12.1469 -1.130 0.26438
M7 -14.1008 12.1205 -1.163 0.25054
M8 -17.2807 12.0988 -1.428 0.15982
M9 -18.6605 12.0819 -1.544 0.12918
M10 -22.8403 12.0699 -1.892 0.06461 .
M11 -21.2202 12.0626 -1.759 0.08506 .
t -0.8202 0.2415 -3.396 0.00140 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.07 on 47 degrees of freedom
Multiple R-squared: 0.5881, Adjusted R-squared: 0.4742
F-statistic: 5.163 on 13 and 47 DF, p-value: 1.371e-05
> 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,] 5.676786e-02 1.135357e-01 0.94323214
[2,] 1.714225e-02 3.428450e-02 0.98285775
[3,] 5.470929e-03 1.094186e-02 0.99452907
[4,] 2.056541e-03 4.113081e-03 0.99794346
[5,] 5.918694e-04 1.183739e-03 0.99940813
[6,] 1.749696e-04 3.499393e-04 0.99982503
[7,] 5.698021e-05 1.139604e-04 0.99994302
[8,] 1.846289e-05 3.692578e-05 0.99998154
[9,] 7.782418e-06 1.556484e-05 0.99999222
[10,] 9.028600e-06 1.805720e-05 0.99999097
[11,] 6.334501e-04 1.266900e-03 0.99936655
[12,] 1.793694e-02 3.587389e-02 0.98206306
[13,] 3.474415e-01 6.948830e-01 0.65255850
[14,] 6.123636e-01 7.752727e-01 0.38763637
[15,] 8.383489e-01 3.233023e-01 0.16165114
[16,] 9.357964e-01 1.284071e-01 0.06420356
[17,] 9.616092e-01 7.678164e-02 0.03839082
[18,] 9.858943e-01 2.821133e-02 0.01410567
[19,] 9.899376e-01 2.012483e-02 0.01006241
[20,] 9.861887e-01 2.762257e-02 0.01381129
[21,] 9.743109e-01 5.137815e-02 0.02568907
[22,] 9.655085e-01 6.898299e-02 0.03449149
[23,] 9.528385e-01 9.432292e-02 0.04716146
[24,] 9.231330e-01 1.537339e-01 0.07686697
[25,] 8.675074e-01 2.649853e-01 0.13249265
[26,] 9.496942e-01 1.006115e-01 0.05030576
[27,] 9.819335e-01 3.613291e-02 0.01806646
[28,] 9.758139e-01 4.837210e-02 0.02418605
> postscript(file="/var/www/html/rcomp/tmp/13l6t1258660840.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/2vz2g1258660840.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/3l0661258660840.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/4vv5w1258660840.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/5orxq1258660840.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
-29.9095238 -12.5667262 2.2332738 1.4332738 1.8332738 -4.0494048
7 8 9 10 11 12
-0.8494048 2.1505952 -1.6494048 2.3505952 2.5505952 -1.8494048
13 14 15 16 17 18
2.9325000 16.2752976 15.0752976 15.2752976 15.6752976 12.7926190
19 20 21 22 23 24
12.9926190 12.9926190 12.1926190 19.1926190 17.3926190 11.9926190
25 26 27 28 29 30
16.7745238 26.1173214 12.9173214 9.1173214 2.5173214 -1.3653571
31 32 33 34 35 36
-6.1653571 -13.1653571 -13.9653571 -16.9653571 -27.7653571 -23.1653571
37 38 39 40 41 42
-15.3834524 -17.0406548 -18.2406548 -20.0406548 -20.6406548 -12.1099405
43 44 45 46 47 48
-16.9099405 -20.9099405 -19.7099405 -26.7099405 -22.5099405 -19.9099405
49 50 51 52 53 54
-15.1280357 -12.7852381 -11.9852381 -5.7852381 0.6147619 4.7320833
55 56 57 58 59 60
10.9320833 18.9320833 23.1320833 22.1320833 30.3320833 32.9320833
61
40.7139881
> postscript(file="/var/www/html/rcomp/tmp/6rhxn1258660840.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 -29.9095238 NA
1 -12.5667262 -29.9095238
2 2.2332738 -12.5667262
3 1.4332738 2.2332738
4 1.8332738 1.4332738
5 -4.0494048 1.8332738
6 -0.8494048 -4.0494048
7 2.1505952 -0.8494048
8 -1.6494048 2.1505952
9 2.3505952 -1.6494048
10 2.5505952 2.3505952
11 -1.8494048 2.5505952
12 2.9325000 -1.8494048
13 16.2752976 2.9325000
14 15.0752976 16.2752976
15 15.2752976 15.0752976
16 15.6752976 15.2752976
17 12.7926190 15.6752976
18 12.9926190 12.7926190
19 12.9926190 12.9926190
20 12.1926190 12.9926190
21 19.1926190 12.1926190
22 17.3926190 19.1926190
23 11.9926190 17.3926190
24 16.7745238 11.9926190
25 26.1173214 16.7745238
26 12.9173214 26.1173214
27 9.1173214 12.9173214
28 2.5173214 9.1173214
29 -1.3653571 2.5173214
30 -6.1653571 -1.3653571
31 -13.1653571 -6.1653571
32 -13.9653571 -13.1653571
33 -16.9653571 -13.9653571
34 -27.7653571 -16.9653571
35 -23.1653571 -27.7653571
36 -15.3834524 -23.1653571
37 -17.0406548 -15.3834524
38 -18.2406548 -17.0406548
39 -20.0406548 -18.2406548
40 -20.6406548 -20.0406548
41 -12.1099405 -20.6406548
42 -16.9099405 -12.1099405
43 -20.9099405 -16.9099405
44 -19.7099405 -20.9099405
45 -26.7099405 -19.7099405
46 -22.5099405 -26.7099405
47 -19.9099405 -22.5099405
48 -15.1280357 -19.9099405
49 -12.7852381 -15.1280357
50 -11.9852381 -12.7852381
51 -5.7852381 -11.9852381
52 0.6147619 -5.7852381
53 4.7320833 0.6147619
54 10.9320833 4.7320833
55 18.9320833 10.9320833
56 23.1320833 18.9320833
57 22.1320833 23.1320833
58 30.3320833 22.1320833
59 32.9320833 30.3320833
60 40.7139881 32.9320833
61 NA 40.7139881
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -12.5667262 -29.9095238
[2,] 2.2332738 -12.5667262
[3,] 1.4332738 2.2332738
[4,] 1.8332738 1.4332738
[5,] -4.0494048 1.8332738
[6,] -0.8494048 -4.0494048
[7,] 2.1505952 -0.8494048
[8,] -1.6494048 2.1505952
[9,] 2.3505952 -1.6494048
[10,] 2.5505952 2.3505952
[11,] -1.8494048 2.5505952
[12,] 2.9325000 -1.8494048
[13,] 16.2752976 2.9325000
[14,] 15.0752976 16.2752976
[15,] 15.2752976 15.0752976
[16,] 15.6752976 15.2752976
[17,] 12.7926190 15.6752976
[18,] 12.9926190 12.7926190
[19,] 12.9926190 12.9926190
[20,] 12.1926190 12.9926190
[21,] 19.1926190 12.1926190
[22,] 17.3926190 19.1926190
[23,] 11.9926190 17.3926190
[24,] 16.7745238 11.9926190
[25,] 26.1173214 16.7745238
[26,] 12.9173214 26.1173214
[27,] 9.1173214 12.9173214
[28,] 2.5173214 9.1173214
[29,] -1.3653571 2.5173214
[30,] -6.1653571 -1.3653571
[31,] -13.1653571 -6.1653571
[32,] -13.9653571 -13.1653571
[33,] -16.9653571 -13.9653571
[34,] -27.7653571 -16.9653571
[35,] -23.1653571 -27.7653571
[36,] -15.3834524 -23.1653571
[37,] -17.0406548 -15.3834524
[38,] -18.2406548 -17.0406548
[39,] -20.0406548 -18.2406548
[40,] -20.6406548 -20.0406548
[41,] -12.1099405 -20.6406548
[42,] -16.9099405 -12.1099405
[43,] -20.9099405 -16.9099405
[44,] -19.7099405 -20.9099405
[45,] -26.7099405 -19.7099405
[46,] -22.5099405 -26.7099405
[47,] -19.9099405 -22.5099405
[48,] -15.1280357 -19.9099405
[49,] -12.7852381 -15.1280357
[50,] -11.9852381 -12.7852381
[51,] -5.7852381 -11.9852381
[52,] 0.6147619 -5.7852381
[53,] 4.7320833 0.6147619
[54,] 10.9320833 4.7320833
[55,] 18.9320833 10.9320833
[56,] 23.1320833 18.9320833
[57,] 22.1320833 23.1320833
[58,] 30.3320833 22.1320833
[59,] 32.9320833 30.3320833
[60,] 40.7139881 32.9320833
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -12.5667262 -29.9095238
2 2.2332738 -12.5667262
3 1.4332738 2.2332738
4 1.8332738 1.4332738
5 -4.0494048 1.8332738
6 -0.8494048 -4.0494048
7 2.1505952 -0.8494048
8 -1.6494048 2.1505952
9 2.3505952 -1.6494048
10 2.5505952 2.3505952
11 -1.8494048 2.5505952
12 2.9325000 -1.8494048
13 16.2752976 2.9325000
14 15.0752976 16.2752976
15 15.2752976 15.0752976
16 15.6752976 15.2752976
17 12.7926190 15.6752976
18 12.9926190 12.7926190
19 12.9926190 12.9926190
20 12.1926190 12.9926190
21 19.1926190 12.1926190
22 17.3926190 19.1926190
23 11.9926190 17.3926190
24 16.7745238 11.9926190
25 26.1173214 16.7745238
26 12.9173214 26.1173214
27 9.1173214 12.9173214
28 2.5173214 9.1173214
29 -1.3653571 2.5173214
30 -6.1653571 -1.3653571
31 -13.1653571 -6.1653571
32 -13.9653571 -13.1653571
33 -16.9653571 -13.9653571
34 -27.7653571 -16.9653571
35 -23.1653571 -27.7653571
36 -15.3834524 -23.1653571
37 -17.0406548 -15.3834524
38 -18.2406548 -17.0406548
39 -20.0406548 -18.2406548
40 -20.6406548 -20.0406548
41 -12.1099405 -20.6406548
42 -16.9099405 -12.1099405
43 -20.9099405 -16.9099405
44 -19.7099405 -20.9099405
45 -26.7099405 -19.7099405
46 -22.5099405 -26.7099405
47 -19.9099405 -22.5099405
48 -15.1280357 -19.9099405
49 -12.7852381 -15.1280357
50 -11.9852381 -12.7852381
51 -5.7852381 -11.9852381
52 0.6147619 -5.7852381
53 4.7320833 0.6147619
54 10.9320833 4.7320833
55 18.9320833 10.9320833
56 23.1320833 18.9320833
57 22.1320833 23.1320833
58 30.3320833 22.1320833
59 32.9320833 30.3320833
60 40.7139881 32.9320833
> 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/7h10p1258660840.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/89cca1258660840.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/98zkv1258660840.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/10rvg41258660840.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/11g9s51258660840.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/12knnw1258660840.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/130klb1258660840.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/14vigv1258660840.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/15ea211258660840.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/165pw51258660840.tab")
+ }
>
> system("convert tmp/13l6t1258660840.ps tmp/13l6t1258660840.png")
> system("convert tmp/2vz2g1258660840.ps tmp/2vz2g1258660840.png")
> system("convert tmp/3l0661258660840.ps tmp/3l0661258660840.png")
> system("convert tmp/4vv5w1258660840.ps tmp/4vv5w1258660840.png")
> system("convert tmp/5orxq1258660840.ps tmp/5orxq1258660840.png")
> system("convert tmp/6rhxn1258660840.ps tmp/6rhxn1258660840.png")
> system("convert tmp/7h10p1258660840.ps tmp/7h10p1258660840.png")
> system("convert tmp/89cca1258660840.ps tmp/89cca1258660840.png")
> system("convert tmp/98zkv1258660840.ps tmp/98zkv1258660840.png")
> system("convert tmp/10rvg41258660840.ps tmp/10rvg41258660840.png")
>
>
> proc.time()
user system elapsed
2.451 1.591 4.847