R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(111078,0,150739,0,159129,0,157928,0,147768,0,137507,0,136919,0,136151,0,133001,0,125554,0,119647,0,114158,0,116193,0,152803,0,161761,0,160942,0,149470,0,139208,0,134588,0,130322,0,126611,0,122401,0,117352,0,112135,0,112879,0,148729,0,157230,0,157221,0,146681,0,136524,0,132111,1,125326,1,122716,1,116615,1,113719,1,110737,1,112093,1,143565,1,149946,1,149147,1,134339,1,122683,1,115614,1,116566,1,111272,1,104609,1,101802,1,94542,1,93051,1,124129,1,130374,1,123946,1,114971,1,105531,1,104919,1,104782,1,101281,1,94545,1,93248,1,84031,1,87486,1),dim=c(2,61),dimnames=list(c('Yt','D'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Yt','D'),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
Yt D M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 111078 0 1 0 0 0 0 0 0 0 0 0 0 1
2 150739 0 0 1 0 0 0 0 0 0 0 0 0 2
3 159129 0 0 0 1 0 0 0 0 0 0 0 0 3
4 157928 0 0 0 0 1 0 0 0 0 0 0 0 4
5 147768 0 0 0 0 0 1 0 0 0 0 0 0 5
6 137507 0 0 0 0 0 0 1 0 0 0 0 0 6
7 136919 0 0 0 0 0 0 0 1 0 0 0 0 7
8 136151 0 0 0 0 0 0 0 0 1 0 0 0 8
9 133001 0 0 0 0 0 0 0 0 0 1 0 0 9
10 125554 0 0 0 0 0 0 0 0 0 0 1 0 10
11 119647 0 0 0 0 0 0 0 0 0 0 0 1 11
12 114158 0 0 0 0 0 0 0 0 0 0 0 0 12
13 116193 0 1 0 0 0 0 0 0 0 0 0 0 13
14 152803 0 0 1 0 0 0 0 0 0 0 0 0 14
15 161761 0 0 0 1 0 0 0 0 0 0 0 0 15
16 160942 0 0 0 0 1 0 0 0 0 0 0 0 16
17 149470 0 0 0 0 0 1 0 0 0 0 0 0 17
18 139208 0 0 0 0 0 0 1 0 0 0 0 0 18
19 134588 0 0 0 0 0 0 0 1 0 0 0 0 19
20 130322 0 0 0 0 0 0 0 0 1 0 0 0 20
21 126611 0 0 0 0 0 0 0 0 0 1 0 0 21
22 122401 0 0 0 0 0 0 0 0 0 0 1 0 22
23 117352 0 0 0 0 0 0 0 0 0 0 0 1 23
24 112135 0 0 0 0 0 0 0 0 0 0 0 0 24
25 112879 0 1 0 0 0 0 0 0 0 0 0 0 25
26 148729 0 0 1 0 0 0 0 0 0 0 0 0 26
27 157230 0 0 0 1 0 0 0 0 0 0 0 0 27
28 157221 0 0 0 0 1 0 0 0 0 0 0 0 28
29 146681 0 0 0 0 0 1 0 0 0 0 0 0 29
30 136524 0 0 0 0 0 0 1 0 0 0 0 0 30
31 132111 1 0 0 0 0 0 0 1 0 0 0 0 31
32 125326 1 0 0 0 0 0 0 0 1 0 0 0 32
33 122716 1 0 0 0 0 0 0 0 0 1 0 0 33
34 116615 1 0 0 0 0 0 0 0 0 0 1 0 34
35 113719 1 0 0 0 0 0 0 0 0 0 0 1 35
36 110737 1 0 0 0 0 0 0 0 0 0 0 0 36
37 112093 1 1 0 0 0 0 0 0 0 0 0 0 37
38 143565 1 0 1 0 0 0 0 0 0 0 0 0 38
39 149946 1 0 0 1 0 0 0 0 0 0 0 0 39
40 149147 1 0 0 0 1 0 0 0 0 0 0 0 40
41 134339 1 0 0 0 0 1 0 0 0 0 0 0 41
42 122683 1 0 0 0 0 0 1 0 0 0 0 0 42
43 115614 1 0 0 0 0 0 0 1 0 0 0 0 43
44 116566 1 0 0 0 0 0 0 0 1 0 0 0 44
45 111272 1 0 0 0 0 0 0 0 0 1 0 0 45
46 104609 1 0 0 0 0 0 0 0 0 0 1 0 46
47 101802 1 0 0 0 0 0 0 0 0 0 0 1 47
48 94542 1 0 0 0 0 0 0 0 0 0 0 0 48
49 93051 1 1 0 0 0 0 0 0 0 0 0 0 49
50 124129 1 0 1 0 0 0 0 0 0 0 0 0 50
51 130374 1 0 0 1 0 0 0 0 0 0 0 0 51
52 123946 1 0 0 0 1 0 0 0 0 0 0 0 52
53 114971 1 0 0 0 0 1 0 0 0 0 0 0 53
54 105531 1 0 0 0 0 0 1 0 0 0 0 0 54
55 104919 1 0 0 0 0 0 0 1 0 0 0 0 55
56 104782 1 0 0 0 0 0 0 0 1 0 0 0 56
57 101281 1 0 0 0 0 0 0 0 0 1 0 0 57
58 94545 1 0 0 0 0 0 0 0 0 0 1 0 58
59 93248 1 0 0 0 0 0 0 0 0 0 0 1 59
60 84031 1 0 0 0 0 0 0 0 0 0 0 0 60
61 87486 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) D M1 M2 M3 M4
124865.4 -577.3 -687.0 34812.9 43102.3 41845.5
M5 M6 M7 M8 M9 M10
31248.9 21488.1 18737.6 17131.2 14072.4 8435.4
M11 t
5438.6 -594.4
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12506.0 -3574.6 143.3 3560.6 10484.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 124865.40 3137.47 39.798 < 2e-16 ***
D -577.33 3028.95 -0.191 0.849656
M1 -687.00 3521.13 -0.195 0.846149
M2 34812.93 3697.62 9.415 2.16e-12 ***
M3 43102.33 3690.91 11.678 1.70e-15 ***
M4 41845.53 3686.17 11.352 4.57e-15 ***
M5 31248.93 3683.43 8.484 4.88e-11 ***
M6 21488.13 3682.67 5.835 4.78e-07 ***
M7 18737.60 3695.13 5.071 6.62e-06 ***
M8 17131.20 3686.17 4.647 2.74e-05 ***
M9 14072.40 3679.19 3.825 0.000385 ***
M10 8435.40 3674.20 2.296 0.026191 *
M11 5438.60 3671.20 1.481 0.145167
t -594.40 85.71 -6.935 1.02e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5803 on 47 degrees of freedom
Multiple R-squared: 0.9359, Adjusted R-squared: 0.9182
F-statistic: 52.82 on 13 and 47 DF, p-value: < 2.2e-16
> 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.025598767 0.051197534 0.97440123
[2,] 0.007570249 0.015140498 0.99242975
[3,] 0.030481153 0.060962306 0.96951885
[4,] 0.133618600 0.267237199 0.86638140
[5,] 0.222532222 0.445064445 0.77746778
[6,] 0.174248974 0.348497948 0.82575103
[7,] 0.149942254 0.299884509 0.85005775
[8,] 0.132368809 0.264737618 0.86763119
[9,] 0.153141363 0.306282726 0.84685864
[10,] 0.132400108 0.264800216 0.86759989
[11,] 0.105797548 0.211595096 0.89420245
[12,] 0.067040863 0.134081726 0.93295914
[13,] 0.039380946 0.078761891 0.96061905
[14,] 0.021841278 0.043682556 0.97815872
[15,] 0.012504199 0.025008399 0.98749580
[16,] 0.012672356 0.025344712 0.98732764
[17,] 0.008056033 0.016112066 0.99194397
[18,] 0.004890050 0.009780101 0.99510995
[19,] 0.003845950 0.007691900 0.99615405
[20,] 0.002904308 0.005808617 0.99709569
[21,] 0.003020171 0.006040341 0.99697983
[22,] 0.002601774 0.005203548 0.99739823
[23,] 0.003927576 0.007855152 0.99607242
[24,] 0.090761209 0.181522418 0.90923879
[25,] 0.444281806 0.888563612 0.55571819
[26,] 0.892354293 0.215291414 0.10764571
[27,] 0.902504153 0.194991694 0.09749585
[28,] 0.876740365 0.246519270 0.12325963
> postscript(file="/var/www/html/rcomp/tmp/1lfe61229626768.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/2l9c21229626768.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/373ry1229626768.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/4j4pv1229626768.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/5r7mi1229626768.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
-12506.00000 -7750.53333 -7055.53333 -6405.33333 -5374.33333 -5280.13333
7 8 9 10 11 12
-2523.20000 -1090.40000 -587.20000 -1802.80000 -4118.60000 -3574.60000
13 14 15 16 17 18
-258.20000 1446.26667 2709.26667 3741.46667 3460.46667 3553.66667
19 20 21 22 23 24
2278.60000 213.40000 155.60000 2177.00000 719.20000 1535.20000
25 26 27 28 29 30
3560.60000 4505.06667 5311.06667 7153.26667 7804.26667 8002.46667
31 32 33 34 35 36
7511.73333 2927.53333 3970.73333 4101.13333 4796.33333 7847.33333
37 38 39 40 41 42
10484.73333 7051.20000 5737.20000 6789.40000 3172.40000 1871.60000
43 44 45 46 47 48
-1852.46667 1300.33333 -340.46667 -772.06667 12.13333 -1214.86667
49 50 51 52 53 54
-1424.46667 -5252.00000 -6702.00000 -11278.80000 -9062.80000 -8147.60000
55 56 57 58 59 60
-5414.66667 -3350.86667 -3198.66667 -3703.26667 -1409.06667 -4593.06667
61
143.33333
> postscript(file="/var/www/html/rcomp/tmp/6ha0m1229626768.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 -12506.00000 NA
1 -7750.53333 -12506.00000
2 -7055.53333 -7750.53333
3 -6405.33333 -7055.53333
4 -5374.33333 -6405.33333
5 -5280.13333 -5374.33333
6 -2523.20000 -5280.13333
7 -1090.40000 -2523.20000
8 -587.20000 -1090.40000
9 -1802.80000 -587.20000
10 -4118.60000 -1802.80000
11 -3574.60000 -4118.60000
12 -258.20000 -3574.60000
13 1446.26667 -258.20000
14 2709.26667 1446.26667
15 3741.46667 2709.26667
16 3460.46667 3741.46667
17 3553.66667 3460.46667
18 2278.60000 3553.66667
19 213.40000 2278.60000
20 155.60000 213.40000
21 2177.00000 155.60000
22 719.20000 2177.00000
23 1535.20000 719.20000
24 3560.60000 1535.20000
25 4505.06667 3560.60000
26 5311.06667 4505.06667
27 7153.26667 5311.06667
28 7804.26667 7153.26667
29 8002.46667 7804.26667
30 7511.73333 8002.46667
31 2927.53333 7511.73333
32 3970.73333 2927.53333
33 4101.13333 3970.73333
34 4796.33333 4101.13333
35 7847.33333 4796.33333
36 10484.73333 7847.33333
37 7051.20000 10484.73333
38 5737.20000 7051.20000
39 6789.40000 5737.20000
40 3172.40000 6789.40000
41 1871.60000 3172.40000
42 -1852.46667 1871.60000
43 1300.33333 -1852.46667
44 -340.46667 1300.33333
45 -772.06667 -340.46667
46 12.13333 -772.06667
47 -1214.86667 12.13333
48 -1424.46667 -1214.86667
49 -5252.00000 -1424.46667
50 -6702.00000 -5252.00000
51 -11278.80000 -6702.00000
52 -9062.80000 -11278.80000
53 -8147.60000 -9062.80000
54 -5414.66667 -8147.60000
55 -3350.86667 -5414.66667
56 -3198.66667 -3350.86667
57 -3703.26667 -3198.66667
58 -1409.06667 -3703.26667
59 -4593.06667 -1409.06667
60 143.33333 -4593.06667
61 NA 143.33333
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -7750.53333 -12506.00000
[2,] -7055.53333 -7750.53333
[3,] -6405.33333 -7055.53333
[4,] -5374.33333 -6405.33333
[5,] -5280.13333 -5374.33333
[6,] -2523.20000 -5280.13333
[7,] -1090.40000 -2523.20000
[8,] -587.20000 -1090.40000
[9,] -1802.80000 -587.20000
[10,] -4118.60000 -1802.80000
[11,] -3574.60000 -4118.60000
[12,] -258.20000 -3574.60000
[13,] 1446.26667 -258.20000
[14,] 2709.26667 1446.26667
[15,] 3741.46667 2709.26667
[16,] 3460.46667 3741.46667
[17,] 3553.66667 3460.46667
[18,] 2278.60000 3553.66667
[19,] 213.40000 2278.60000
[20,] 155.60000 213.40000
[21,] 2177.00000 155.60000
[22,] 719.20000 2177.00000
[23,] 1535.20000 719.20000
[24,] 3560.60000 1535.20000
[25,] 4505.06667 3560.60000
[26,] 5311.06667 4505.06667
[27,] 7153.26667 5311.06667
[28,] 7804.26667 7153.26667
[29,] 8002.46667 7804.26667
[30,] 7511.73333 8002.46667
[31,] 2927.53333 7511.73333
[32,] 3970.73333 2927.53333
[33,] 4101.13333 3970.73333
[34,] 4796.33333 4101.13333
[35,] 7847.33333 4796.33333
[36,] 10484.73333 7847.33333
[37,] 7051.20000 10484.73333
[38,] 5737.20000 7051.20000
[39,] 6789.40000 5737.20000
[40,] 3172.40000 6789.40000
[41,] 1871.60000 3172.40000
[42,] -1852.46667 1871.60000
[43,] 1300.33333 -1852.46667
[44,] -340.46667 1300.33333
[45,] -772.06667 -340.46667
[46,] 12.13333 -772.06667
[47,] -1214.86667 12.13333
[48,] -1424.46667 -1214.86667
[49,] -5252.00000 -1424.46667
[50,] -6702.00000 -5252.00000
[51,] -11278.80000 -6702.00000
[52,] -9062.80000 -11278.80000
[53,] -8147.60000 -9062.80000
[54,] -5414.66667 -8147.60000
[55,] -3350.86667 -5414.66667
[56,] -3198.66667 -3350.86667
[57,] -3703.26667 -3198.66667
[58,] -1409.06667 -3703.26667
[59,] -4593.06667 -1409.06667
[60,] 143.33333 -4593.06667
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -7750.53333 -12506.00000
2 -7055.53333 -7750.53333
3 -6405.33333 -7055.53333
4 -5374.33333 -6405.33333
5 -5280.13333 -5374.33333
6 -2523.20000 -5280.13333
7 -1090.40000 -2523.20000
8 -587.20000 -1090.40000
9 -1802.80000 -587.20000
10 -4118.60000 -1802.80000
11 -3574.60000 -4118.60000
12 -258.20000 -3574.60000
13 1446.26667 -258.20000
14 2709.26667 1446.26667
15 3741.46667 2709.26667
16 3460.46667 3741.46667
17 3553.66667 3460.46667
18 2278.60000 3553.66667
19 213.40000 2278.60000
20 155.60000 213.40000
21 2177.00000 155.60000
22 719.20000 2177.00000
23 1535.20000 719.20000
24 3560.60000 1535.20000
25 4505.06667 3560.60000
26 5311.06667 4505.06667
27 7153.26667 5311.06667
28 7804.26667 7153.26667
29 8002.46667 7804.26667
30 7511.73333 8002.46667
31 2927.53333 7511.73333
32 3970.73333 2927.53333
33 4101.13333 3970.73333
34 4796.33333 4101.13333
35 7847.33333 4796.33333
36 10484.73333 7847.33333
37 7051.20000 10484.73333
38 5737.20000 7051.20000
39 6789.40000 5737.20000
40 3172.40000 6789.40000
41 1871.60000 3172.40000
42 -1852.46667 1871.60000
43 1300.33333 -1852.46667
44 -340.46667 1300.33333
45 -772.06667 -340.46667
46 12.13333 -772.06667
47 -1214.86667 12.13333
48 -1424.46667 -1214.86667
49 -5252.00000 -1424.46667
50 -6702.00000 -5252.00000
51 -11278.80000 -6702.00000
52 -9062.80000 -11278.80000
53 -8147.60000 -9062.80000
54 -5414.66667 -8147.60000
55 -3350.86667 -5414.66667
56 -3198.66667 -3350.86667
57 -3703.26667 -3198.66667
58 -1409.06667 -3703.26667
59 -4593.06667 -1409.06667
60 143.33333 -4593.06667
> 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/79vjn1229626768.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/8y60y1229626768.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/937vy1229626768.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/10vi0k1229626768.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/11xhjg1229626768.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/12h21n1229626768.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/13rh5x1229626769.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/145jrv1229626769.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/15g4er1229626769.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/16ti3g1229626769.tab")
+ }
>
> system("convert tmp/1lfe61229626768.ps tmp/1lfe61229626768.png")
> system("convert tmp/2l9c21229626768.ps tmp/2l9c21229626768.png")
> system("convert tmp/373ry1229626768.ps tmp/373ry1229626768.png")
> system("convert tmp/4j4pv1229626768.ps tmp/4j4pv1229626768.png")
> system("convert tmp/5r7mi1229626768.ps tmp/5r7mi1229626768.png")
> system("convert tmp/6ha0m1229626768.ps tmp/6ha0m1229626768.png")
> system("convert tmp/79vjn1229626768.ps tmp/79vjn1229626768.png")
> system("convert tmp/8y60y1229626768.ps tmp/8y60y1229626768.png")
> system("convert tmp/937vy1229626768.ps tmp/937vy1229626768.png")
> system("convert tmp/10vi0k1229626768.ps tmp/10vi0k1229626768.png")
>
>
> proc.time()
user system elapsed
2.447 1.579 3.074