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(921365,0,987921,0,1132614,0,1332224,0,1418133,0,1411549,0,1695920,0,1636173,0,1539653,0,1395314,0,1127575,0,1036076,0,989236,0,1008380,0,1207763,0,1368839,0,1469798,0,1498721,0,1761769,0,1653214,0,1599104,0,1421179,0,1163995,0,1037735,0,1015407,0,1039210,0,1258049,0,1469445,0,1552346,0,1549144,0,1785895,0,1662335,0,1629440,0,1467430,0,1202209,0,1076982,0,1039367,1,1063449,1,1335135,1,1491602,1,1591972,1,1641248,1,1898849,1,1798580,1,1762444,1,1622044,1,1368955,1,1262973,1,1195650,1,1269530,1,1479279,1,1607819,1,1712466,1,1721766,1,1949843,1,1821326,1,1757802,1,1590367,1,1260647,1,1149235,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 = '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 t
1 921365 0 1 0 0 0 0 0 0 0 0 0 0 1
2 987921 0 0 1 0 0 0 0 0 0 0 0 0 2
3 1132614 0 0 0 1 0 0 0 0 0 0 0 0 3
4 1332224 0 0 0 0 1 0 0 0 0 0 0 0 4
5 1418133 0 0 0 0 0 1 0 0 0 0 0 0 5
6 1411549 0 0 0 0 0 0 1 0 0 0 0 0 6
7 1695920 0 0 0 0 0 0 0 1 0 0 0 0 7
8 1636173 0 0 0 0 0 0 0 0 1 0 0 0 8
9 1539653 0 0 0 0 0 0 0 0 0 1 0 0 9
10 1395314 0 0 0 0 0 0 0 0 0 0 1 0 10
11 1127575 0 0 0 0 0 0 0 0 0 0 0 1 11
12 1036076 0 0 0 0 0 0 0 0 0 0 0 0 12
13 989236 0 1 0 0 0 0 0 0 0 0 0 0 13
14 1008380 0 0 1 0 0 0 0 0 0 0 0 0 14
15 1207763 0 0 0 1 0 0 0 0 0 0 0 0 15
16 1368839 0 0 0 0 1 0 0 0 0 0 0 0 16
17 1469798 0 0 0 0 0 1 0 0 0 0 0 0 17
18 1498721 0 0 0 0 0 0 1 0 0 0 0 0 18
19 1761769 0 0 0 0 0 0 0 1 0 0 0 0 19
20 1653214 0 0 0 0 0 0 0 0 1 0 0 0 20
21 1599104 0 0 0 0 0 0 0 0 0 1 0 0 21
22 1421179 0 0 0 0 0 0 0 0 0 0 1 0 22
23 1163995 0 0 0 0 0 0 0 0 0 0 0 1 23
24 1037735 0 0 0 0 0 0 0 0 0 0 0 0 24
25 1015407 0 1 0 0 0 0 0 0 0 0 0 0 25
26 1039210 0 0 1 0 0 0 0 0 0 0 0 0 26
27 1258049 0 0 0 1 0 0 0 0 0 0 0 0 27
28 1469445 0 0 0 0 1 0 0 0 0 0 0 0 28
29 1552346 0 0 0 0 0 1 0 0 0 0 0 0 29
30 1549144 0 0 0 0 0 0 1 0 0 0 0 0 30
31 1785895 0 0 0 0 0 0 0 1 0 0 0 0 31
32 1662335 0 0 0 0 0 0 0 0 1 0 0 0 32
33 1629440 0 0 0 0 0 0 0 0 0 1 0 0 33
34 1467430 0 0 0 0 0 0 0 0 0 0 1 0 34
35 1202209 0 0 0 0 0 0 0 0 0 0 0 1 35
36 1076982 0 0 0 0 0 0 0 0 0 0 0 0 36
37 1039367 1 1 0 0 0 0 0 0 0 0 0 0 37
38 1063449 1 0 1 0 0 0 0 0 0 0 0 0 38
39 1335135 1 0 0 1 0 0 0 0 0 0 0 0 39
40 1491602 1 0 0 0 1 0 0 0 0 0 0 0 40
41 1591972 1 0 0 0 0 1 0 0 0 0 0 0 41
42 1641248 1 0 0 0 0 0 1 0 0 0 0 0 42
43 1898849 1 0 0 0 0 0 0 1 0 0 0 0 43
44 1798580 1 0 0 0 0 0 0 0 1 0 0 0 44
45 1762444 1 0 0 0 0 0 0 0 0 1 0 0 45
46 1622044 1 0 0 0 0 0 0 0 0 0 1 0 46
47 1368955 1 0 0 0 0 0 0 0 0 0 0 1 47
48 1262973 1 0 0 0 0 0 0 0 0 0 0 0 48
49 1195650 1 1 0 0 0 0 0 0 0 0 0 0 49
50 1269530 1 0 1 0 0 0 0 0 0 0 0 0 50
51 1479279 1 0 0 1 0 0 0 0 0 0 0 0 51
52 1607819 1 0 0 0 1 0 0 0 0 0 0 0 52
53 1712466 1 0 0 0 0 1 0 0 0 0 0 0 53
54 1721766 1 0 0 0 0 0 1 0 0 0 0 0 54
55 1949843 1 0 0 0 0 0 0 1 0 0 0 0 55
56 1821326 1 0 0 0 0 0 0 0 1 0 0 0 56
57 1757802 1 0 0 0 0 0 0 0 0 1 0 0 57
58 1590367 1 0 0 0 0 0 0 0 0 0 1 0 58
59 1260647 1 0 0 0 0 0 0 0 0 0 0 1 59
60 1149235 1 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
952333.9 52433.0 -37833.4 -209.7 204791.1 372339.6
M5 M6 M7 M8 M9 M10
463427.6 475100.9 725201.3 617202.4 556696.2 394405.1
M11 t
115945.3 3869.3
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-88351.0 -15227.4 255.6 24390.9 72482.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 952333.9 23529.0 40.475 < 2e-16 ***
X 52433.0 21129.7 2.481 0.0168 *
M1 -37833.4 26228.3 -1.442 0.1559
M2 -209.7 26078.9 -0.008 0.9936
M3 204791.1 25943.1 7.894 4.24e-10 ***
M4 372339.6 25820.9 14.420 < 2e-16 ***
M5 463427.6 25712.6 18.023 < 2e-16 ***
M6 475100.9 25618.4 18.545 < 2e-16 ***
M7 725201.3 25538.4 28.397 < 2e-16 ***
M8 617202.4 25472.7 24.230 < 2e-16 ***
M9 556696.2 25421.5 21.899 < 2e-16 ***
M10 394405.1 25384.9 15.537 < 2e-16 ***
M11 115945.3 25362.9 4.571 3.64e-05 ***
t 3869.3 610.0 6.343 8.83e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 40090 on 46 degrees of freedom
Multiple R-squared: 0.9832, Adjusted R-squared: 0.9784
F-statistic: 206.5 on 13 and 46 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,] 7.391226e-02 1.478245e-01 0.9260877
[2,] 5.054057e-02 1.010811e-01 0.9494594
[3,] 1.815742e-02 3.631483e-02 0.9818426
[4,] 1.546405e-02 3.092810e-02 0.9845360
[5,] 5.493936e-03 1.098787e-02 0.9945061
[6,] 2.829100e-03 5.658200e-03 0.9971709
[7,] 1.042422e-03 2.084844e-03 0.9989576
[8,] 1.251261e-03 2.502522e-03 0.9987487
[9,] 4.400530e-04 8.801060e-04 0.9995599
[10,] 2.160736e-04 4.321472e-04 0.9997839
[11,] 1.091745e-04 2.183489e-04 0.9998908
[12,] 2.222142e-04 4.444284e-04 0.9997778
[13,] 1.575300e-04 3.150599e-04 0.9998425
[14,] 6.138863e-05 1.227773e-04 0.9999386
[15,] 2.504689e-05 5.009379e-05 0.9999750
[16,] 4.640129e-05 9.280257e-05 0.9999536
[17,] 1.576943e-05 3.153886e-05 0.9999842
[18,] 5.154096e-06 1.030819e-05 0.9999948
[19,] 1.618973e-06 3.237946e-06 0.9999984
[20,] 6.735490e-07 1.347098e-06 0.9999993
[21,] 4.405646e-07 8.811293e-07 0.9999996
[22,] 2.182265e-06 4.364530e-06 0.9999978
[23,] 3.545836e-05 7.091672e-05 0.9999645
[24,] 6.411985e-05 1.282397e-04 0.9999359
[25,] 4.228176e-04 8.456352e-04 0.9995772
[26,] 3.796498e-03 7.592996e-03 0.9962035
[27,] 1.530460e-02 3.060921e-02 0.9846954
> postscript(file="/var/www/html/rcomp/tmp/1qlfe1261310877.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/25j7v1261310877.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/36wmg1261310877.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/429en1261310877.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/57q9v1261310877.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 7
2995.256 28058.256 -36118.744 -7926.544 -16974.744 -39101.344 -8699.944
8 9 10 11 12 13 14
35682.656 -4200.344 9882.456 16734.056 37311.056 24435.222 2086.222
15 16 17 18 19 20 21
-7400.778 -17742.578 -11740.778 1639.622 10718.022 6292.622 8819.622
22 23 24 25 26 27 28
-10683.578 6723.022 -7460.978 4175.189 -13514.811 -3545.811 36432.389
29 30 31 32 33 34 35
24376.189 5631.589 -11587.011 -31017.411 -7275.411 -10863.611 -1494.011
36 37 38 39 40 41 42
-14645.011 -70728.817 -88139.817 -25323.817 -40274.617 -34861.817 -1128.417
43 44 45 46 47 48 49
2502.983 6363.583 26864.583 44886.383 66387.983 72481.983 39123.150
50 51 52 53 54 55 56
71510.150 72389.150 29511.350 39201.150 32958.550 7065.950 -17321.450
57 58 59 60
-24208.450 -33221.650 -88351.050 -87687.050
> postscript(file="/var/www/html/rcomp/tmp/6w3fd1261310877.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 2995.256 NA
1 28058.256 2995.256
2 -36118.744 28058.256
3 -7926.544 -36118.744
4 -16974.744 -7926.544
5 -39101.344 -16974.744
6 -8699.944 -39101.344
7 35682.656 -8699.944
8 -4200.344 35682.656
9 9882.456 -4200.344
10 16734.056 9882.456
11 37311.056 16734.056
12 24435.222 37311.056
13 2086.222 24435.222
14 -7400.778 2086.222
15 -17742.578 -7400.778
16 -11740.778 -17742.578
17 1639.622 -11740.778
18 10718.022 1639.622
19 6292.622 10718.022
20 8819.622 6292.622
21 -10683.578 8819.622
22 6723.022 -10683.578
23 -7460.978 6723.022
24 4175.189 -7460.978
25 -13514.811 4175.189
26 -3545.811 -13514.811
27 36432.389 -3545.811
28 24376.189 36432.389
29 5631.589 24376.189
30 -11587.011 5631.589
31 -31017.411 -11587.011
32 -7275.411 -31017.411
33 -10863.611 -7275.411
34 -1494.011 -10863.611
35 -14645.011 -1494.011
36 -70728.817 -14645.011
37 -88139.817 -70728.817
38 -25323.817 -88139.817
39 -40274.617 -25323.817
40 -34861.817 -40274.617
41 -1128.417 -34861.817
42 2502.983 -1128.417
43 6363.583 2502.983
44 26864.583 6363.583
45 44886.383 26864.583
46 66387.983 44886.383
47 72481.983 66387.983
48 39123.150 72481.983
49 71510.150 39123.150
50 72389.150 71510.150
51 29511.350 72389.150
52 39201.150 29511.350
53 32958.550 39201.150
54 7065.950 32958.550
55 -17321.450 7065.950
56 -24208.450 -17321.450
57 -33221.650 -24208.450
58 -88351.050 -33221.650
59 -87687.050 -88351.050
60 NA -87687.050
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 28058.256 2995.256
[2,] -36118.744 28058.256
[3,] -7926.544 -36118.744
[4,] -16974.744 -7926.544
[5,] -39101.344 -16974.744
[6,] -8699.944 -39101.344
[7,] 35682.656 -8699.944
[8,] -4200.344 35682.656
[9,] 9882.456 -4200.344
[10,] 16734.056 9882.456
[11,] 37311.056 16734.056
[12,] 24435.222 37311.056
[13,] 2086.222 24435.222
[14,] -7400.778 2086.222
[15,] -17742.578 -7400.778
[16,] -11740.778 -17742.578
[17,] 1639.622 -11740.778
[18,] 10718.022 1639.622
[19,] 6292.622 10718.022
[20,] 8819.622 6292.622
[21,] -10683.578 8819.622
[22,] 6723.022 -10683.578
[23,] -7460.978 6723.022
[24,] 4175.189 -7460.978
[25,] -13514.811 4175.189
[26,] -3545.811 -13514.811
[27,] 36432.389 -3545.811
[28,] 24376.189 36432.389
[29,] 5631.589 24376.189
[30,] -11587.011 5631.589
[31,] -31017.411 -11587.011
[32,] -7275.411 -31017.411
[33,] -10863.611 -7275.411
[34,] -1494.011 -10863.611
[35,] -14645.011 -1494.011
[36,] -70728.817 -14645.011
[37,] -88139.817 -70728.817
[38,] -25323.817 -88139.817
[39,] -40274.617 -25323.817
[40,] -34861.817 -40274.617
[41,] -1128.417 -34861.817
[42,] 2502.983 -1128.417
[43,] 6363.583 2502.983
[44,] 26864.583 6363.583
[45,] 44886.383 26864.583
[46,] 66387.983 44886.383
[47,] 72481.983 66387.983
[48,] 39123.150 72481.983
[49,] 71510.150 39123.150
[50,] 72389.150 71510.150
[51,] 29511.350 72389.150
[52,] 39201.150 29511.350
[53,] 32958.550 39201.150
[54,] 7065.950 32958.550
[55,] -17321.450 7065.950
[56,] -24208.450 -17321.450
[57,] -33221.650 -24208.450
[58,] -88351.050 -33221.650
[59,] -87687.050 -88351.050
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 28058.256 2995.256
2 -36118.744 28058.256
3 -7926.544 -36118.744
4 -16974.744 -7926.544
5 -39101.344 -16974.744
6 -8699.944 -39101.344
7 35682.656 -8699.944
8 -4200.344 35682.656
9 9882.456 -4200.344
10 16734.056 9882.456
11 37311.056 16734.056
12 24435.222 37311.056
13 2086.222 24435.222
14 -7400.778 2086.222
15 -17742.578 -7400.778
16 -11740.778 -17742.578
17 1639.622 -11740.778
18 10718.022 1639.622
19 6292.622 10718.022
20 8819.622 6292.622
21 -10683.578 8819.622
22 6723.022 -10683.578
23 -7460.978 6723.022
24 4175.189 -7460.978
25 -13514.811 4175.189
26 -3545.811 -13514.811
27 36432.389 -3545.811
28 24376.189 36432.389
29 5631.589 24376.189
30 -11587.011 5631.589
31 -31017.411 -11587.011
32 -7275.411 -31017.411
33 -10863.611 -7275.411
34 -1494.011 -10863.611
35 -14645.011 -1494.011
36 -70728.817 -14645.011
37 -88139.817 -70728.817
38 -25323.817 -88139.817
39 -40274.617 -25323.817
40 -34861.817 -40274.617
41 -1128.417 -34861.817
42 2502.983 -1128.417
43 6363.583 2502.983
44 26864.583 6363.583
45 44886.383 26864.583
46 66387.983 44886.383
47 72481.983 66387.983
48 39123.150 72481.983
49 71510.150 39123.150
50 72389.150 71510.150
51 29511.350 72389.150
52 39201.150 29511.350
53 32958.550 39201.150
54 7065.950 32958.550
55 -17321.450 7065.950
56 -24208.450 -17321.450
57 -33221.650 -24208.450
58 -88351.050 -33221.650
59 -87687.050 -88351.050
> 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/7rt371261310877.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/8svvl1261310877.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/92wkh1261310877.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/103a801261310877.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/11o99g1261310877.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/122rpx1261310877.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/13tq3x1261310877.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/14318t1261310877.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/15afo81261310877.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/16nqkp1261310877.tab")
+ }
>
> try(system("convert tmp/1qlfe1261310877.ps tmp/1qlfe1261310877.png",intern=TRUE))
character(0)
> try(system("convert tmp/25j7v1261310877.ps tmp/25j7v1261310877.png",intern=TRUE))
character(0)
> try(system("convert tmp/36wmg1261310877.ps tmp/36wmg1261310877.png",intern=TRUE))
character(0)
> try(system("convert tmp/429en1261310877.ps tmp/429en1261310877.png",intern=TRUE))
character(0)
> try(system("convert tmp/57q9v1261310877.ps tmp/57q9v1261310877.png",intern=TRUE))
character(0)
> try(system("convert tmp/6w3fd1261310877.ps tmp/6w3fd1261310877.png",intern=TRUE))
character(0)
> try(system("convert tmp/7rt371261310877.ps tmp/7rt371261310877.png",intern=TRUE))
character(0)
> try(system("convert tmp/8svvl1261310877.ps tmp/8svvl1261310877.png",intern=TRUE))
character(0)
> try(system("convert tmp/92wkh1261310877.ps tmp/92wkh1261310877.png",intern=TRUE))
character(0)
> try(system("convert tmp/103a801261310877.ps tmp/103a801261310877.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.323 1.516 3.476