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(15044.5,1,14944.2,1,16754.8,1,14254,1,15454.9,1,15644.8,1,14568.3,1,12520.2,1,14803,1,15873.2,1,14755.3,1,12875.1,1,14291.1,1,14205.3,1,15859.4,1,15258.9,1,15498.6,1,15106.5,1,15023.6,1,12083,1,15761.3,1,16943,1,15070.3,1,13659.6,1,14768.9,0,14725.1,0,15998.1,0,15370.6,0,14956.9,0,15469.7,0,15101.8,0,11703.7,0,16283.6,0,16726.5,0,14968.9,0,14861,0,14583.3,0,15305.8,0,17903.9,0,16379.4,0,15420.3,0,17870.5,0,15912.8,0,13866.5,0,17823.2,0,17872,0,17420.4,0,16704.4,0,15991.2,0,16583.6,0,19123.5,0,17838.7,0,17209.4,0,18586.5,0,16258.1,0,15141.6,0,19202.1,0,17746.5,0,19090.1,0,18040.3,0,17515.5,0,17751.8,0,21072.4,0,17170,0,19439.5,0,19795.4,0,17574.9,0,16165.4,0,19464.6,0,19932.1,0,19961.2,0,17343.4,0,18924.2,0,18574.1,0,21350.6,0,18594.6,0,19823.1,0,20844.4,0,19640.2,0,17735.4,0,19813.6,0,22160,0,20664.3,0,17877.4,0,21211.2,0,21423.1,0,21688.7,0,23243.2,0,21490.2,0,22925.8,0,23184.8,0,18562.2,0),dim=c(2,92),dimnames=list(c('Y','X'),1:92))
> y <- array(NA,dim=c(2,92),dimnames=list(c('Y','X'),1:92))
> 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 15044.5 1 1 0 0 0 0 0 0 0 0 0 0
2 14944.2 1 0 1 0 0 0 0 0 0 0 0 0
3 16754.8 1 0 0 1 0 0 0 0 0 0 0 0
4 14254.0 1 0 0 0 1 0 0 0 0 0 0 0
5 15454.9 1 0 0 0 0 1 0 0 0 0 0 0
6 15644.8 1 0 0 0 0 0 1 0 0 0 0 0
7 14568.3 1 0 0 0 0 0 0 1 0 0 0 0
8 12520.2 1 0 0 0 0 0 0 0 1 0 0 0
9 14803.0 1 0 0 0 0 0 0 0 0 1 0 0
10 15873.2 1 0 0 0 0 0 0 0 0 0 1 0
11 14755.3 1 0 0 0 0 0 0 0 0 0 0 1
12 12875.1 1 0 0 0 0 0 0 0 0 0 0 0
13 14291.1 1 1 0 0 0 0 0 0 0 0 0 0
14 14205.3 1 0 1 0 0 0 0 0 0 0 0 0
15 15859.4 1 0 0 1 0 0 0 0 0 0 0 0
16 15258.9 1 0 0 0 1 0 0 0 0 0 0 0
17 15498.6 1 0 0 0 0 1 0 0 0 0 0 0
18 15106.5 1 0 0 0 0 0 1 0 0 0 0 0
19 15023.6 1 0 0 0 0 0 0 1 0 0 0 0
20 12083.0 1 0 0 0 0 0 0 0 1 0 0 0
21 15761.3 1 0 0 0 0 0 0 0 0 1 0 0
22 16943.0 1 0 0 0 0 0 0 0 0 0 1 0
23 15070.3 1 0 0 0 0 0 0 0 0 0 0 1
24 13659.6 1 0 0 0 0 0 0 0 0 0 0 0
25 14768.9 0 1 0 0 0 0 0 0 0 0 0 0
26 14725.1 0 0 1 0 0 0 0 0 0 0 0 0
27 15998.1 0 0 0 1 0 0 0 0 0 0 0 0
28 15370.6 0 0 0 0 1 0 0 0 0 0 0 0
29 14956.9 0 0 0 0 0 1 0 0 0 0 0 0
30 15469.7 0 0 0 0 0 0 1 0 0 0 0 0
31 15101.8 0 0 0 0 0 0 0 1 0 0 0 0
32 11703.7 0 0 0 0 0 0 0 0 1 0 0 0
33 16283.6 0 0 0 0 0 0 0 0 0 1 0 0
34 16726.5 0 0 0 0 0 0 0 0 0 0 1 0
35 14968.9 0 0 0 0 0 0 0 0 0 0 0 1
36 14861.0 0 0 0 0 0 0 0 0 0 0 0 0
37 14583.3 0 1 0 0 0 0 0 0 0 0 0 0
38 15305.8 0 0 1 0 0 0 0 0 0 0 0 0
39 17903.9 0 0 0 1 0 0 0 0 0 0 0 0
40 16379.4 0 0 0 0 1 0 0 0 0 0 0 0
41 15420.3 0 0 0 0 0 1 0 0 0 0 0 0
42 17870.5 0 0 0 0 0 0 1 0 0 0 0 0
43 15912.8 0 0 0 0 0 0 0 1 0 0 0 0
44 13866.5 0 0 0 0 0 0 0 0 1 0 0 0
45 17823.2 0 0 0 0 0 0 0 0 0 1 0 0
46 17872.0 0 0 0 0 0 0 0 0 0 0 1 0
47 17420.4 0 0 0 0 0 0 0 0 0 0 0 1
48 16704.4 0 0 0 0 0 0 0 0 0 0 0 0
49 15991.2 0 1 0 0 0 0 0 0 0 0 0 0
50 16583.6 0 0 1 0 0 0 0 0 0 0 0 0
51 19123.5 0 0 0 1 0 0 0 0 0 0 0 0
52 17838.7 0 0 0 0 1 0 0 0 0 0 0 0
53 17209.4 0 0 0 0 0 1 0 0 0 0 0 0
54 18586.5 0 0 0 0 0 0 1 0 0 0 0 0
55 16258.1 0 0 0 0 0 0 0 1 0 0 0 0
56 15141.6 0 0 0 0 0 0 0 0 1 0 0 0
57 19202.1 0 0 0 0 0 0 0 0 0 1 0 0
58 17746.5 0 0 0 0 0 0 0 0 0 0 1 0
59 19090.1 0 0 0 0 0 0 0 0 0 0 0 1
60 18040.3 0 0 0 0 0 0 0 0 0 0 0 0
61 17515.5 0 1 0 0 0 0 0 0 0 0 0 0
62 17751.8 0 0 1 0 0 0 0 0 0 0 0 0
63 21072.4 0 0 0 1 0 0 0 0 0 0 0 0
64 17170.0 0 0 0 0 1 0 0 0 0 0 0 0
65 19439.5 0 0 0 0 0 1 0 0 0 0 0 0
66 19795.4 0 0 0 0 0 0 1 0 0 0 0 0
67 17574.9 0 0 0 0 0 0 0 1 0 0 0 0
68 16165.4 0 0 0 0 0 0 0 0 1 0 0 0
69 19464.6 0 0 0 0 0 0 0 0 0 1 0 0
70 19932.1 0 0 0 0 0 0 0 0 0 0 1 0
71 19961.2 0 0 0 0 0 0 0 0 0 0 0 1
72 17343.4 0 0 0 0 0 0 0 0 0 0 0 0
73 18924.2 0 1 0 0 0 0 0 0 0 0 0 0
74 18574.1 0 0 1 0 0 0 0 0 0 0 0 0
75 21350.6 0 0 0 1 0 0 0 0 0 0 0 0
76 18594.6 0 0 0 0 1 0 0 0 0 0 0 0
77 19823.1 0 0 0 0 0 1 0 0 0 0 0 0
78 20844.4 0 0 0 0 0 0 1 0 0 0 0 0
79 19640.2 0 0 0 0 0 0 0 1 0 0 0 0
80 17735.4 0 0 0 0 0 0 0 0 1 0 0 0
81 19813.6 0 0 0 0 0 0 0 0 0 1 0 0
82 22160.0 0 0 0 0 0 0 0 0 0 0 1 0
83 20664.3 0 0 0 0 0 0 0 0 0 0 0 1
84 17877.4 0 0 0 0 0 0 0 0 0 0 0 0
85 21211.2 0 1 0 0 0 0 0 0 0 0 0 0
86 21423.1 0 0 1 0 0 0 0 0 0 0 0 0
87 21688.7 0 0 0 1 0 0 0 0 0 0 0 0
88 23243.2 0 0 0 0 1 0 0 0 0 0 0 0
89 21490.2 0 0 0 0 0 1 0 0 0 0 0 0
90 22925.8 0 0 0 0 0 0 1 0 0 0 0 0
91 23184.8 0 0 0 0 0 0 0 1 0 0 0 0
92 18562.2 0 0 0 0 0 0 0 0 1 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
16804.2 -3134.0 520.6 668.5 2698.3 1243.0
M5 M6 M7 M8 M9 M10
1390.9 2259.8 1137.4 -1298.4 1684.3 2270.3
M11
1509.9
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3802.05 -1195.71 -83.52 1020.21 5243.24
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16804.2 788.3 21.317 < 2e-16 ***
X -3134.0 487.7 -6.426 9.15e-09 ***
M1 520.6 1062.6 0.490 0.6255
M2 668.5 1062.6 0.629 0.5311
M3 2698.3 1062.6 2.539 0.0131 *
M4 1243.0 1062.6 1.170 0.2456
M5 1390.9 1062.6 1.309 0.1943
M6 2259.8 1062.6 2.127 0.0366 *
M7 1137.4 1062.6 1.070 0.2877
M8 -1298.4 1062.6 -1.222 0.2253
M9 1684.3 1097.3 1.535 0.1288
M10 2270.3 1097.3 2.069 0.0418 *
M11 1509.9 1097.3 1.376 0.1727
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2053 on 79 degrees of freedom
Multiple R-squared: 0.4537, Adjusted R-squared: 0.3707
F-statistic: 5.467 on 12 and 79 DF, p-value: 1.101e-06
> 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,] 3.657494e-02 7.314988e-02 0.9634251
[2,] 8.924929e-03 1.784986e-02 0.9910751
[3,] 2.471591e-03 4.943183e-03 0.9975284
[4,] 6.224860e-04 1.244972e-03 0.9993775
[5,] 1.503643e-04 3.007286e-04 0.9998496
[6,] 6.894792e-05 1.378958e-04 0.9999311
[7,] 3.710896e-05 7.421793e-05 0.9999629
[8,] 8.599587e-06 1.719917e-05 0.9999914
[9,] 2.982167e-06 5.964333e-06 0.9999970
[10,] 6.872567e-07 1.374513e-06 0.9999993
[11,] 1.601551e-07 3.203102e-07 0.9999998
[12,] 5.044658e-08 1.008932e-07 0.9999999
[13,] 1.821918e-08 3.643835e-08 1.0000000
[14,] 7.746401e-09 1.549280e-08 1.0000000
[15,] 2.534128e-09 5.068257e-09 1.0000000
[16,] 7.495362e-10 1.499072e-09 1.0000000
[17,] 4.962985e-10 9.925970e-10 1.0000000
[18,] 4.867006e-10 9.734013e-10 1.0000000
[19,] 1.386984e-10 2.773968e-10 1.0000000
[20,] 5.847876e-11 1.169575e-10 1.0000000
[21,] 2.401518e-10 4.803037e-10 1.0000000
[22,] 1.174465e-10 2.348930e-10 1.0000000
[23,] 6.722616e-11 1.344523e-10 1.0000000
[24,] 4.107834e-10 8.215669e-10 1.0000000
[25,] 7.161710e-10 1.432342e-09 1.0000000
[26,] 5.487830e-10 1.097566e-09 1.0000000
[27,] 1.435487e-08 2.870974e-08 1.0000000
[28,] 1.320008e-08 2.640016e-08 1.0000000
[29,] 3.227897e-08 6.455794e-08 1.0000000
[30,] 9.311501e-08 1.862300e-07 0.9999999
[31,] 7.136990e-08 1.427398e-07 0.9999999
[32,] 3.100630e-07 6.201260e-07 0.9999997
[33,] 1.309631e-06 2.619261e-06 0.9999987
[34,] 1.657591e-06 3.315182e-06 0.9999983
[35,] 2.457318e-06 4.914636e-06 0.9999975
[36,] 6.620944e-06 1.324189e-05 0.9999934
[37,] 1.347761e-05 2.695523e-05 0.9999865
[38,] 2.351028e-05 4.702057e-05 0.9999765
[39,] 5.816733e-05 1.163347e-04 0.9999418
[40,] 1.375445e-04 2.750889e-04 0.9998625
[41,] 2.937311e-04 5.874622e-04 0.9997063
[42,] 4.607299e-04 9.214599e-04 0.9995393
[43,] 7.007962e-04 1.401592e-03 0.9992992
[44,] 1.565120e-03 3.130239e-03 0.9984349
[45,] 2.448773e-03 4.897546e-03 0.9975512
[46,] 4.011871e-03 8.023741e-03 0.9959881
[47,] 5.618010e-03 1.123602e-02 0.9943820
[48,] 9.172534e-03 1.834507e-02 0.9908275
[49,] 2.023192e-02 4.046385e-02 0.9797681
[50,] 2.682617e-02 5.365233e-02 0.9731738
[51,] 3.595221e-02 7.190441e-02 0.9640478
[52,] 9.399841e-02 1.879968e-01 0.9060016
[53,] 1.109574e-01 2.219148e-01 0.8890426
[54,] 8.513762e-02 1.702752e-01 0.9148624
[55,] 9.143902e-02 1.828780e-01 0.9085610
[56,] 7.870639e-02 1.574128e-01 0.9212936
[57,] 5.096580e-02 1.019316e-01 0.9490342
[58,] 5.736710e-02 1.147342e-01 0.9426329
[59,] 7.180509e-02 1.436102e-01 0.9281949
[60,] 4.642693e-02 9.285386e-02 0.9535731
[61,] 1.860644e-01 3.721289e-01 0.8139356
> postscript(file="/var/www/html/rcomp/tmp/1kpdx1229016028.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/2tpd11229016028.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/3a9hj1229016028.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/4aq6i1229016028.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/560401229016028.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 = 92
Frequency = 1
1 2 3 4 5 6
853.76991 605.58241 386.38241 -659.16759 393.79491 -285.14259
7 8 9 10 11 12
-239.25509 148.45741 -551.47866 -67.26437 -424.76437 -795.06437
13 14 15 16 17 18
100.36991 -133.31759 -509.01759 345.73241 437.49491 -823.44259
19 20 21 22 23 24
216.04491 -288.74259 406.82134 1002.53563 -109.76437 -10.56437
25 26 27 28 29 30
-2555.83997 -2747.52747 -3504.32747 -2676.57747 -3238.21497 -3594.25247
31 32 33 34 35 36
-2839.76497 -3802.05247 -2204.88854 -2347.97425 -3345.17425 -1943.17425
37 38 39 40 41 42
-2741.43997 -2166.82747 -1598.52747 -1667.77747 -2774.81497 -1193.45247
43 44 45 46 47 48
-2028.76497 -1639.25247 -665.28854 -1202.47425 -893.67425 -99.77425
49 50 51 52 53 54
-1333.53997 -889.02747 -378.92747 -208.47747 -985.71497 -477.45247
55 56 57 58 59 60
-1683.46497 -364.15247 713.61146 -1327.97425 776.02575 1236.12575
61 62 63 64 65 66
190.76003 279.17253 1569.97253 -877.17747 1244.38503 731.44753
67 68 69 70 71 72
-366.66497 659.64753 976.11146 857.62575 1647.12575 539.22575
73 74 75 76 77 78
1599.46003 1101.47253 1848.17253 547.42253 1627.98503 1780.44753
79 80 81 82 83 84
1698.63503 2229.64753 1325.11146 3085.52575 2350.22575 1073.22575
85 86 87 88 89 90
3886.46003 3950.47253 2186.27253 5196.02253 3295.08503 3861.84753
91 92
5243.23503 3056.44753
> postscript(file="/var/www/html/rcomp/tmp/6e22b1229016028.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 = 92
Frequency = 1
lag(myerror, k = 1) myerror
0 853.76991 NA
1 605.58241 853.76991
2 386.38241 605.58241
3 -659.16759 386.38241
4 393.79491 -659.16759
5 -285.14259 393.79491
6 -239.25509 -285.14259
7 148.45741 -239.25509
8 -551.47866 148.45741
9 -67.26437 -551.47866
10 -424.76437 -67.26437
11 -795.06437 -424.76437
12 100.36991 -795.06437
13 -133.31759 100.36991
14 -509.01759 -133.31759
15 345.73241 -509.01759
16 437.49491 345.73241
17 -823.44259 437.49491
18 216.04491 -823.44259
19 -288.74259 216.04491
20 406.82134 -288.74259
21 1002.53563 406.82134
22 -109.76437 1002.53563
23 -10.56437 -109.76437
24 -2555.83997 -10.56437
25 -2747.52747 -2555.83997
26 -3504.32747 -2747.52747
27 -2676.57747 -3504.32747
28 -3238.21497 -2676.57747
29 -3594.25247 -3238.21497
30 -2839.76497 -3594.25247
31 -3802.05247 -2839.76497
32 -2204.88854 -3802.05247
33 -2347.97425 -2204.88854
34 -3345.17425 -2347.97425
35 -1943.17425 -3345.17425
36 -2741.43997 -1943.17425
37 -2166.82747 -2741.43997
38 -1598.52747 -2166.82747
39 -1667.77747 -1598.52747
40 -2774.81497 -1667.77747
41 -1193.45247 -2774.81497
42 -2028.76497 -1193.45247
43 -1639.25247 -2028.76497
44 -665.28854 -1639.25247
45 -1202.47425 -665.28854
46 -893.67425 -1202.47425
47 -99.77425 -893.67425
48 -1333.53997 -99.77425
49 -889.02747 -1333.53997
50 -378.92747 -889.02747
51 -208.47747 -378.92747
52 -985.71497 -208.47747
53 -477.45247 -985.71497
54 -1683.46497 -477.45247
55 -364.15247 -1683.46497
56 713.61146 -364.15247
57 -1327.97425 713.61146
58 776.02575 -1327.97425
59 1236.12575 776.02575
60 190.76003 1236.12575
61 279.17253 190.76003
62 1569.97253 279.17253
63 -877.17747 1569.97253
64 1244.38503 -877.17747
65 731.44753 1244.38503
66 -366.66497 731.44753
67 659.64753 -366.66497
68 976.11146 659.64753
69 857.62575 976.11146
70 1647.12575 857.62575
71 539.22575 1647.12575
72 1599.46003 539.22575
73 1101.47253 1599.46003
74 1848.17253 1101.47253
75 547.42253 1848.17253
76 1627.98503 547.42253
77 1780.44753 1627.98503
78 1698.63503 1780.44753
79 2229.64753 1698.63503
80 1325.11146 2229.64753
81 3085.52575 1325.11146
82 2350.22575 3085.52575
83 1073.22575 2350.22575
84 3886.46003 1073.22575
85 3950.47253 3886.46003
86 2186.27253 3950.47253
87 5196.02253 2186.27253
88 3295.08503 5196.02253
89 3861.84753 3295.08503
90 5243.23503 3861.84753
91 3056.44753 5243.23503
92 NA 3056.44753
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 605.58241 853.76991
[2,] 386.38241 605.58241
[3,] -659.16759 386.38241
[4,] 393.79491 -659.16759
[5,] -285.14259 393.79491
[6,] -239.25509 -285.14259
[7,] 148.45741 -239.25509
[8,] -551.47866 148.45741
[9,] -67.26437 -551.47866
[10,] -424.76437 -67.26437
[11,] -795.06437 -424.76437
[12,] 100.36991 -795.06437
[13,] -133.31759 100.36991
[14,] -509.01759 -133.31759
[15,] 345.73241 -509.01759
[16,] 437.49491 345.73241
[17,] -823.44259 437.49491
[18,] 216.04491 -823.44259
[19,] -288.74259 216.04491
[20,] 406.82134 -288.74259
[21,] 1002.53563 406.82134
[22,] -109.76437 1002.53563
[23,] -10.56437 -109.76437
[24,] -2555.83997 -10.56437
[25,] -2747.52747 -2555.83997
[26,] -3504.32747 -2747.52747
[27,] -2676.57747 -3504.32747
[28,] -3238.21497 -2676.57747
[29,] -3594.25247 -3238.21497
[30,] -2839.76497 -3594.25247
[31,] -3802.05247 -2839.76497
[32,] -2204.88854 -3802.05247
[33,] -2347.97425 -2204.88854
[34,] -3345.17425 -2347.97425
[35,] -1943.17425 -3345.17425
[36,] -2741.43997 -1943.17425
[37,] -2166.82747 -2741.43997
[38,] -1598.52747 -2166.82747
[39,] -1667.77747 -1598.52747
[40,] -2774.81497 -1667.77747
[41,] -1193.45247 -2774.81497
[42,] -2028.76497 -1193.45247
[43,] -1639.25247 -2028.76497
[44,] -665.28854 -1639.25247
[45,] -1202.47425 -665.28854
[46,] -893.67425 -1202.47425
[47,] -99.77425 -893.67425
[48,] -1333.53997 -99.77425
[49,] -889.02747 -1333.53997
[50,] -378.92747 -889.02747
[51,] -208.47747 -378.92747
[52,] -985.71497 -208.47747
[53,] -477.45247 -985.71497
[54,] -1683.46497 -477.45247
[55,] -364.15247 -1683.46497
[56,] 713.61146 -364.15247
[57,] -1327.97425 713.61146
[58,] 776.02575 -1327.97425
[59,] 1236.12575 776.02575
[60,] 190.76003 1236.12575
[61,] 279.17253 190.76003
[62,] 1569.97253 279.17253
[63,] -877.17747 1569.97253
[64,] 1244.38503 -877.17747
[65,] 731.44753 1244.38503
[66,] -366.66497 731.44753
[67,] 659.64753 -366.66497
[68,] 976.11146 659.64753
[69,] 857.62575 976.11146
[70,] 1647.12575 857.62575
[71,] 539.22575 1647.12575
[72,] 1599.46003 539.22575
[73,] 1101.47253 1599.46003
[74,] 1848.17253 1101.47253
[75,] 547.42253 1848.17253
[76,] 1627.98503 547.42253
[77,] 1780.44753 1627.98503
[78,] 1698.63503 1780.44753
[79,] 2229.64753 1698.63503
[80,] 1325.11146 2229.64753
[81,] 3085.52575 1325.11146
[82,] 2350.22575 3085.52575
[83,] 1073.22575 2350.22575
[84,] 3886.46003 1073.22575
[85,] 3950.47253 3886.46003
[86,] 2186.27253 3950.47253
[87,] 5196.02253 2186.27253
[88,] 3295.08503 5196.02253
[89,] 3861.84753 3295.08503
[90,] 5243.23503 3861.84753
[91,] 3056.44753 5243.23503
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 605.58241 853.76991
2 386.38241 605.58241
3 -659.16759 386.38241
4 393.79491 -659.16759
5 -285.14259 393.79491
6 -239.25509 -285.14259
7 148.45741 -239.25509
8 -551.47866 148.45741
9 -67.26437 -551.47866
10 -424.76437 -67.26437
11 -795.06437 -424.76437
12 100.36991 -795.06437
13 -133.31759 100.36991
14 -509.01759 -133.31759
15 345.73241 -509.01759
16 437.49491 345.73241
17 -823.44259 437.49491
18 216.04491 -823.44259
19 -288.74259 216.04491
20 406.82134 -288.74259
21 1002.53563 406.82134
22 -109.76437 1002.53563
23 -10.56437 -109.76437
24 -2555.83997 -10.56437
25 -2747.52747 -2555.83997
26 -3504.32747 -2747.52747
27 -2676.57747 -3504.32747
28 -3238.21497 -2676.57747
29 -3594.25247 -3238.21497
30 -2839.76497 -3594.25247
31 -3802.05247 -2839.76497
32 -2204.88854 -3802.05247
33 -2347.97425 -2204.88854
34 -3345.17425 -2347.97425
35 -1943.17425 -3345.17425
36 -2741.43997 -1943.17425
37 -2166.82747 -2741.43997
38 -1598.52747 -2166.82747
39 -1667.77747 -1598.52747
40 -2774.81497 -1667.77747
41 -1193.45247 -2774.81497
42 -2028.76497 -1193.45247
43 -1639.25247 -2028.76497
44 -665.28854 -1639.25247
45 -1202.47425 -665.28854
46 -893.67425 -1202.47425
47 -99.77425 -893.67425
48 -1333.53997 -99.77425
49 -889.02747 -1333.53997
50 -378.92747 -889.02747
51 -208.47747 -378.92747
52 -985.71497 -208.47747
53 -477.45247 -985.71497
54 -1683.46497 -477.45247
55 -364.15247 -1683.46497
56 713.61146 -364.15247
57 -1327.97425 713.61146
58 776.02575 -1327.97425
59 1236.12575 776.02575
60 190.76003 1236.12575
61 279.17253 190.76003
62 1569.97253 279.17253
63 -877.17747 1569.97253
64 1244.38503 -877.17747
65 731.44753 1244.38503
66 -366.66497 731.44753
67 659.64753 -366.66497
68 976.11146 659.64753
69 857.62575 976.11146
70 1647.12575 857.62575
71 539.22575 1647.12575
72 1599.46003 539.22575
73 1101.47253 1599.46003
74 1848.17253 1101.47253
75 547.42253 1848.17253
76 1627.98503 547.42253
77 1780.44753 1627.98503
78 1698.63503 1780.44753
79 2229.64753 1698.63503
80 1325.11146 2229.64753
81 3085.52575 1325.11146
82 2350.22575 3085.52575
83 1073.22575 2350.22575
84 3886.46003 1073.22575
85 3950.47253 3886.46003
86 2186.27253 3950.47253
87 5196.02253 2186.27253
88 3295.08503 5196.02253
89 3861.84753 3295.08503
90 5243.23503 3861.84753
91 3056.44753 5243.23503
> 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/7in0m1229016028.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/8g5xq1229016028.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/9w6y71229016028.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/10a9h01229016028.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/116kc61229016028.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/12ik4w1229016028.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/13ojut1229016028.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/14yilh1229016028.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/15bgtt1229016028.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/16eodu1229016028.tab")
+ }
>
> system("convert tmp/1kpdx1229016028.ps tmp/1kpdx1229016028.png")
> system("convert tmp/2tpd11229016028.ps tmp/2tpd11229016028.png")
> system("convert tmp/3a9hj1229016028.ps tmp/3a9hj1229016028.png")
> system("convert tmp/4aq6i1229016028.ps tmp/4aq6i1229016028.png")
> system("convert tmp/560401229016028.ps tmp/560401229016028.png")
> system("convert tmp/6e22b1229016028.ps tmp/6e22b1229016028.png")
> system("convert tmp/7in0m1229016028.ps tmp/7in0m1229016028.png")
> system("convert tmp/8g5xq1229016028.ps tmp/8g5xq1229016028.png")
> system("convert tmp/9w6y71229016028.ps tmp/9w6y71229016028.png")
> system("convert tmp/10a9h01229016028.ps tmp/10a9h01229016028.png")
>
>
> proc.time()
user system elapsed
2.876 1.647 3.584