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(519164,0.9,517009,1.3,509933,1.4,509127,1.5,500857,1.1,506971,1.6,569323,1.5,579714,1.6,577992,1.7,565464,1.6,547344,1.7,554788,1.6,562325,1.6,560854,1.3,555332,1.1,543599,1.6,536662,1.9,542722,1.6,593530,1.7,610763,1.6,612613,1.4,611324,2.1,594167,1.9,595454,1.7,590865,1.8,589379,2,584428,2.5,573100,2.1,567456,2.1,569028,2.3,620735,2.4,628884,2.4,628232,2.3,612117,1.7,595404,2,597141,2.3,593408,2,590072,2,579799,1.3,574205,1.7,572775,1.9,572942,1.7,619567,1.6,625809,1.7,619916,1.8,587625,1.9,565742,1.9,557274,1.9,560576,2,548854,2.1,531673,1.9,525919,1.9,511038,1.3,498662,1.3,555362,1.4,564591,1.2,541657,1.3,527070,1.8,509846,2.2,514258,2.6,516922,2.8,507561,3.1,492622,3.9,490243,3.7,469357,4.6,477580,5.1,528379,5.2,533590,4.9,517945,5.1,506174,4.8,501866,3.9,516141,3.5),dim=c(2,72),dimnames=list(c('TWIB','GI
'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('TWIB','GI
'),1:72))
> 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 = 'Do not include Seasonal 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
TWIB GI\r
1 519164 0.9
2 517009 1.3
3 509933 1.4
4 509127 1.5
5 500857 1.1
6 506971 1.6
7 569323 1.5
8 579714 1.6
9 577992 1.7
10 565464 1.6
11 547344 1.7
12 554788 1.6
13 562325 1.6
14 560854 1.3
15 555332 1.1
16 543599 1.6
17 536662 1.9
18 542722 1.6
19 593530 1.7
20 610763 1.6
21 612613 1.4
22 611324 2.1
23 594167 1.9
24 595454 1.7
25 590865 1.8
26 589379 2.0
27 584428 2.5
28 573100 2.1
29 567456 2.1
30 569028 2.3
31 620735 2.4
32 628884 2.4
33 628232 2.3
34 612117 1.7
35 595404 2.0
36 597141 2.3
37 593408 2.0
38 590072 2.0
39 579799 1.3
40 574205 1.7
41 572775 1.9
42 572942 1.7
43 619567 1.6
44 625809 1.7
45 619916 1.8
46 587625 1.9
47 565742 1.9
48 557274 1.9
49 560576 2.0
50 548854 2.1
51 531673 1.9
52 525919 1.9
53 511038 1.3
54 498662 1.3
55 555362 1.4
56 564591 1.2
57 541657 1.3
58 527070 1.8
59 509846 2.2
60 514258 2.6
61 516922 2.8
62 507561 3.1
63 492622 3.9
64 490243 3.7
65 469357 4.6
66 477580 5.1
67 528379 5.2
68 533590 4.9
69 517945 5.1
70 506174 4.8
71 501866 3.9
72 516141 3.5
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `GI\r`
590737 -15986
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-72295 -28781 1058 30131 76514
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 590737 10340 57.129 < 2e-16 ***
`GI\r` -15986 4313 -3.707 0.000416 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 37680 on 70 degrees of freedom
Multiple R-squared: 0.1641, Adjusted R-squared: 0.1521
F-statistic: 13.74 on 1 and 70 DF, p-value: 0.0004162
> 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.01691351 0.03382702 0.98308649
[2,] 0.00358412 0.00716824 0.99641588
[3,] 0.23767205 0.47534410 0.76232795
[4,] 0.38453804 0.76907608 0.61546196
[5,] 0.35422487 0.70844975 0.64577513
[6,] 0.27819485 0.55638971 0.72180515
[7,] 0.19999774 0.39999547 0.80000226
[8,] 0.13693036 0.27386072 0.86306964
[9,] 0.09525962 0.19051924 0.90474038
[10,] 0.09695437 0.19390873 0.90304563
[11,] 0.10578079 0.21156158 0.89421921
[12,] 0.07466779 0.14933557 0.92533221
[13,] 0.06757542 0.13515084 0.93242458
[14,] 0.04705913 0.09411827 0.95294087
[15,] 0.06124036 0.12248071 0.93875964
[16,] 0.13237985 0.26475971 0.86762015
[17,] 0.26226807 0.52453614 0.73773193
[18,] 0.25472368 0.50944737 0.74527632
[19,] 0.21361946 0.42723892 0.78638054
[20,] 0.19384690 0.38769381 0.80615310
[21,] 0.15770203 0.31540407 0.84229797
[22,] 0.12183209 0.24366418 0.87816791
[23,] 0.11754186 0.23508372 0.88245814
[24,] 0.08998374 0.17996748 0.91001626
[25,] 0.06906461 0.13812922 0.93093539
[26,] 0.05570324 0.11140648 0.94429676
[27,] 0.06663176 0.13326351 0.93336824
[28,] 0.09910579 0.19821158 0.90089421
[29,] 0.15247216 0.30494431 0.84752784
[30,] 0.20641221 0.41282442 0.79358779
[31,] 0.19271967 0.38543934 0.80728033
[32,] 0.19506727 0.39013453 0.80493273
[33,] 0.18600565 0.37201130 0.81399435
[34,] 0.17434011 0.34868022 0.82565989
[35,] 0.16674212 0.33348424 0.83325788
[36,] 0.13481309 0.26962618 0.86518691
[37,] 0.11001358 0.22002717 0.88998642
[38,] 0.08676045 0.17352089 0.91323955
[39,] 0.21114051 0.42228102 0.78885949
[40,] 0.51045726 0.97908547 0.48954274
[41,] 0.83658597 0.32682805 0.16341403
[42,] 0.91321926 0.17356148 0.08678074
[43,] 0.92742297 0.14515406 0.07257703
[44,] 0.93229042 0.13541916 0.06770958
[45,] 0.94958110 0.10083781 0.05041890
[46,] 0.95741371 0.08517259 0.04258629
[47,] 0.95182333 0.09635335 0.04817667
[48,] 0.94523786 0.10952427 0.05476214
[49,] 0.93990860 0.12018280 0.06009140
[50,] 0.96471192 0.07057617 0.03528808
[51,] 0.95779903 0.08440195 0.04220097
[52,] 0.97437543 0.05124913 0.02562457
[53,] 0.97049071 0.05901859 0.02950929
[54,] 0.96419062 0.07161876 0.03580938
[55,] 0.95939888 0.08120224 0.04060112
[56,] 0.95534969 0.08930063 0.04465031
[57,] 0.95039520 0.09920960 0.04960480
[58,] 0.93631825 0.12736350 0.06368175
[59,] 0.91155854 0.17688292 0.08844146
[60,] 0.86925462 0.26149076 0.13074538
[61,] 0.92496820 0.15006360 0.07503180
[62,] 0.98471615 0.03056769 0.01528385
[63,] 0.95155226 0.09689548 0.04844774
> postscript(file="/var/www/html/rcomp/tmp/10sgb1258743279.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/2i1l71258743279.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/336ck1258743279.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/4bzr11258743279.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/55dpe1258743279.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 = 72
Frequency = 1
1 2 3 4 5 6
-57185.6672 -52946.1291 -58423.4946 -57630.8601 -72295.3981 -58188.2256
7 8 9 10 11 12
2565.1399 14554.7744 14431.4089 304.7744 -16216.5911 -10371.2256
13 14 15 16 17 18
-2834.2256 -9101.1291 -17820.3981 -21560.2256 -23701.3221 -22437.2256
19 20 21 22 23 24
29969.4089 45603.7744 44256.5054 54157.9469 33803.6779 31893.4089
25 26 27 28 29 30
28903.0434 30614.3124 33656.4849 15933.9469 10289.9469 15059.2159
31 32 33 34 35 36
68364.8504 76513.8504 74263.2159 48556.4089 36639.3124 43172.2159
37 38 39 40 41 42
34643.3124 31307.3124 9843.8709 10644.4089 12411.6779 9381.4089
43 44 45 46 47 48
54407.7744 62248.4089 57954.0434 27261.6779 5378.6779 -3089.3221
49 50 51 52 53 54
1811.3124 -8312.0531 -28690.3221 -34444.3221 -58917.1291 -71293.1291
55 56 57 58 59 60
-12994.4946 -6962.7636 -28298.1291 -34891.9566 -45721.4186 -34914.8805
61 62 63 64 65 66
-29053.6115 -33618.7080 -35768.6320 -41344.9010 -47843.1904 -31627.0179
67 68 69 70 71 72
20770.6166 21185.7131 8737.9821 -7828.9214 -26524.6320 -18644.1700
> postscript(file="/var/www/html/rcomp/tmp/6ci1m1258743279.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -57185.6672 NA
1 -52946.1291 -57185.6672
2 -58423.4946 -52946.1291
3 -57630.8601 -58423.4946
4 -72295.3981 -57630.8601
5 -58188.2256 -72295.3981
6 2565.1399 -58188.2256
7 14554.7744 2565.1399
8 14431.4089 14554.7744
9 304.7744 14431.4089
10 -16216.5911 304.7744
11 -10371.2256 -16216.5911
12 -2834.2256 -10371.2256
13 -9101.1291 -2834.2256
14 -17820.3981 -9101.1291
15 -21560.2256 -17820.3981
16 -23701.3221 -21560.2256
17 -22437.2256 -23701.3221
18 29969.4089 -22437.2256
19 45603.7744 29969.4089
20 44256.5054 45603.7744
21 54157.9469 44256.5054
22 33803.6779 54157.9469
23 31893.4089 33803.6779
24 28903.0434 31893.4089
25 30614.3124 28903.0434
26 33656.4849 30614.3124
27 15933.9469 33656.4849
28 10289.9469 15933.9469
29 15059.2159 10289.9469
30 68364.8504 15059.2159
31 76513.8504 68364.8504
32 74263.2159 76513.8504
33 48556.4089 74263.2159
34 36639.3124 48556.4089
35 43172.2159 36639.3124
36 34643.3124 43172.2159
37 31307.3124 34643.3124
38 9843.8709 31307.3124
39 10644.4089 9843.8709
40 12411.6779 10644.4089
41 9381.4089 12411.6779
42 54407.7744 9381.4089
43 62248.4089 54407.7744
44 57954.0434 62248.4089
45 27261.6779 57954.0434
46 5378.6779 27261.6779
47 -3089.3221 5378.6779
48 1811.3124 -3089.3221
49 -8312.0531 1811.3124
50 -28690.3221 -8312.0531
51 -34444.3221 -28690.3221
52 -58917.1291 -34444.3221
53 -71293.1291 -58917.1291
54 -12994.4946 -71293.1291
55 -6962.7636 -12994.4946
56 -28298.1291 -6962.7636
57 -34891.9566 -28298.1291
58 -45721.4186 -34891.9566
59 -34914.8805 -45721.4186
60 -29053.6115 -34914.8805
61 -33618.7080 -29053.6115
62 -35768.6320 -33618.7080
63 -41344.9010 -35768.6320
64 -47843.1904 -41344.9010
65 -31627.0179 -47843.1904
66 20770.6166 -31627.0179
67 21185.7131 20770.6166
68 8737.9821 21185.7131
69 -7828.9214 8737.9821
70 -26524.6320 -7828.9214
71 -18644.1700 -26524.6320
72 NA -18644.1700
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -52946.1291 -57185.6672
[2,] -58423.4946 -52946.1291
[3,] -57630.8601 -58423.4946
[4,] -72295.3981 -57630.8601
[5,] -58188.2256 -72295.3981
[6,] 2565.1399 -58188.2256
[7,] 14554.7744 2565.1399
[8,] 14431.4089 14554.7744
[9,] 304.7744 14431.4089
[10,] -16216.5911 304.7744
[11,] -10371.2256 -16216.5911
[12,] -2834.2256 -10371.2256
[13,] -9101.1291 -2834.2256
[14,] -17820.3981 -9101.1291
[15,] -21560.2256 -17820.3981
[16,] -23701.3221 -21560.2256
[17,] -22437.2256 -23701.3221
[18,] 29969.4089 -22437.2256
[19,] 45603.7744 29969.4089
[20,] 44256.5054 45603.7744
[21,] 54157.9469 44256.5054
[22,] 33803.6779 54157.9469
[23,] 31893.4089 33803.6779
[24,] 28903.0434 31893.4089
[25,] 30614.3124 28903.0434
[26,] 33656.4849 30614.3124
[27,] 15933.9469 33656.4849
[28,] 10289.9469 15933.9469
[29,] 15059.2159 10289.9469
[30,] 68364.8504 15059.2159
[31,] 76513.8504 68364.8504
[32,] 74263.2159 76513.8504
[33,] 48556.4089 74263.2159
[34,] 36639.3124 48556.4089
[35,] 43172.2159 36639.3124
[36,] 34643.3124 43172.2159
[37,] 31307.3124 34643.3124
[38,] 9843.8709 31307.3124
[39,] 10644.4089 9843.8709
[40,] 12411.6779 10644.4089
[41,] 9381.4089 12411.6779
[42,] 54407.7744 9381.4089
[43,] 62248.4089 54407.7744
[44,] 57954.0434 62248.4089
[45,] 27261.6779 57954.0434
[46,] 5378.6779 27261.6779
[47,] -3089.3221 5378.6779
[48,] 1811.3124 -3089.3221
[49,] -8312.0531 1811.3124
[50,] -28690.3221 -8312.0531
[51,] -34444.3221 -28690.3221
[52,] -58917.1291 -34444.3221
[53,] -71293.1291 -58917.1291
[54,] -12994.4946 -71293.1291
[55,] -6962.7636 -12994.4946
[56,] -28298.1291 -6962.7636
[57,] -34891.9566 -28298.1291
[58,] -45721.4186 -34891.9566
[59,] -34914.8805 -45721.4186
[60,] -29053.6115 -34914.8805
[61,] -33618.7080 -29053.6115
[62,] -35768.6320 -33618.7080
[63,] -41344.9010 -35768.6320
[64,] -47843.1904 -41344.9010
[65,] -31627.0179 -47843.1904
[66,] 20770.6166 -31627.0179
[67,] 21185.7131 20770.6166
[68,] 8737.9821 21185.7131
[69,] -7828.9214 8737.9821
[70,] -26524.6320 -7828.9214
[71,] -18644.1700 -26524.6320
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -52946.1291 -57185.6672
2 -58423.4946 -52946.1291
3 -57630.8601 -58423.4946
4 -72295.3981 -57630.8601
5 -58188.2256 -72295.3981
6 2565.1399 -58188.2256
7 14554.7744 2565.1399
8 14431.4089 14554.7744
9 304.7744 14431.4089
10 -16216.5911 304.7744
11 -10371.2256 -16216.5911
12 -2834.2256 -10371.2256
13 -9101.1291 -2834.2256
14 -17820.3981 -9101.1291
15 -21560.2256 -17820.3981
16 -23701.3221 -21560.2256
17 -22437.2256 -23701.3221
18 29969.4089 -22437.2256
19 45603.7744 29969.4089
20 44256.5054 45603.7744
21 54157.9469 44256.5054
22 33803.6779 54157.9469
23 31893.4089 33803.6779
24 28903.0434 31893.4089
25 30614.3124 28903.0434
26 33656.4849 30614.3124
27 15933.9469 33656.4849
28 10289.9469 15933.9469
29 15059.2159 10289.9469
30 68364.8504 15059.2159
31 76513.8504 68364.8504
32 74263.2159 76513.8504
33 48556.4089 74263.2159
34 36639.3124 48556.4089
35 43172.2159 36639.3124
36 34643.3124 43172.2159
37 31307.3124 34643.3124
38 9843.8709 31307.3124
39 10644.4089 9843.8709
40 12411.6779 10644.4089
41 9381.4089 12411.6779
42 54407.7744 9381.4089
43 62248.4089 54407.7744
44 57954.0434 62248.4089
45 27261.6779 57954.0434
46 5378.6779 27261.6779
47 -3089.3221 5378.6779
48 1811.3124 -3089.3221
49 -8312.0531 1811.3124
50 -28690.3221 -8312.0531
51 -34444.3221 -28690.3221
52 -58917.1291 -34444.3221
53 -71293.1291 -58917.1291
54 -12994.4946 -71293.1291
55 -6962.7636 -12994.4946
56 -28298.1291 -6962.7636
57 -34891.9566 -28298.1291
58 -45721.4186 -34891.9566
59 -34914.8805 -45721.4186
60 -29053.6115 -34914.8805
61 -33618.7080 -29053.6115
62 -35768.6320 -33618.7080
63 -41344.9010 -35768.6320
64 -47843.1904 -41344.9010
65 -31627.0179 -47843.1904
66 20770.6166 -31627.0179
67 21185.7131 20770.6166
68 8737.9821 21185.7131
69 -7828.9214 8737.9821
70 -26524.6320 -7828.9214
71 -18644.1700 -26524.6320
> 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/7ethf1258743279.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/8rszu1258743279.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/9t3041258743279.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/10nqkz1258743279.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/11q80w1258743279.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/12uxq21258743279.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/13fjob1258743279.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/14cnj51258743279.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/15y1fe1258743279.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/16ur4y1258743279.tab")
+ }
>
> system("convert tmp/10sgb1258743279.ps tmp/10sgb1258743279.png")
> system("convert tmp/2i1l71258743279.ps tmp/2i1l71258743279.png")
> system("convert tmp/336ck1258743279.ps tmp/336ck1258743279.png")
> system("convert tmp/4bzr11258743279.ps tmp/4bzr11258743279.png")
> system("convert tmp/55dpe1258743279.ps tmp/55dpe1258743279.png")
> system("convert tmp/6ci1m1258743279.ps tmp/6ci1m1258743279.png")
> system("convert tmp/7ethf1258743279.ps tmp/7ethf1258743279.png")
> system("convert tmp/8rszu1258743279.ps tmp/8rszu1258743279.png")
> system("convert tmp/9t3041258743279.ps tmp/9t3041258743279.png")
> system("convert tmp/10nqkz1258743279.ps tmp/10nqkz1258743279.png")
>
>
> proc.time()
user system elapsed
2.576 1.578 3.043