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(5.81,0,5.76,0,5.99,0,6.12,0,6.03,0,6.25,0,5.80,0,5.67,0,5.89,0,5.91,0,5.86,0,6.07,0,6.27,0,6.68,0,6.77,0,6.71,0,6.62,0,6.50,0,5.89,0,6.05,0,6.43,0,6.47,0,6.62,0,6.77,0,6.70,0,6.95,0,6.73,0,7.07,0,7.28,0,7.32,0,6.76,0,6.93,0,6.99,0,7.16,0,7.28,0,7.08,0,7.34,0,7.87,0,6.28,1,6.30,1,6.36,1,6.28,1,5.89,1,6.04,1,5.96,1,6.10,1,6.26,1,6.02,1,6.25,1,6.41,1,6.22,1,6.57,1,6.18,1,6.26,1,6.10,1,6.02,1,6.06,1,6.35,1,6.21,1,6.48,1,6.74,1,6.53,1,6.80,1,6.75,1,6.56,1,6.66,1,6.18,1,6.40,1,6.43,1,6.54,1,6.44,1,6.64,1,6.82,1,6.97,1,7.00,1,6.91,1,6.74,1,6.98,1,6.37,1,6.56,1,6.63,1,6.87,1,6.68,1,6.75,1,6.84,1,7.15,1,7.09,1,6.97,1,7.15,1),dim=c(2,89),dimnames=list(c('Y','X'),1:89))
> y <- array(NA,dim=c(2,89),dimnames=list(c('Y','X'),1:89))
> 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 5.81 0 1 0 0 0 0 0 0 0 0 0 0
2 5.76 0 0 1 0 0 0 0 0 0 0 0 0
3 5.99 0 0 0 1 0 0 0 0 0 0 0 0
4 6.12 0 0 0 0 1 0 0 0 0 0 0 0
5 6.03 0 0 0 0 0 1 0 0 0 0 0 0
6 6.25 0 0 0 0 0 0 1 0 0 0 0 0
7 5.80 0 0 0 0 0 0 0 1 0 0 0 0
8 5.67 0 0 0 0 0 0 0 0 1 0 0 0
9 5.89 0 0 0 0 0 0 0 0 0 1 0 0
10 5.91 0 0 0 0 0 0 0 0 0 0 1 0
11 5.86 0 0 0 0 0 0 0 0 0 0 0 1
12 6.07 0 0 0 0 0 0 0 0 0 0 0 0
13 6.27 0 1 0 0 0 0 0 0 0 0 0 0
14 6.68 0 0 1 0 0 0 0 0 0 0 0 0
15 6.77 0 0 0 1 0 0 0 0 0 0 0 0
16 6.71 0 0 0 0 1 0 0 0 0 0 0 0
17 6.62 0 0 0 0 0 1 0 0 0 0 0 0
18 6.50 0 0 0 0 0 0 1 0 0 0 0 0
19 5.89 0 0 0 0 0 0 0 1 0 0 0 0
20 6.05 0 0 0 0 0 0 0 0 1 0 0 0
21 6.43 0 0 0 0 0 0 0 0 0 1 0 0
22 6.47 0 0 0 0 0 0 0 0 0 0 1 0
23 6.62 0 0 0 0 0 0 0 0 0 0 0 1
24 6.77 0 0 0 0 0 0 0 0 0 0 0 0
25 6.70 0 1 0 0 0 0 0 0 0 0 0 0
26 6.95 0 0 1 0 0 0 0 0 0 0 0 0
27 6.73 0 0 0 1 0 0 0 0 0 0 0 0
28 7.07 0 0 0 0 1 0 0 0 0 0 0 0
29 7.28 0 0 0 0 0 1 0 0 0 0 0 0
30 7.32 0 0 0 0 0 0 1 0 0 0 0 0
31 6.76 0 0 0 0 0 0 0 1 0 0 0 0
32 6.93 0 0 0 0 0 0 0 0 1 0 0 0
33 6.99 0 0 0 0 0 0 0 0 0 1 0 0
34 7.16 0 0 0 0 0 0 0 0 0 0 1 0
35 7.28 0 0 0 0 0 0 0 0 0 0 0 1
36 7.08 0 0 0 0 0 0 0 0 0 0 0 0
37 7.34 0 1 0 0 0 0 0 0 0 0 0 0
38 7.87 0 0 1 0 0 0 0 0 0 0 0 0
39 6.28 1 0 0 1 0 0 0 0 0 0 0 0
40 6.30 1 0 0 0 1 0 0 0 0 0 0 0
41 6.36 1 0 0 0 0 1 0 0 0 0 0 0
42 6.28 1 0 0 0 0 0 1 0 0 0 0 0
43 5.89 1 0 0 0 0 0 0 1 0 0 0 0
44 6.04 1 0 0 0 0 0 0 0 1 0 0 0
45 5.96 1 0 0 0 0 0 0 0 0 1 0 0
46 6.10 1 0 0 0 0 0 0 0 0 0 1 0
47 6.26 1 0 0 0 0 0 0 0 0 0 0 1
48 6.02 1 0 0 0 0 0 0 0 0 0 0 0
49 6.25 1 1 0 0 0 0 0 0 0 0 0 0
50 6.41 1 0 1 0 0 0 0 0 0 0 0 0
51 6.22 1 0 0 1 0 0 0 0 0 0 0 0
52 6.57 1 0 0 0 1 0 0 0 0 0 0 0
53 6.18 1 0 0 0 0 1 0 0 0 0 0 0
54 6.26 1 0 0 0 0 0 1 0 0 0 0 0
55 6.10 1 0 0 0 0 0 0 1 0 0 0 0
56 6.02 1 0 0 0 0 0 0 0 1 0 0 0
57 6.06 1 0 0 0 0 0 0 0 0 1 0 0
58 6.35 1 0 0 0 0 0 0 0 0 0 1 0
59 6.21 1 0 0 0 0 0 0 0 0 0 0 1
60 6.48 1 0 0 0 0 0 0 0 0 0 0 0
61 6.74 1 1 0 0 0 0 0 0 0 0 0 0
62 6.53 1 0 1 0 0 0 0 0 0 0 0 0
63 6.80 1 0 0 1 0 0 0 0 0 0 0 0
64 6.75 1 0 0 0 1 0 0 0 0 0 0 0
65 6.56 1 0 0 0 0 1 0 0 0 0 0 0
66 6.66 1 0 0 0 0 0 1 0 0 0 0 0
67 6.18 1 0 0 0 0 0 0 1 0 0 0 0
68 6.40 1 0 0 0 0 0 0 0 1 0 0 0
69 6.43 1 0 0 0 0 0 0 0 0 1 0 0
70 6.54 1 0 0 0 0 0 0 0 0 0 1 0
71 6.44 1 0 0 0 0 0 0 0 0 0 0 1
72 6.64 1 0 0 0 0 0 0 0 0 0 0 0
73 6.82 1 1 0 0 0 0 0 0 0 0 0 0
74 6.97 1 0 1 0 0 0 0 0 0 0 0 0
75 7.00 1 0 0 1 0 0 0 0 0 0 0 0
76 6.91 1 0 0 0 1 0 0 0 0 0 0 0
77 6.74 1 0 0 0 0 1 0 0 0 0 0 0
78 6.98 1 0 0 0 0 0 1 0 0 0 0 0
79 6.37 1 0 0 0 0 0 0 1 0 0 0 0
80 6.56 1 0 0 0 0 0 0 0 1 0 0 0
81 6.63 1 0 0 0 0 0 0 0 0 1 0 0
82 6.87 1 0 0 0 0 0 0 0 0 0 1 0
83 6.68 1 0 0 0 0 0 0 0 0 0 0 1
84 6.75 1 0 0 0 0 0 0 0 0 0 0 0
85 6.84 1 1 0 0 0 0 0 0 0 0 0 0
86 7.15 1 0 1 0 0 0 0 0 0 0 0 0
87 7.09 1 0 0 1 0 0 0 0 0 0 0 0
88 6.97 1 0 0 0 1 0 0 0 0 0 0 0
89 7.15 1 0 0 0 0 1 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
6.56129 -0.02976 0.04984 0.24359 0.06731 0.13231
M5 M6 M7 M8 M9 M10
0.07231 0.06286 -0.40286 -0.30571 -0.20286 -0.05857
M11
-0.06571
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.0449 -0.3188 0.0164 0.2413 1.0651
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.56129 0.17224 38.093 <2e-16 ***
X -0.02976 0.09320 -0.319 0.750
M1 0.04984 0.22440 0.222 0.825
M2 0.24359 0.22440 1.086 0.281
M3 0.06731 0.22435 0.300 0.765
M4 0.13231 0.22435 0.590 0.557
M5 0.07231 0.22435 0.322 0.748
M6 0.06286 0.23165 0.271 0.787
M7 -0.40286 0.23165 -1.739 0.086 .
M8 -0.30571 0.23165 -1.320 0.191
M9 -0.20286 0.23165 -0.876 0.384
M10 -0.05857 0.23165 -0.253 0.801
M11 -0.06571 0.23165 -0.284 0.777
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4334 on 76 degrees of freedom
Multiple R-squared: 0.1658, Adjusted R-squared: 0.03404
F-statistic: 1.258 on 12 and 76 DF, p-value: 0.2610
> 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.9611677 0.0776645141 3.883226e-02
[2,] 0.9565243 0.0869513298 4.347566e-02
[3,] 0.9329820 0.1340359273 6.701796e-02
[4,] 0.9091209 0.1817581970 9.087910e-02
[5,] 0.9058438 0.1883124280 9.415621e-02
[6,] 0.9066010 0.1867980807 9.339904e-02
[7,] 0.9192264 0.1615472724 8.077364e-02
[8,] 0.9451981 0.1096038776 5.480194e-02
[9,] 0.9548613 0.0902773297 4.513866e-02
[10,] 0.9732070 0.0535859972 2.679300e-02
[11,] 0.9862502 0.0274996684 1.374983e-02
[12,] 0.9885170 0.0229660658 1.148303e-02
[13,] 0.9916797 0.0166405789 8.320289e-03
[14,] 0.9962808 0.0074383469 3.719173e-03
[15,] 0.9980460 0.0039079122 1.953956e-03
[16,] 0.9988034 0.0023931390 1.196570e-03
[17,] 0.9994281 0.0011438569 5.719284e-04
[18,] 0.9994876 0.0010248190 5.124095e-04
[19,] 0.9996295 0.0007410928 3.705464e-04
[20,] 0.9997413 0.0005173757 2.586879e-04
[21,] 0.9996896 0.0006208698 3.104349e-04
[22,] 0.9998289 0.0003421923 1.710961e-04
[23,] 0.9999346 0.0001308787 6.543935e-05
[24,] 0.9999222 0.0001555693 7.778465e-05
[25,] 0.9999111 0.0001777992 8.889960e-05
[26,] 0.9998524 0.0002951722 1.475861e-04
[27,] 0.9997772 0.0004456300 2.228150e-04
[28,] 0.9996525 0.0006950761 3.475380e-04
[29,] 0.9994413 0.0011173692 5.586846e-04
[30,] 0.9992796 0.0014408539 7.204270e-04
[31,] 0.9991981 0.0016037096 8.018548e-04
[32,] 0.9986046 0.0027907975 1.395399e-03
[33,] 0.9989245 0.0021510533 1.075527e-03
[34,] 0.9990507 0.0018985475 9.492737e-04
[35,] 0.9990087 0.0019826125 9.913062e-04
[36,] 0.9996555 0.0006889961 3.444981e-04
[37,] 0.9995176 0.0009648198 4.824099e-04
[38,] 0.9998015 0.0003970171 1.985085e-04
[39,] 0.9998853 0.0002293149 1.146575e-04
[40,] 0.9997834 0.0004331285 2.165642e-04
[41,] 0.9998208 0.0003584065 1.792033e-04
[42,] 0.9998714 0.0002572967 1.286484e-04
[43,] 0.9998551 0.0002898389 1.449194e-04
[44,] 0.9998361 0.0003277345 1.638672e-04
[45,] 0.9997075 0.0005849881 2.924941e-04
[46,] 0.9994014 0.0011971251 5.985625e-04
[47,] 0.9997968 0.0004064139 2.032069e-04
[48,] 0.9997115 0.0005770023 2.885012e-04
[49,] 0.9994486 0.0011027977 5.513989e-04
[50,] 0.9996256 0.0007488557 3.744279e-04
[51,] 0.9995348 0.0009304020 4.652010e-04
[52,] 0.9989701 0.0020598937 1.029947e-03
[53,] 0.9975865 0.0048269896 2.413495e-03
[54,] 0.9949017 0.0101966191 5.098310e-03
[55,] 0.9939892 0.0120216947 6.010847e-03
[56,] 0.9888854 0.0222292645 1.111463e-02
[57,] 0.9683532 0.0632936256 3.164681e-02
[58,] 0.9083829 0.1832341401 9.161707e-02
> postscript(file="/var/www/html/rcomp/tmp/1fm0w1290878920.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/28v0h1290878920.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/38v0h1290878920.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/4i4zk1290878920.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/5i4zk1290878920.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 = 89
Frequency = 1
1 2 3 4 5 6
-0.80113026 -1.04488026 -0.63860033 -0.57360033 -0.60360033 -0.37414887
7 8 9 10 11 12
-0.35843459 -0.58557744 -0.46843459 -0.59272030 -0.63557744 -0.49129173
13 14 15 16 17 18
-0.34113026 -0.12488026 0.14139967 0.01639967 -0.01360033 -0.12414887
19 20 21 22 23 24
-0.26843459 -0.20557744 0.07156541 -0.03272030 0.12442256 0.20870827
25 26 27 28 29 30
0.08886974 0.14511974 0.10139967 0.37639967 0.64639967 0.69585113
31 32 33 34 35 36
0.60156541 0.67442256 0.63156541 0.65727970 0.78442256 0.51870827
37 38 39 40 41 42
0.72886974 1.06511974 -0.31883980 -0.36383980 -0.24383980 -0.31438834
43 44 45 46 47 48
-0.23867406 -0.18581692 -0.36867406 -0.37295977 -0.20581692 -0.51153120
49 50 51 52 53 54
-0.33136974 -0.36511974 -0.37883980 -0.09383980 -0.42383980 -0.33438834
55 56 57 58 59 60
-0.02867406 -0.20581692 -0.26867406 -0.12295977 -0.25581692 -0.05153120
61 62 63 64 65 66
0.15863026 -0.24511974 0.20116020 0.08616020 -0.04383980 0.06561166
67 68 69 70 71 72
0.05132594 0.17418308 0.10132594 0.06704023 -0.02581692 0.10846880
73 74 75 76 77 78
0.23863026 0.19488026 0.40116020 0.24616020 0.13616020 0.38561166
79 80 81 82 83 84
0.24132594 0.33418308 0.30132594 0.39704023 0.21418308 0.21846880
85 86 87 88 89
0.25863026 0.37488026 0.49116020 0.30616020 0.54616020
> postscript(file="/var/www/html/rcomp/tmp/6i4zk1290878920.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 = 89
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.80113026 NA
1 -1.04488026 -0.80113026
2 -0.63860033 -1.04488026
3 -0.57360033 -0.63860033
4 -0.60360033 -0.57360033
5 -0.37414887 -0.60360033
6 -0.35843459 -0.37414887
7 -0.58557744 -0.35843459
8 -0.46843459 -0.58557744
9 -0.59272030 -0.46843459
10 -0.63557744 -0.59272030
11 -0.49129173 -0.63557744
12 -0.34113026 -0.49129173
13 -0.12488026 -0.34113026
14 0.14139967 -0.12488026
15 0.01639967 0.14139967
16 -0.01360033 0.01639967
17 -0.12414887 -0.01360033
18 -0.26843459 -0.12414887
19 -0.20557744 -0.26843459
20 0.07156541 -0.20557744
21 -0.03272030 0.07156541
22 0.12442256 -0.03272030
23 0.20870827 0.12442256
24 0.08886974 0.20870827
25 0.14511974 0.08886974
26 0.10139967 0.14511974
27 0.37639967 0.10139967
28 0.64639967 0.37639967
29 0.69585113 0.64639967
30 0.60156541 0.69585113
31 0.67442256 0.60156541
32 0.63156541 0.67442256
33 0.65727970 0.63156541
34 0.78442256 0.65727970
35 0.51870827 0.78442256
36 0.72886974 0.51870827
37 1.06511974 0.72886974
38 -0.31883980 1.06511974
39 -0.36383980 -0.31883980
40 -0.24383980 -0.36383980
41 -0.31438834 -0.24383980
42 -0.23867406 -0.31438834
43 -0.18581692 -0.23867406
44 -0.36867406 -0.18581692
45 -0.37295977 -0.36867406
46 -0.20581692 -0.37295977
47 -0.51153120 -0.20581692
48 -0.33136974 -0.51153120
49 -0.36511974 -0.33136974
50 -0.37883980 -0.36511974
51 -0.09383980 -0.37883980
52 -0.42383980 -0.09383980
53 -0.33438834 -0.42383980
54 -0.02867406 -0.33438834
55 -0.20581692 -0.02867406
56 -0.26867406 -0.20581692
57 -0.12295977 -0.26867406
58 -0.25581692 -0.12295977
59 -0.05153120 -0.25581692
60 0.15863026 -0.05153120
61 -0.24511974 0.15863026
62 0.20116020 -0.24511974
63 0.08616020 0.20116020
64 -0.04383980 0.08616020
65 0.06561166 -0.04383980
66 0.05132594 0.06561166
67 0.17418308 0.05132594
68 0.10132594 0.17418308
69 0.06704023 0.10132594
70 -0.02581692 0.06704023
71 0.10846880 -0.02581692
72 0.23863026 0.10846880
73 0.19488026 0.23863026
74 0.40116020 0.19488026
75 0.24616020 0.40116020
76 0.13616020 0.24616020
77 0.38561166 0.13616020
78 0.24132594 0.38561166
79 0.33418308 0.24132594
80 0.30132594 0.33418308
81 0.39704023 0.30132594
82 0.21418308 0.39704023
83 0.21846880 0.21418308
84 0.25863026 0.21846880
85 0.37488026 0.25863026
86 0.49116020 0.37488026
87 0.30616020 0.49116020
88 0.54616020 0.30616020
89 NA 0.54616020
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.04488026 -0.80113026
[2,] -0.63860033 -1.04488026
[3,] -0.57360033 -0.63860033
[4,] -0.60360033 -0.57360033
[5,] -0.37414887 -0.60360033
[6,] -0.35843459 -0.37414887
[7,] -0.58557744 -0.35843459
[8,] -0.46843459 -0.58557744
[9,] -0.59272030 -0.46843459
[10,] -0.63557744 -0.59272030
[11,] -0.49129173 -0.63557744
[12,] -0.34113026 -0.49129173
[13,] -0.12488026 -0.34113026
[14,] 0.14139967 -0.12488026
[15,] 0.01639967 0.14139967
[16,] -0.01360033 0.01639967
[17,] -0.12414887 -0.01360033
[18,] -0.26843459 -0.12414887
[19,] -0.20557744 -0.26843459
[20,] 0.07156541 -0.20557744
[21,] -0.03272030 0.07156541
[22,] 0.12442256 -0.03272030
[23,] 0.20870827 0.12442256
[24,] 0.08886974 0.20870827
[25,] 0.14511974 0.08886974
[26,] 0.10139967 0.14511974
[27,] 0.37639967 0.10139967
[28,] 0.64639967 0.37639967
[29,] 0.69585113 0.64639967
[30,] 0.60156541 0.69585113
[31,] 0.67442256 0.60156541
[32,] 0.63156541 0.67442256
[33,] 0.65727970 0.63156541
[34,] 0.78442256 0.65727970
[35,] 0.51870827 0.78442256
[36,] 0.72886974 0.51870827
[37,] 1.06511974 0.72886974
[38,] -0.31883980 1.06511974
[39,] -0.36383980 -0.31883980
[40,] -0.24383980 -0.36383980
[41,] -0.31438834 -0.24383980
[42,] -0.23867406 -0.31438834
[43,] -0.18581692 -0.23867406
[44,] -0.36867406 -0.18581692
[45,] -0.37295977 -0.36867406
[46,] -0.20581692 -0.37295977
[47,] -0.51153120 -0.20581692
[48,] -0.33136974 -0.51153120
[49,] -0.36511974 -0.33136974
[50,] -0.37883980 -0.36511974
[51,] -0.09383980 -0.37883980
[52,] -0.42383980 -0.09383980
[53,] -0.33438834 -0.42383980
[54,] -0.02867406 -0.33438834
[55,] -0.20581692 -0.02867406
[56,] -0.26867406 -0.20581692
[57,] -0.12295977 -0.26867406
[58,] -0.25581692 -0.12295977
[59,] -0.05153120 -0.25581692
[60,] 0.15863026 -0.05153120
[61,] -0.24511974 0.15863026
[62,] 0.20116020 -0.24511974
[63,] 0.08616020 0.20116020
[64,] -0.04383980 0.08616020
[65,] 0.06561166 -0.04383980
[66,] 0.05132594 0.06561166
[67,] 0.17418308 0.05132594
[68,] 0.10132594 0.17418308
[69,] 0.06704023 0.10132594
[70,] -0.02581692 0.06704023
[71,] 0.10846880 -0.02581692
[72,] 0.23863026 0.10846880
[73,] 0.19488026 0.23863026
[74,] 0.40116020 0.19488026
[75,] 0.24616020 0.40116020
[76,] 0.13616020 0.24616020
[77,] 0.38561166 0.13616020
[78,] 0.24132594 0.38561166
[79,] 0.33418308 0.24132594
[80,] 0.30132594 0.33418308
[81,] 0.39704023 0.30132594
[82,] 0.21418308 0.39704023
[83,] 0.21846880 0.21418308
[84,] 0.25863026 0.21846880
[85,] 0.37488026 0.25863026
[86,] 0.49116020 0.37488026
[87,] 0.30616020 0.49116020
[88,] 0.54616020 0.30616020
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.04488026 -0.80113026
2 -0.63860033 -1.04488026
3 -0.57360033 -0.63860033
4 -0.60360033 -0.57360033
5 -0.37414887 -0.60360033
6 -0.35843459 -0.37414887
7 -0.58557744 -0.35843459
8 -0.46843459 -0.58557744
9 -0.59272030 -0.46843459
10 -0.63557744 -0.59272030
11 -0.49129173 -0.63557744
12 -0.34113026 -0.49129173
13 -0.12488026 -0.34113026
14 0.14139967 -0.12488026
15 0.01639967 0.14139967
16 -0.01360033 0.01639967
17 -0.12414887 -0.01360033
18 -0.26843459 -0.12414887
19 -0.20557744 -0.26843459
20 0.07156541 -0.20557744
21 -0.03272030 0.07156541
22 0.12442256 -0.03272030
23 0.20870827 0.12442256
24 0.08886974 0.20870827
25 0.14511974 0.08886974
26 0.10139967 0.14511974
27 0.37639967 0.10139967
28 0.64639967 0.37639967
29 0.69585113 0.64639967
30 0.60156541 0.69585113
31 0.67442256 0.60156541
32 0.63156541 0.67442256
33 0.65727970 0.63156541
34 0.78442256 0.65727970
35 0.51870827 0.78442256
36 0.72886974 0.51870827
37 1.06511974 0.72886974
38 -0.31883980 1.06511974
39 -0.36383980 -0.31883980
40 -0.24383980 -0.36383980
41 -0.31438834 -0.24383980
42 -0.23867406 -0.31438834
43 -0.18581692 -0.23867406
44 -0.36867406 -0.18581692
45 -0.37295977 -0.36867406
46 -0.20581692 -0.37295977
47 -0.51153120 -0.20581692
48 -0.33136974 -0.51153120
49 -0.36511974 -0.33136974
50 -0.37883980 -0.36511974
51 -0.09383980 -0.37883980
52 -0.42383980 -0.09383980
53 -0.33438834 -0.42383980
54 -0.02867406 -0.33438834
55 -0.20581692 -0.02867406
56 -0.26867406 -0.20581692
57 -0.12295977 -0.26867406
58 -0.25581692 -0.12295977
59 -0.05153120 -0.25581692
60 0.15863026 -0.05153120
61 -0.24511974 0.15863026
62 0.20116020 -0.24511974
63 0.08616020 0.20116020
64 -0.04383980 0.08616020
65 0.06561166 -0.04383980
66 0.05132594 0.06561166
67 0.17418308 0.05132594
68 0.10132594 0.17418308
69 0.06704023 0.10132594
70 -0.02581692 0.06704023
71 0.10846880 -0.02581692
72 0.23863026 0.10846880
73 0.19488026 0.23863026
74 0.40116020 0.19488026
75 0.24616020 0.40116020
76 0.13616020 0.24616020
77 0.38561166 0.13616020
78 0.24132594 0.38561166
79 0.33418308 0.24132594
80 0.30132594 0.33418308
81 0.39704023 0.30132594
82 0.21418308 0.39704023
83 0.21846880 0.21418308
84 0.25863026 0.21846880
85 0.37488026 0.25863026
86 0.49116020 0.37488026
87 0.30616020 0.49116020
88 0.54616020 0.30616020
> 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/7tdgn1290878920.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/8m5gq1290878920.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/9m5gq1290878920.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/10m5gq1290878920.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/110wvh1290878920.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/12aov21290878920.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/13h7av1290878920.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/143pqj1290878920.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/15oq7p1290878920.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/1698nd1290878920.tab")
+ }
>
> try(system("convert tmp/1fm0w1290878920.ps tmp/1fm0w1290878920.png",intern=TRUE))
character(0)
> try(system("convert tmp/28v0h1290878920.ps tmp/28v0h1290878920.png",intern=TRUE))
character(0)
> try(system("convert tmp/38v0h1290878920.ps tmp/38v0h1290878920.png",intern=TRUE))
character(0)
> try(system("convert tmp/4i4zk1290878920.ps tmp/4i4zk1290878920.png",intern=TRUE))
character(0)
> try(system("convert tmp/5i4zk1290878920.ps tmp/5i4zk1290878920.png",intern=TRUE))
character(0)
> try(system("convert tmp/6i4zk1290878920.ps tmp/6i4zk1290878920.png",intern=TRUE))
character(0)
> try(system("convert tmp/7tdgn1290878920.ps tmp/7tdgn1290878920.png",intern=TRUE))
character(0)
> try(system("convert tmp/8m5gq1290878920.ps tmp/8m5gq1290878920.png",intern=TRUE))
character(0)
> try(system("convert tmp/9m5gq1290878920.ps tmp/9m5gq1290878920.png",intern=TRUE))
character(0)
> try(system("convert tmp/10m5gq1290878920.ps tmp/10m5gq1290878920.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.786 1.615 6.399