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(565464,0,547344,0,554788,0,562325,0,560854,0,555332,0,543599,0,536662,0,542722,0,593530,0,610763,0,612613,0,611324,0,594167,0,595454,0,590865,0,589379,0,584428,0,573100,0,567456,0,569028,0,620735,0,628884,0,628232,0,612117,0,595404,0,597141,0,593408,0,590072,0,579799,0,574205,0,572775,0,572942,0,619567,0,625809,0,619916,0,587625,0,565742,0,557274,0,560576,1,548854,1,531673,1,525919,1,511038,1,498662,1,555362,1,564591,1,541657,1,527070,1,509846,1,514258,1,516922,1,507561,1,492622,1,490243,1,469357,1,477580,1,528379,1,533590,1,517945,1),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> 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
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 565464 0 1 0 0 0 0 0 0 0 0 0 0
2 547344 0 0 1 0 0 0 0 0 0 0 0 0
3 554788 0 0 0 1 0 0 0 0 0 0 0 0
4 562325 0 0 0 0 1 0 0 0 0 0 0 0
5 560854 0 0 0 0 0 1 0 0 0 0 0 0
6 555332 0 0 0 0 0 0 1 0 0 0 0 0
7 543599 0 0 0 0 0 0 0 1 0 0 0 0
8 536662 0 0 0 0 0 0 0 0 1 0 0 0
9 542722 0 0 0 0 0 0 0 0 0 1 0 0
10 593530 0 0 0 0 0 0 0 0 0 0 1 0
11 610763 0 0 0 0 0 0 0 0 0 0 0 1
12 612613 0 0 0 0 0 0 0 0 0 0 0 0
13 611324 0 1 0 0 0 0 0 0 0 0 0 0
14 594167 0 0 1 0 0 0 0 0 0 0 0 0
15 595454 0 0 0 1 0 0 0 0 0 0 0 0
16 590865 0 0 0 0 1 0 0 0 0 0 0 0
17 589379 0 0 0 0 0 1 0 0 0 0 0 0
18 584428 0 0 0 0 0 0 1 0 0 0 0 0
19 573100 0 0 0 0 0 0 0 1 0 0 0 0
20 567456 0 0 0 0 0 0 0 0 1 0 0 0
21 569028 0 0 0 0 0 0 0 0 0 1 0 0
22 620735 0 0 0 0 0 0 0 0 0 0 1 0
23 628884 0 0 0 0 0 0 0 0 0 0 0 1
24 628232 0 0 0 0 0 0 0 0 0 0 0 0
25 612117 0 1 0 0 0 0 0 0 0 0 0 0
26 595404 0 0 1 0 0 0 0 0 0 0 0 0
27 597141 0 0 0 1 0 0 0 0 0 0 0 0
28 593408 0 0 0 0 1 0 0 0 0 0 0 0
29 590072 0 0 0 0 0 1 0 0 0 0 0 0
30 579799 0 0 0 0 0 0 1 0 0 0 0 0
31 574205 0 0 0 0 0 0 0 1 0 0 0 0
32 572775 0 0 0 0 0 0 0 0 1 0 0 0
33 572942 0 0 0 0 0 0 0 0 0 1 0 0
34 619567 0 0 0 0 0 0 0 0 0 0 1 0
35 625809 0 0 0 0 0 0 0 0 0 0 0 1
36 619916 0 0 0 0 0 0 0 0 0 0 0 0
37 587625 0 1 0 0 0 0 0 0 0 0 0 0
38 565742 0 0 1 0 0 0 0 0 0 0 0 0
39 557274 0 0 0 1 0 0 0 0 0 0 0 0
40 560576 1 0 0 0 1 0 0 0 0 0 0 0
41 548854 1 0 0 0 0 1 0 0 0 0 0 0
42 531673 1 0 0 0 0 0 1 0 0 0 0 0
43 525919 1 0 0 0 0 0 0 1 0 0 0 0
44 511038 1 0 0 0 0 0 0 0 1 0 0 0
45 498662 1 0 0 0 0 0 0 0 0 1 0 0
46 555362 1 0 0 0 0 0 0 0 0 0 1 0
47 564591 1 0 0 0 0 0 0 0 0 0 0 1
48 541657 1 0 0 0 0 0 0 0 0 0 0 0
49 527070 1 1 0 0 0 0 0 0 0 0 0 0
50 509846 1 0 1 0 0 0 0 0 0 0 0 0
51 514258 1 0 0 1 0 0 0 0 0 0 0 0
52 516922 1 0 0 0 1 0 0 0 0 0 0 0
53 507561 1 0 0 0 0 1 0 0 0 0 0 0
54 492622 1 0 0 0 0 0 1 0 0 0 0 0
55 490243 1 0 0 0 0 0 0 1 0 0 0 0
56 469357 1 0 0 0 0 0 0 0 1 0 0 0
57 477580 1 0 0 0 0 0 0 0 0 1 0 0
58 528379 1 0 0 0 0 0 0 0 0 0 1 0
59 533590 1 0 0 0 0 0 0 0 0 0 0 1
60 517945 1 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
610131 -65145 -16382 -34601 -33319 -19253
M5 M6 M7 M8 M9 M10
-24729 -35302 -42659 -52615 -51886 -558
M11
8655
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-28552 -15653 3284 11003 34844
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 610130 8594 70.994 < 2e-16 ***
X -65145 5136 -12.684 < 2e-16 ***
M1 -16382 11846 -1.383 0.17325
M2 -34601 11846 -2.921 0.00535 **
M3 -33319 11846 -2.813 0.00715 **
M4 -19253 11802 -1.631 0.10948
M5 -24729 11802 -2.095 0.04155 *
M6 -35302 11802 -2.991 0.00441 **
M7 -42659 11802 -3.615 0.00073 ***
M8 -52615 11802 -4.458 5.11e-05 ***
M9 -51886 11802 -4.397 6.25e-05 ***
M10 -558 11802 -0.047 0.96249
M11 8655 11802 0.733 0.46698
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18660 on 47 degrees of freedom
Multiple R-squared: 0.8292, Adjusted R-squared: 0.7856
F-statistic: 19.01 on 12 and 47 DF, p-value: 4.079e-14
> 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.98774319 0.02451362 0.01225681
[2,] 0.98157838 0.03684323 0.01842162
[3,] 0.97384277 0.05231446 0.02615723
[4,] 0.96520333 0.06959334 0.03479667
[5,] 0.95544899 0.08910201 0.04455101
[6,] 0.93853669 0.12292661 0.06146331
[7,] 0.91876210 0.16247580 0.08123790
[8,] 0.88252895 0.23494209 0.11747105
[9,] 0.85033938 0.29932124 0.14966062
[10,] 0.83146005 0.33707990 0.16853995
[11,] 0.82015743 0.35968513 0.17984257
[12,] 0.81364826 0.37270348 0.18635174
[13,] 0.77122676 0.45754649 0.22877324
[14,] 0.70994035 0.58011930 0.29005965
[15,] 0.62789954 0.74420091 0.37210046
[16,] 0.55205885 0.89588231 0.44794115
[17,] 0.48993278 0.97986556 0.51006722
[18,] 0.42614590 0.85229180 0.57385410
[19,] 0.34590468 0.69180936 0.65409532
[20,] 0.26687432 0.53374864 0.73312568
[21,] 0.25204346 0.50408692 0.74795654
[22,] 0.18269659 0.36539317 0.81730341
[23,] 0.12907611 0.25815222 0.87092389
[24,] 0.09594183 0.19188366 0.90405817
[25,] 0.10879801 0.21759602 0.89120199
[26,] 0.12466826 0.24933652 0.87533174
[27,] 0.14513572 0.29027144 0.85486428
[28,] 0.15504650 0.31009300 0.84495350
[29,] 0.23756473 0.47512945 0.76243527
> postscript(file="/var/www/html/rcomp/tmp/1qylw1229683131.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/2w3o31229683131.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/308rg1229683131.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/4di861229683131.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/503xf1229683131.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 = 60
Frequency = 1
1 2 3 4 5 6
-28284.95455 -28185.55455 -22023.95455 -28552.10909 -24547.90909 -19496.70909
7 8 9 10 11 12
-23872.10909 -20853.50909 -15522.70909 -16042.50909 -8022.30909 2482.49091
13 14 15 16 17 18
17575.04545 18637.44545 18642.04545 -12.10909 3977.09091 9599.29091
19 20 21 22 23 24
5628.89091 9940.49091 10783.29091 11162.49091 10098.69091 18101.49091
25 26 27 28 29 30
18368.04545 19874.44545 20329.04545 2530.89091 4670.09091 4970.29091
31 32 33 34 35 36
6733.89091 15259.49091 14697.29091 9994.49091 7023.69091 9785.49091
37 38 39 40 41 42
-6123.95455 -9787.55455 -19537.95455 34843.66364 28596.86364 21989.06364
43 44 45 46 47 48
23592.66364 18667.26364 5562.06364 10934.26364 10950.46364 -3328.73636
49 50 51 52 53 54
-1534.18182 -538.78182 2590.81818 -8810.33636 -12696.13636 -17061.93636
55 56 57 58 59 60
-12083.33636 -23013.73636 -15519.93636 -16048.73636 -20050.53636 -27040.73636
> postscript(file="/var/www/html/rcomp/tmp/6jpoi1229683131.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -28284.95455 NA
1 -28185.55455 -28284.95455
2 -22023.95455 -28185.55455
3 -28552.10909 -22023.95455
4 -24547.90909 -28552.10909
5 -19496.70909 -24547.90909
6 -23872.10909 -19496.70909
7 -20853.50909 -23872.10909
8 -15522.70909 -20853.50909
9 -16042.50909 -15522.70909
10 -8022.30909 -16042.50909
11 2482.49091 -8022.30909
12 17575.04545 2482.49091
13 18637.44545 17575.04545
14 18642.04545 18637.44545
15 -12.10909 18642.04545
16 3977.09091 -12.10909
17 9599.29091 3977.09091
18 5628.89091 9599.29091
19 9940.49091 5628.89091
20 10783.29091 9940.49091
21 11162.49091 10783.29091
22 10098.69091 11162.49091
23 18101.49091 10098.69091
24 18368.04545 18101.49091
25 19874.44545 18368.04545
26 20329.04545 19874.44545
27 2530.89091 20329.04545
28 4670.09091 2530.89091
29 4970.29091 4670.09091
30 6733.89091 4970.29091
31 15259.49091 6733.89091
32 14697.29091 15259.49091
33 9994.49091 14697.29091
34 7023.69091 9994.49091
35 9785.49091 7023.69091
36 -6123.95455 9785.49091
37 -9787.55455 -6123.95455
38 -19537.95455 -9787.55455
39 34843.66364 -19537.95455
40 28596.86364 34843.66364
41 21989.06364 28596.86364
42 23592.66364 21989.06364
43 18667.26364 23592.66364
44 5562.06364 18667.26364
45 10934.26364 5562.06364
46 10950.46364 10934.26364
47 -3328.73636 10950.46364
48 -1534.18182 -3328.73636
49 -538.78182 -1534.18182
50 2590.81818 -538.78182
51 -8810.33636 2590.81818
52 -12696.13636 -8810.33636
53 -17061.93636 -12696.13636
54 -12083.33636 -17061.93636
55 -23013.73636 -12083.33636
56 -15519.93636 -23013.73636
57 -16048.73636 -15519.93636
58 -20050.53636 -16048.73636
59 -27040.73636 -20050.53636
60 NA -27040.73636
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -28185.55455 -28284.95455
[2,] -22023.95455 -28185.55455
[3,] -28552.10909 -22023.95455
[4,] -24547.90909 -28552.10909
[5,] -19496.70909 -24547.90909
[6,] -23872.10909 -19496.70909
[7,] -20853.50909 -23872.10909
[8,] -15522.70909 -20853.50909
[9,] -16042.50909 -15522.70909
[10,] -8022.30909 -16042.50909
[11,] 2482.49091 -8022.30909
[12,] 17575.04545 2482.49091
[13,] 18637.44545 17575.04545
[14,] 18642.04545 18637.44545
[15,] -12.10909 18642.04545
[16,] 3977.09091 -12.10909
[17,] 9599.29091 3977.09091
[18,] 5628.89091 9599.29091
[19,] 9940.49091 5628.89091
[20,] 10783.29091 9940.49091
[21,] 11162.49091 10783.29091
[22,] 10098.69091 11162.49091
[23,] 18101.49091 10098.69091
[24,] 18368.04545 18101.49091
[25,] 19874.44545 18368.04545
[26,] 20329.04545 19874.44545
[27,] 2530.89091 20329.04545
[28,] 4670.09091 2530.89091
[29,] 4970.29091 4670.09091
[30,] 6733.89091 4970.29091
[31,] 15259.49091 6733.89091
[32,] 14697.29091 15259.49091
[33,] 9994.49091 14697.29091
[34,] 7023.69091 9994.49091
[35,] 9785.49091 7023.69091
[36,] -6123.95455 9785.49091
[37,] -9787.55455 -6123.95455
[38,] -19537.95455 -9787.55455
[39,] 34843.66364 -19537.95455
[40,] 28596.86364 34843.66364
[41,] 21989.06364 28596.86364
[42,] 23592.66364 21989.06364
[43,] 18667.26364 23592.66364
[44,] 5562.06364 18667.26364
[45,] 10934.26364 5562.06364
[46,] 10950.46364 10934.26364
[47,] -3328.73636 10950.46364
[48,] -1534.18182 -3328.73636
[49,] -538.78182 -1534.18182
[50,] 2590.81818 -538.78182
[51,] -8810.33636 2590.81818
[52,] -12696.13636 -8810.33636
[53,] -17061.93636 -12696.13636
[54,] -12083.33636 -17061.93636
[55,] -23013.73636 -12083.33636
[56,] -15519.93636 -23013.73636
[57,] -16048.73636 -15519.93636
[58,] -20050.53636 -16048.73636
[59,] -27040.73636 -20050.53636
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -28185.55455 -28284.95455
2 -22023.95455 -28185.55455
3 -28552.10909 -22023.95455
4 -24547.90909 -28552.10909
5 -19496.70909 -24547.90909
6 -23872.10909 -19496.70909
7 -20853.50909 -23872.10909
8 -15522.70909 -20853.50909
9 -16042.50909 -15522.70909
10 -8022.30909 -16042.50909
11 2482.49091 -8022.30909
12 17575.04545 2482.49091
13 18637.44545 17575.04545
14 18642.04545 18637.44545
15 -12.10909 18642.04545
16 3977.09091 -12.10909
17 9599.29091 3977.09091
18 5628.89091 9599.29091
19 9940.49091 5628.89091
20 10783.29091 9940.49091
21 11162.49091 10783.29091
22 10098.69091 11162.49091
23 18101.49091 10098.69091
24 18368.04545 18101.49091
25 19874.44545 18368.04545
26 20329.04545 19874.44545
27 2530.89091 20329.04545
28 4670.09091 2530.89091
29 4970.29091 4670.09091
30 6733.89091 4970.29091
31 15259.49091 6733.89091
32 14697.29091 15259.49091
33 9994.49091 14697.29091
34 7023.69091 9994.49091
35 9785.49091 7023.69091
36 -6123.95455 9785.49091
37 -9787.55455 -6123.95455
38 -19537.95455 -9787.55455
39 34843.66364 -19537.95455
40 28596.86364 34843.66364
41 21989.06364 28596.86364
42 23592.66364 21989.06364
43 18667.26364 23592.66364
44 5562.06364 18667.26364
45 10934.26364 5562.06364
46 10950.46364 10934.26364
47 -3328.73636 10950.46364
48 -1534.18182 -3328.73636
49 -538.78182 -1534.18182
50 2590.81818 -538.78182
51 -8810.33636 2590.81818
52 -12696.13636 -8810.33636
53 -17061.93636 -12696.13636
54 -12083.33636 -17061.93636
55 -23013.73636 -12083.33636
56 -15519.93636 -23013.73636
57 -16048.73636 -15519.93636
58 -20050.53636 -16048.73636
59 -27040.73636 -20050.53636
> 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/706121229683131.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/8ro8t1229683131.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/91ybh1229683131.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/10vycg1229683131.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/119qx61229683131.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/123mdm1229683131.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/131j6x1229683131.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/14z2gj1229683131.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/15cqjg1229683131.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/16iqwt1229683131.tab")
+ }
>
> system("convert tmp/1qylw1229683131.ps tmp/1qylw1229683131.png")
> system("convert tmp/2w3o31229683131.ps tmp/2w3o31229683131.png")
> system("convert tmp/308rg1229683131.ps tmp/308rg1229683131.png")
> system("convert tmp/4di861229683131.ps tmp/4di861229683131.png")
> system("convert tmp/503xf1229683131.ps tmp/503xf1229683131.png")
> system("convert tmp/6jpoi1229683131.ps tmp/6jpoi1229683131.png")
> system("convert tmp/706121229683131.ps tmp/706121229683131.png")
> system("convert tmp/8ro8t1229683131.ps tmp/8ro8t1229683131.png")
> system("convert tmp/91ybh1229683131.ps tmp/91ybh1229683131.png")
> system("convert tmp/10vycg1229683131.ps tmp/10vycg1229683131.png")
>
>
> proc.time()
user system elapsed
2.426 1.567 3.031