R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(4143,4429,5219,4929,5761,5592,4163,4962,5208,4755,4491,5732,5731,5040,6102,4904,5369,5578,4619,4731,5011,5299,4146,4625,4736,4219,5116,4205,4121,5103,4300,4578,3809,5657,4248,3830,4736,4839,4411,4570,4104,4801,3953,3828,4440,4026,4109,4785,3224,3552,3940,3913,3681,4309,3830,4143,4087,3818,3380,3430,3458,3970,5260,5024,5634,6549,4676),dim=c(1,67),dimnames=list(c('nb'),1:67))
> y <- array(NA,dim=c(1,67),dimnames=list(c('nb'),1:67))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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
nb M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 4143 1 0 0 0 0 0 0 0 0 0 0
2 4429 0 1 0 0 0 0 0 0 0 0 0
3 5219 0 0 1 0 0 0 0 0 0 0 0
4 4929 0 0 0 1 0 0 0 0 0 0 0
5 5761 0 0 0 0 1 0 0 0 0 0 0
6 5592 0 0 0 0 0 1 0 0 0 0 0
7 4163 0 0 0 0 0 0 1 0 0 0 0
8 4962 0 0 0 0 0 0 0 1 0 0 0
9 5208 0 0 0 0 0 0 0 0 1 0 0
10 4755 0 0 0 0 0 0 0 0 0 1 0
11 4491 0 0 0 0 0 0 0 0 0 0 1
12 5732 0 0 0 0 0 0 0 0 0 0 0
13 5731 1 0 0 0 0 0 0 0 0 0 0
14 5040 0 1 0 0 0 0 0 0 0 0 0
15 6102 0 0 1 0 0 0 0 0 0 0 0
16 4904 0 0 0 1 0 0 0 0 0 0 0
17 5369 0 0 0 0 1 0 0 0 0 0 0
18 5578 0 0 0 0 0 1 0 0 0 0 0
19 4619 0 0 0 0 0 0 1 0 0 0 0
20 4731 0 0 0 0 0 0 0 1 0 0 0
21 5011 0 0 0 0 0 0 0 0 1 0 0
22 5299 0 0 0 0 0 0 0 0 0 1 0
23 4146 0 0 0 0 0 0 0 0 0 0 1
24 4625 0 0 0 0 0 0 0 0 0 0 0
25 4736 1 0 0 0 0 0 0 0 0 0 0
26 4219 0 1 0 0 0 0 0 0 0 0 0
27 5116 0 0 1 0 0 0 0 0 0 0 0
28 4205 0 0 0 1 0 0 0 0 0 0 0
29 4121 0 0 0 0 1 0 0 0 0 0 0
30 5103 0 0 0 0 0 1 0 0 0 0 0
31 4300 0 0 0 0 0 0 1 0 0 0 0
32 4578 0 0 0 0 0 0 0 1 0 0 0
33 3809 0 0 0 0 0 0 0 0 1 0 0
34 5657 0 0 0 0 0 0 0 0 0 1 0
35 4248 0 0 0 0 0 0 0 0 0 0 1
36 3830 0 0 0 0 0 0 0 0 0 0 0
37 4736 1 0 0 0 0 0 0 0 0 0 0
38 4839 0 1 0 0 0 0 0 0 0 0 0
39 4411 0 0 1 0 0 0 0 0 0 0 0
40 4570 0 0 0 1 0 0 0 0 0 0 0
41 4104 0 0 0 0 1 0 0 0 0 0 0
42 4801 0 0 0 0 0 1 0 0 0 0 0
43 3953 0 0 0 0 0 0 1 0 0 0 0
44 3828 0 0 0 0 0 0 0 1 0 0 0
45 4440 0 0 0 0 0 0 0 0 1 0 0
46 4026 0 0 0 0 0 0 0 0 0 1 0
47 4109 0 0 0 0 0 0 0 0 0 0 1
48 4785 0 0 0 0 0 0 0 0 0 0 0
49 3224 1 0 0 0 0 0 0 0 0 0 0
50 3552 0 1 0 0 0 0 0 0 0 0 0
51 3940 0 0 1 0 0 0 0 0 0 0 0
52 3913 0 0 0 1 0 0 0 0 0 0 0
53 3681 0 0 0 0 1 0 0 0 0 0 0
54 4309 0 0 0 0 0 1 0 0 0 0 0
55 3830 0 0 0 0 0 0 1 0 0 0 0
56 4143 0 0 0 0 0 0 0 1 0 0 0
57 4087 0 0 0 0 0 0 0 0 1 0 0
58 3818 0 0 0 0 0 0 0 0 0 1 0
59 3380 0 0 0 0 0 0 0 0 0 0 1
60 3430 0 0 0 0 0 0 0 0 0 0 0
61 3458 1 0 0 0 0 0 0 0 0 0 0
62 3970 0 1 0 0 0 0 0 0 0 0 0
63 5260 0 0 1 0 0 0 0 0 0 0 0
64 5024 0 0 0 1 0 0 0 0 0 0 0
65 5634 0 0 0 0 1 0 0 0 0 0 0
66 6549 0 0 0 0 0 1 0 0 0 0 0
67 4676 0 0 0 0 0 0 1 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
4480.4 -142.4 -138.9 527.6 110.4 297.9
M6 M7 M8 M9 M10 M11
841.6 -223.6 -32.0 30.6 230.6 -405.6
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1114.0 -559.0 71.2 407.1 1393.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4480.4 307.7 14.562 <2e-16 ***
M1 -142.4 416.6 -0.342 0.7338
M2 -138.9 416.6 -0.333 0.7401
M3 527.6 416.6 1.266 0.2107
M4 110.4 416.6 0.265 0.7919
M5 297.9 416.6 0.715 0.4775
M6 841.6 416.6 2.020 0.0482 *
M7 -223.6 416.6 -0.537 0.5937
M8 -32.0 435.1 -0.074 0.9416
M9 30.6 435.1 0.070 0.9442
M10 230.6 435.1 0.530 0.5983
M11 -405.6 435.1 -0.932 0.3553
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 688 on 55 degrees of freedom
Multiple R-squared: 0.2232, Adjusted R-squared: 0.06786
F-statistic: 1.437 on 11 and 55 DF, p-value: 0.1832
> 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.75794641 0.48410719 0.2420536
[2,] 0.61735819 0.76528362 0.3826418
[3,] 0.51787227 0.96425545 0.4821277
[4,] 0.38485641 0.76971282 0.6151436
[5,] 0.29897786 0.59795573 0.7010221
[6,] 0.21302076 0.42604152 0.7869792
[7,] 0.15507280 0.31014560 0.8449272
[8,] 0.12888346 0.25776691 0.8711165
[9,] 0.08730328 0.17460655 0.9126967
[10,] 0.12277499 0.24554999 0.8772250
[11,] 0.09477161 0.18954323 0.9052284
[12,] 0.07370042 0.14740084 0.9262996
[13,] 0.06227505 0.12455010 0.9377249
[14,] 0.05811623 0.11623246 0.9418838
[15,] 0.13222662 0.26445324 0.8677734
[16,] 0.10054868 0.20109737 0.8994513
[17,] 0.06671936 0.13343872 0.9332806
[18,] 0.04779400 0.09558801 0.9522060
[19,] 0.07553783 0.15107566 0.9244622
[20,] 0.12015460 0.24030920 0.8798454
[21,] 0.08738318 0.17476635 0.9126168
[22,] 0.11838728 0.23677456 0.8816127
[23,] 0.13723643 0.27447285 0.8627636
[24,] 0.13724683 0.27449365 0.8627532
[25,] 0.14115344 0.28230688 0.8588466
[26,] 0.09746562 0.19493125 0.9025344
[27,] 0.09878172 0.19756345 0.9012183
[28,] 0.08429902 0.16859803 0.9157010
[29,] 0.05873565 0.11747131 0.9412643
[30,] 0.04947901 0.09895802 0.9505210
[31,] 0.03096911 0.06193821 0.9690309
[32,] 0.02891663 0.05783327 0.9710834
[33,] 0.01906339 0.03812679 0.9809366
[34,] 0.01967351 0.03934703 0.9803265
[35,] 0.02249687 0.04499374 0.9775031
[36,] 0.01581010 0.03162019 0.9841899
[37,] 0.02069118 0.04138236 0.9793088
[38,] 0.01612043 0.03224086 0.9838796
> postscript(file="/var/www/rcomp/tmp/1e4md1292715654.ps",horizontal=F,onefile=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/rcomp/tmp/2e4md1292715654.ps",horizontal=F,onefile=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/rcomp/tmp/37vlg1292715654.ps",horizontal=F,onefile=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/rcomp/tmp/47vlg1292715654.ps",horizontal=F,onefile=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/rcomp/tmp/57vlg1292715654.ps",horizontal=F,onefile=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 = 67
Frequency = 1
1 2 3 4 5 6
-195.00000 87.50000 211.00000 338.16667 982.66667 270.00000
7 8 9 10 11 12
-93.83333 513.60000 697.00000 44.00000 416.20000 1251.60000
13 14 15 16 17 18
1393.00000 698.50000 1094.00000 313.16667 590.66667 256.00000
19 20 21 22 23 24
362.16667 282.60000 500.00000 588.00000 71.20000 144.60000
25 26 27 28 29 30
398.00000 -122.50000 108.00000 -385.83333 -657.33333 -219.00000
31 32 33 34 35 36
43.16667 129.60000 -702.00000 946.00000 173.20000 -650.40000
37 38 39 40 41 42
398.00000 497.50000 -597.00000 -20.83333 -674.33333 -521.00000
43 44 45 46 47 48
-303.83333 -620.40000 -71.00000 -685.00000 34.20000 304.60000
49 50 51 52 53 54
-1114.00000 -789.50000 -1068.00000 -677.83333 -1097.33333 -1013.00000
55 56 57 58 59 60
-426.83333 -305.40000 -424.00000 -893.00000 -694.80000 -1050.40000
61 62 63 64 65 66
-880.00000 -371.50000 252.00000 433.16667 855.66667 1227.00000
67
419.16667
> postscript(file="/var/www/rcomp/tmp/6i5lj1292715654.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 -195.00000 NA
1 87.50000 -195.00000
2 211.00000 87.50000
3 338.16667 211.00000
4 982.66667 338.16667
5 270.00000 982.66667
6 -93.83333 270.00000
7 513.60000 -93.83333
8 697.00000 513.60000
9 44.00000 697.00000
10 416.20000 44.00000
11 1251.60000 416.20000
12 1393.00000 1251.60000
13 698.50000 1393.00000
14 1094.00000 698.50000
15 313.16667 1094.00000
16 590.66667 313.16667
17 256.00000 590.66667
18 362.16667 256.00000
19 282.60000 362.16667
20 500.00000 282.60000
21 588.00000 500.00000
22 71.20000 588.00000
23 144.60000 71.20000
24 398.00000 144.60000
25 -122.50000 398.00000
26 108.00000 -122.50000
27 -385.83333 108.00000
28 -657.33333 -385.83333
29 -219.00000 -657.33333
30 43.16667 -219.00000
31 129.60000 43.16667
32 -702.00000 129.60000
33 946.00000 -702.00000
34 173.20000 946.00000
35 -650.40000 173.20000
36 398.00000 -650.40000
37 497.50000 398.00000
38 -597.00000 497.50000
39 -20.83333 -597.00000
40 -674.33333 -20.83333
41 -521.00000 -674.33333
42 -303.83333 -521.00000
43 -620.40000 -303.83333
44 -71.00000 -620.40000
45 -685.00000 -71.00000
46 34.20000 -685.00000
47 304.60000 34.20000
48 -1114.00000 304.60000
49 -789.50000 -1114.00000
50 -1068.00000 -789.50000
51 -677.83333 -1068.00000
52 -1097.33333 -677.83333
53 -1013.00000 -1097.33333
54 -426.83333 -1013.00000
55 -305.40000 -426.83333
56 -424.00000 -305.40000
57 -893.00000 -424.00000
58 -694.80000 -893.00000
59 -1050.40000 -694.80000
60 -880.00000 -1050.40000
61 -371.50000 -880.00000
62 252.00000 -371.50000
63 433.16667 252.00000
64 855.66667 433.16667
65 1227.00000 855.66667
66 419.16667 1227.00000
67 NA 419.16667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 87.50000 -195.00000
[2,] 211.00000 87.50000
[3,] 338.16667 211.00000
[4,] 982.66667 338.16667
[5,] 270.00000 982.66667
[6,] -93.83333 270.00000
[7,] 513.60000 -93.83333
[8,] 697.00000 513.60000
[9,] 44.00000 697.00000
[10,] 416.20000 44.00000
[11,] 1251.60000 416.20000
[12,] 1393.00000 1251.60000
[13,] 698.50000 1393.00000
[14,] 1094.00000 698.50000
[15,] 313.16667 1094.00000
[16,] 590.66667 313.16667
[17,] 256.00000 590.66667
[18,] 362.16667 256.00000
[19,] 282.60000 362.16667
[20,] 500.00000 282.60000
[21,] 588.00000 500.00000
[22,] 71.20000 588.00000
[23,] 144.60000 71.20000
[24,] 398.00000 144.60000
[25,] -122.50000 398.00000
[26,] 108.00000 -122.50000
[27,] -385.83333 108.00000
[28,] -657.33333 -385.83333
[29,] -219.00000 -657.33333
[30,] 43.16667 -219.00000
[31,] 129.60000 43.16667
[32,] -702.00000 129.60000
[33,] 946.00000 -702.00000
[34,] 173.20000 946.00000
[35,] -650.40000 173.20000
[36,] 398.00000 -650.40000
[37,] 497.50000 398.00000
[38,] -597.00000 497.50000
[39,] -20.83333 -597.00000
[40,] -674.33333 -20.83333
[41,] -521.00000 -674.33333
[42,] -303.83333 -521.00000
[43,] -620.40000 -303.83333
[44,] -71.00000 -620.40000
[45,] -685.00000 -71.00000
[46,] 34.20000 -685.00000
[47,] 304.60000 34.20000
[48,] -1114.00000 304.60000
[49,] -789.50000 -1114.00000
[50,] -1068.00000 -789.50000
[51,] -677.83333 -1068.00000
[52,] -1097.33333 -677.83333
[53,] -1013.00000 -1097.33333
[54,] -426.83333 -1013.00000
[55,] -305.40000 -426.83333
[56,] -424.00000 -305.40000
[57,] -893.00000 -424.00000
[58,] -694.80000 -893.00000
[59,] -1050.40000 -694.80000
[60,] -880.00000 -1050.40000
[61,] -371.50000 -880.00000
[62,] 252.00000 -371.50000
[63,] 433.16667 252.00000
[64,] 855.66667 433.16667
[65,] 1227.00000 855.66667
[66,] 419.16667 1227.00000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 87.50000 -195.00000
2 211.00000 87.50000
3 338.16667 211.00000
4 982.66667 338.16667
5 270.00000 982.66667
6 -93.83333 270.00000
7 513.60000 -93.83333
8 697.00000 513.60000
9 44.00000 697.00000
10 416.20000 44.00000
11 1251.60000 416.20000
12 1393.00000 1251.60000
13 698.50000 1393.00000
14 1094.00000 698.50000
15 313.16667 1094.00000
16 590.66667 313.16667
17 256.00000 590.66667
18 362.16667 256.00000
19 282.60000 362.16667
20 500.00000 282.60000
21 588.00000 500.00000
22 71.20000 588.00000
23 144.60000 71.20000
24 398.00000 144.60000
25 -122.50000 398.00000
26 108.00000 -122.50000
27 -385.83333 108.00000
28 -657.33333 -385.83333
29 -219.00000 -657.33333
30 43.16667 -219.00000
31 129.60000 43.16667
32 -702.00000 129.60000
33 946.00000 -702.00000
34 173.20000 946.00000
35 -650.40000 173.20000
36 398.00000 -650.40000
37 497.50000 398.00000
38 -597.00000 497.50000
39 -20.83333 -597.00000
40 -674.33333 -20.83333
41 -521.00000 -674.33333
42 -303.83333 -521.00000
43 -620.40000 -303.83333
44 -71.00000 -620.40000
45 -685.00000 -71.00000
46 34.20000 -685.00000
47 304.60000 34.20000
48 -1114.00000 304.60000
49 -789.50000 -1114.00000
50 -1068.00000 -789.50000
51 -677.83333 -1068.00000
52 -1097.33333 -677.83333
53 -1013.00000 -1097.33333
54 -426.83333 -1013.00000
55 -305.40000 -426.83333
56 -424.00000 -305.40000
57 -893.00000 -424.00000
58 -694.80000 -893.00000
59 -1050.40000 -694.80000
60 -880.00000 -1050.40000
61 -371.50000 -880.00000
62 252.00000 -371.50000
63 433.16667 252.00000
64 855.66667 433.16667
65 1227.00000 855.66667
66 419.16667 1227.00000
> 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/rcomp/tmp/7ae241292715654.ps",horizontal=F,onefile=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/rcomp/tmp/8ae241292715654.ps",horizontal=F,onefile=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/rcomp/tmp/9ae241292715654.ps",horizontal=F,onefile=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/rcomp/tmp/1035j71292715654.ps",horizontal=F,onefile=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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/1176hd1292715654.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/rcomp/tmp/12soyi1292715654.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/rcomp/tmp/136yer1292715654.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/rcomp/tmp/149gux1292715654.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/rcomp/tmp/15vhbl1292715654.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/rcomp/tmp/16g0rr1292715654.tab")
+ }
>
> try(system("convert tmp/1e4md1292715654.ps tmp/1e4md1292715654.png",intern=TRUE))
character(0)
> try(system("convert tmp/2e4md1292715654.ps tmp/2e4md1292715654.png",intern=TRUE))
character(0)
> try(system("convert tmp/37vlg1292715654.ps tmp/37vlg1292715654.png",intern=TRUE))
character(0)
> try(system("convert tmp/47vlg1292715654.ps tmp/47vlg1292715654.png",intern=TRUE))
character(0)
> try(system("convert tmp/57vlg1292715654.ps tmp/57vlg1292715654.png",intern=TRUE))
character(0)
> try(system("convert tmp/6i5lj1292715654.ps tmp/6i5lj1292715654.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ae241292715654.ps tmp/7ae241292715654.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ae241292715654.ps tmp/8ae241292715654.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ae241292715654.ps tmp/9ae241292715654.png",intern=TRUE))
character(0)
> try(system("convert tmp/1035j71292715654.ps tmp/1035j71292715654.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.120 1.400 4.499