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(17192.4,0,15386.1,0,14287.1,0,17526.6,0,14497,0,14398.3,0,16629.6,0,16670.7,0,16614.8,0,16869.2,0,15663.9,0,16359.9,0,18447.7,0,16889,0,16505,0,18320.9,0,15052.1,0,15699.8,0,18135.3,0,16768.7,0,18883,0,19021,0,18101.9,0,17776.1,0,21489.9,0,17065.3,0,18690,0,18953.1,0,16398.9,0,16895.6,0,18553,0,19270,0,19422.1,0,17579.4,0,18637.3,0,18076.7,0,20438.6,0,18075.2,0,19563,0,19899.2,0,19227.5,0,17789.6,0,19220.8,0,21968.9,0,21131.5,0,19484.6,0,22168.7,1,20866.8,1,22176.2,1,23533.8,1,21479.6,1,24347.7,1,22751.6,1,20328.3,1,23650.4,1,23335.7,1,19614.9,1,18042.3,1,17282.5,1,16847.2,1,18159.5,1),dim=c(2,61),dimnames=list(c('Invoer','Dummy'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Invoer','Dummy'),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
Invoer Dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 17192.4 0 1 0 0 0 0 0 0 0 0 0 0 1
2 15386.1 0 0 1 0 0 0 0 0 0 0 0 0 2
3 14287.1 0 0 0 1 0 0 0 0 0 0 0 0 3
4 17526.6 0 0 0 0 1 0 0 0 0 0 0 0 4
5 14497.0 0 0 0 0 0 1 0 0 0 0 0 0 5
6 14398.3 0 0 0 0 0 0 1 0 0 0 0 0 6
7 16629.6 0 0 0 0 0 0 0 1 0 0 0 0 7
8 16670.7 0 0 0 0 0 0 0 0 1 0 0 0 8
9 16614.8 0 0 0 0 0 0 0 0 0 1 0 0 9
10 16869.2 0 0 0 0 0 0 0 0 0 0 1 0 10
11 15663.9 0 0 0 0 0 0 0 0 0 0 0 1 11
12 16359.9 0 0 0 0 0 0 0 0 0 0 0 0 12
13 18447.7 0 1 0 0 0 0 0 0 0 0 0 0 13
14 16889.0 0 0 1 0 0 0 0 0 0 0 0 0 14
15 16505.0 0 0 0 1 0 0 0 0 0 0 0 0 15
16 18320.9 0 0 0 0 1 0 0 0 0 0 0 0 16
17 15052.1 0 0 0 0 0 1 0 0 0 0 0 0 17
18 15699.8 0 0 0 0 0 0 1 0 0 0 0 0 18
19 18135.3 0 0 0 0 0 0 0 1 0 0 0 0 19
20 16768.7 0 0 0 0 0 0 0 0 1 0 0 0 20
21 18883.0 0 0 0 0 0 0 0 0 0 1 0 0 21
22 19021.0 0 0 0 0 0 0 0 0 0 0 1 0 22
23 18101.9 0 0 0 0 0 0 0 0 0 0 0 1 23
24 17776.1 0 0 0 0 0 0 0 0 0 0 0 0 24
25 21489.9 0 1 0 0 0 0 0 0 0 0 0 0 25
26 17065.3 0 0 1 0 0 0 0 0 0 0 0 0 26
27 18690.0 0 0 0 1 0 0 0 0 0 0 0 0 27
28 18953.1 0 0 0 0 1 0 0 0 0 0 0 0 28
29 16398.9 0 0 0 0 0 1 0 0 0 0 0 0 29
30 16895.6 0 0 0 0 0 0 1 0 0 0 0 0 30
31 18553.0 0 0 0 0 0 0 0 1 0 0 0 0 31
32 19270.0 0 0 0 0 0 0 0 0 1 0 0 0 32
33 19422.1 0 0 0 0 0 0 0 0 0 1 0 0 33
34 17579.4 0 0 0 0 0 0 0 0 0 0 1 0 34
35 18637.3 0 0 0 0 0 0 0 0 0 0 0 1 35
36 18076.7 0 0 0 0 0 0 0 0 0 0 0 0 36
37 20438.6 0 1 0 0 0 0 0 0 0 0 0 0 37
38 18075.2 0 0 1 0 0 0 0 0 0 0 0 0 38
39 19563.0 0 0 0 1 0 0 0 0 0 0 0 0 39
40 19899.2 0 0 0 0 1 0 0 0 0 0 0 0 40
41 19227.5 0 0 0 0 0 1 0 0 0 0 0 0 41
42 17789.6 0 0 0 0 0 0 1 0 0 0 0 0 42
43 19220.8 0 0 0 0 0 0 0 1 0 0 0 0 43
44 21968.9 0 0 0 0 0 0 0 0 1 0 0 0 44
45 21131.5 0 0 0 0 0 0 0 0 0 1 0 0 45
46 19484.6 0 0 0 0 0 0 0 0 0 0 1 0 46
47 22168.7 1 0 0 0 0 0 0 0 0 0 0 1 47
48 20866.8 1 0 0 0 0 0 0 0 0 0 0 0 48
49 22176.2 1 1 0 0 0 0 0 0 0 0 0 0 49
50 23533.8 1 0 1 0 0 0 0 0 0 0 0 0 50
51 21479.6 1 0 0 1 0 0 0 0 0 0 0 0 51
52 24347.7 1 0 0 0 1 0 0 0 0 0 0 0 52
53 22751.6 1 0 0 0 0 1 0 0 0 0 0 0 53
54 20328.3 1 0 0 0 0 0 1 0 0 0 0 0 54
55 23650.4 1 0 0 0 0 0 0 1 0 0 0 0 55
56 23335.7 1 0 0 0 0 0 0 0 1 0 0 0 56
57 19614.9 1 0 0 0 0 0 0 0 0 1 0 0 57
58 18042.3 1 0 0 0 0 0 0 0 0 0 1 0 58
59 17282.5 1 0 0 0 0 0 0 0 0 0 0 1 59
60 16847.2 1 0 0 0 0 0 0 0 0 0 0 0 60
61 18159.5 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) Dummy M1 M2 M3 M4
14662.87 649.50 2134.05 1185.18 1015.17 2634.65
M5 M6 M7 M8 M9 M10
325.50 -322.67 1807.75 2087.66 1533.04 514.01
M11 t
470.59 85.07
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4476.441 -726.499 3.181 813.279 2782.539
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14662.87 874.70 16.763 < 2e-16 ***
Dummy 649.50 749.61 0.866 0.3906
M1 2134.05 994.86 2.145 0.0371 *
M2 1185.18 1044.11 1.135 0.2621
M3 1015.17 1043.03 0.973 0.3354
M4 2634.65 1042.28 2.528 0.0149 *
M5 325.50 1041.84 0.312 0.7561
M6 -322.67 1041.73 -0.310 0.7581
M7 1807.75 1041.94 1.735 0.0893 .
M8 2087.66 1042.47 2.003 0.0510 .
M9 1533.04 1043.32 1.469 0.1484
M10 514.01 1044.49 0.492 0.6249
M11 470.59 1037.16 0.454 0.6521
t 85.07 18.29 4.651 2.71e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1640 on 47 degrees of freedom
Multiple R-squared: 0.6348, Adjusted R-squared: 0.5338
F-statistic: 6.285 on 13 and 47 DF, p-value: 1.198e-06
> 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.0301077646 0.0602155291 0.96989224
[2,] 0.0069645088 0.0139290176 0.99303549
[3,] 0.0015491187 0.0030982373 0.99845088
[4,] 0.0020441542 0.0040883084 0.99795585
[5,] 0.0013803663 0.0027607326 0.99861963
[6,] 0.0006433189 0.0012866378 0.99935668
[7,] 0.0003741890 0.0007483780 0.99962581
[8,] 0.0001003505 0.0002007010 0.99989965
[9,] 0.0001817240 0.0003634481 0.99981828
[10,] 0.0003054504 0.0006109008 0.99969455
[11,] 0.0001864693 0.0003729387 0.99981353
[12,] 0.0001998510 0.0003997020 0.99980015
[13,] 0.0002194436 0.0004388872 0.99978056
[14,] 0.0001071262 0.0002142523 0.99989287
[15,] 0.0001298131 0.0002596262 0.99987019
[16,] 0.0003303100 0.0006606200 0.99966969
[17,] 0.0003013794 0.0006027587 0.99969862
[18,] 0.0067652710 0.0135305420 0.99323473
[19,] 0.0035827525 0.0071655051 0.99641725
[20,] 0.0021816968 0.0043633937 0.99781830
[21,] 0.0013686672 0.0027373345 0.99863133
[22,] 0.0032423520 0.0064847040 0.99675765
[23,] 0.0017455510 0.0034911019 0.99825445
[24,] 0.0030805692 0.0061611384 0.99691943
[25,] 0.0070154413 0.0140308827 0.99298456
[26,] 0.0052396224 0.0104792449 0.99476038
[27,] 0.2966072476 0.5932144952 0.70339275
[28,] 0.9425593466 0.1148813068 0.05744065
> postscript(file="/var/www/html/rcomp/tmp/17a4l1260617728.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/22n681260617728.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/36q311260617728.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/4updx1260617728.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/5sqwa1260617728.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
310.407295 -632.101000 -1646.161000 -111.221000 -916.741000 -452.341000
7 8 9 10 11 12
-436.541000 -760.421000 -346.781000 841.579000 -405.380492 676.139508
13 14 15 16 17 18
544.818049 -150.090246 -449.150246 -337.810246 -1382.530246 -171.730246
19 20 21 22 23 24
48.269754 -1683.310246 900.529754 1972.489754 1011.730262 1071.450262
25 26 27 28 29 30
2566.128803 -994.679492 714.960508 -726.499492 -1056.619492 3.180508
31 32 33 34 35 36
-554.919492 -202.899492 418.740508 -489.999492 526.241016 351.161016
37 38 39 40 41 42
493.939557 -1005.668738 567.071262 -801.288738 751.091262 -123.708738
43 44 45 46 47 48
-908.008738 1475.111262 1107.251262 394.311262 2387.249230 1470.869230
49 50 51 52 53 54
561.147770 2782.539475 813.279475 1976.819475 2604.799475 744.599475
55 56 57 58 59 60
1851.199475 1171.519475 -2079.740525 -2718.380525 -3519.840016 -3569.620016
61
-4476.441475
> postscript(file="/var/www/html/rcomp/tmp/6m8ts1260617728.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 310.407295 NA
1 -632.101000 310.407295
2 -1646.161000 -632.101000
3 -111.221000 -1646.161000
4 -916.741000 -111.221000
5 -452.341000 -916.741000
6 -436.541000 -452.341000
7 -760.421000 -436.541000
8 -346.781000 -760.421000
9 841.579000 -346.781000
10 -405.380492 841.579000
11 676.139508 -405.380492
12 544.818049 676.139508
13 -150.090246 544.818049
14 -449.150246 -150.090246
15 -337.810246 -449.150246
16 -1382.530246 -337.810246
17 -171.730246 -1382.530246
18 48.269754 -171.730246
19 -1683.310246 48.269754
20 900.529754 -1683.310246
21 1972.489754 900.529754
22 1011.730262 1972.489754
23 1071.450262 1011.730262
24 2566.128803 1071.450262
25 -994.679492 2566.128803
26 714.960508 -994.679492
27 -726.499492 714.960508
28 -1056.619492 -726.499492
29 3.180508 -1056.619492
30 -554.919492 3.180508
31 -202.899492 -554.919492
32 418.740508 -202.899492
33 -489.999492 418.740508
34 526.241016 -489.999492
35 351.161016 526.241016
36 493.939557 351.161016
37 -1005.668738 493.939557
38 567.071262 -1005.668738
39 -801.288738 567.071262
40 751.091262 -801.288738
41 -123.708738 751.091262
42 -908.008738 -123.708738
43 1475.111262 -908.008738
44 1107.251262 1475.111262
45 394.311262 1107.251262
46 2387.249230 394.311262
47 1470.869230 2387.249230
48 561.147770 1470.869230
49 2782.539475 561.147770
50 813.279475 2782.539475
51 1976.819475 813.279475
52 2604.799475 1976.819475
53 744.599475 2604.799475
54 1851.199475 744.599475
55 1171.519475 1851.199475
56 -2079.740525 1171.519475
57 -2718.380525 -2079.740525
58 -3519.840016 -2718.380525
59 -3569.620016 -3519.840016
60 -4476.441475 -3569.620016
61 NA -4476.441475
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -632.101000 310.407295
[2,] -1646.161000 -632.101000
[3,] -111.221000 -1646.161000
[4,] -916.741000 -111.221000
[5,] -452.341000 -916.741000
[6,] -436.541000 -452.341000
[7,] -760.421000 -436.541000
[8,] -346.781000 -760.421000
[9,] 841.579000 -346.781000
[10,] -405.380492 841.579000
[11,] 676.139508 -405.380492
[12,] 544.818049 676.139508
[13,] -150.090246 544.818049
[14,] -449.150246 -150.090246
[15,] -337.810246 -449.150246
[16,] -1382.530246 -337.810246
[17,] -171.730246 -1382.530246
[18,] 48.269754 -171.730246
[19,] -1683.310246 48.269754
[20,] 900.529754 -1683.310246
[21,] 1972.489754 900.529754
[22,] 1011.730262 1972.489754
[23,] 1071.450262 1011.730262
[24,] 2566.128803 1071.450262
[25,] -994.679492 2566.128803
[26,] 714.960508 -994.679492
[27,] -726.499492 714.960508
[28,] -1056.619492 -726.499492
[29,] 3.180508 -1056.619492
[30,] -554.919492 3.180508
[31,] -202.899492 -554.919492
[32,] 418.740508 -202.899492
[33,] -489.999492 418.740508
[34,] 526.241016 -489.999492
[35,] 351.161016 526.241016
[36,] 493.939557 351.161016
[37,] -1005.668738 493.939557
[38,] 567.071262 -1005.668738
[39,] -801.288738 567.071262
[40,] 751.091262 -801.288738
[41,] -123.708738 751.091262
[42,] -908.008738 -123.708738
[43,] 1475.111262 -908.008738
[44,] 1107.251262 1475.111262
[45,] 394.311262 1107.251262
[46,] 2387.249230 394.311262
[47,] 1470.869230 2387.249230
[48,] 561.147770 1470.869230
[49,] 2782.539475 561.147770
[50,] 813.279475 2782.539475
[51,] 1976.819475 813.279475
[52,] 2604.799475 1976.819475
[53,] 744.599475 2604.799475
[54,] 1851.199475 744.599475
[55,] 1171.519475 1851.199475
[56,] -2079.740525 1171.519475
[57,] -2718.380525 -2079.740525
[58,] -3519.840016 -2718.380525
[59,] -3569.620016 -3519.840016
[60,] -4476.441475 -3569.620016
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -632.101000 310.407295
2 -1646.161000 -632.101000
3 -111.221000 -1646.161000
4 -916.741000 -111.221000
5 -452.341000 -916.741000
6 -436.541000 -452.341000
7 -760.421000 -436.541000
8 -346.781000 -760.421000
9 841.579000 -346.781000
10 -405.380492 841.579000
11 676.139508 -405.380492
12 544.818049 676.139508
13 -150.090246 544.818049
14 -449.150246 -150.090246
15 -337.810246 -449.150246
16 -1382.530246 -337.810246
17 -171.730246 -1382.530246
18 48.269754 -171.730246
19 -1683.310246 48.269754
20 900.529754 -1683.310246
21 1972.489754 900.529754
22 1011.730262 1972.489754
23 1071.450262 1011.730262
24 2566.128803 1071.450262
25 -994.679492 2566.128803
26 714.960508 -994.679492
27 -726.499492 714.960508
28 -1056.619492 -726.499492
29 3.180508 -1056.619492
30 -554.919492 3.180508
31 -202.899492 -554.919492
32 418.740508 -202.899492
33 -489.999492 418.740508
34 526.241016 -489.999492
35 351.161016 526.241016
36 493.939557 351.161016
37 -1005.668738 493.939557
38 567.071262 -1005.668738
39 -801.288738 567.071262
40 751.091262 -801.288738
41 -123.708738 751.091262
42 -908.008738 -123.708738
43 1475.111262 -908.008738
44 1107.251262 1475.111262
45 394.311262 1107.251262
46 2387.249230 394.311262
47 1470.869230 2387.249230
48 561.147770 1470.869230
49 2782.539475 561.147770
50 813.279475 2782.539475
51 1976.819475 813.279475
52 2604.799475 1976.819475
53 744.599475 2604.799475
54 1851.199475 744.599475
55 1171.519475 1851.199475
56 -2079.740525 1171.519475
57 -2718.380525 -2079.740525
58 -3519.840016 -2718.380525
59 -3569.620016 -3519.840016
60 -4476.441475 -3569.620016
> 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/73kky1260617728.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/8gej01260617728.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/9xcvg1260617728.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/105hu71260617728.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/11fktb1260617728.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/12ir831260617728.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/13e17r1260617728.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/14m1dc1260617728.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/15xj0y1260617728.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/167gp41260617728.tab")
+ }
>
> system("convert tmp/17a4l1260617728.ps tmp/17a4l1260617728.png")
> system("convert tmp/22n681260617728.ps tmp/22n681260617728.png")
> system("convert tmp/36q311260617728.ps tmp/36q311260617728.png")
> system("convert tmp/4updx1260617728.ps tmp/4updx1260617728.png")
> system("convert tmp/5sqwa1260617728.ps tmp/5sqwa1260617728.png")
> system("convert tmp/6m8ts1260617728.ps tmp/6m8ts1260617728.png")
> system("convert tmp/73kky1260617728.ps tmp/73kky1260617728.png")
> system("convert tmp/8gej01260617728.ps tmp/8gej01260617728.png")
> system("convert tmp/9xcvg1260617728.ps tmp/9xcvg1260617728.png")
> system("convert tmp/105hu71260617728.ps tmp/105hu71260617728.png")
>
>
> proc.time()
user system elapsed
2.374 1.549 4.149