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(1.4,2,1.2,2,1,2,1.7,2,2.4,2,2,2,2.1,2,2,2,1.8,2,2.7,2,2.3,2,1.9,2,2,2,2.3,2,2.8,2,2.4,2,2.3,2,2.7,2,2.7,2,2.9,2,3,2,2.2,2,2.3,2,2.8,2.21,2.8,2.25,2.8,2.25,2.2,2.45,2.6,2.5,2.8,2.5,2.5,2.64,2.4,2.75,2.3,2.93,1.9,3,1.7,3.17,2,3.25,2.1,3.39,1.7,3.5,1.8,3.5,1.8,3.65,1.8,3.75,1.3,3.75,1.3,3.9,1.3,4,1.2,4,1.4,4,2.2,4,2.9,4,3.1,4,3.5,4,3.6,4,4.4,4,4.1,4,5.1,4,5.8,4,5.9,4.18,5.4,4.25,5.5,4.25,4.8,3.97,3.2,3.42,2.7,2.75,2.1,2.31,1.9,2,0.6,1.66,0.7,1.31,-0.2,1.09,-1,1,-1.7,1,-0.7,1,-1,1),dim=c(2,69),dimnames=list(c('Y','X'),1:69))
> y <- array(NA,dim=c(2,69),dimnames=list(c('Y','X'),1:69))
> 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 1.4 2.00 1 0 0 0 0 0 0 0 0 0 0
2 1.2 2.00 0 1 0 0 0 0 0 0 0 0 0
3 1.0 2.00 0 0 1 0 0 0 0 0 0 0 0
4 1.7 2.00 0 0 0 1 0 0 0 0 0 0 0
5 2.4 2.00 0 0 0 0 1 0 0 0 0 0 0
6 2.0 2.00 0 0 0 0 0 1 0 0 0 0 0
7 2.1 2.00 0 0 0 0 0 0 1 0 0 0 0
8 2.0 2.00 0 0 0 0 0 0 0 1 0 0 0
9 1.8 2.00 0 0 0 0 0 0 0 0 1 0 0
10 2.7 2.00 0 0 0 0 0 0 0 0 0 1 0
11 2.3 2.00 0 0 0 0 0 0 0 0 0 0 1
12 1.9 2.00 0 0 0 0 0 0 0 0 0 0 0
13 2.0 2.00 1 0 0 0 0 0 0 0 0 0 0
14 2.3 2.00 0 1 0 0 0 0 0 0 0 0 0
15 2.8 2.00 0 0 1 0 0 0 0 0 0 0 0
16 2.4 2.00 0 0 0 1 0 0 0 0 0 0 0
17 2.3 2.00 0 0 0 0 1 0 0 0 0 0 0
18 2.7 2.00 0 0 0 0 0 1 0 0 0 0 0
19 2.7 2.00 0 0 0 0 0 0 1 0 0 0 0
20 2.9 2.00 0 0 0 0 0 0 0 1 0 0 0
21 3.0 2.00 0 0 0 0 0 0 0 0 1 0 0
22 2.2 2.00 0 0 0 0 0 0 0 0 0 1 0
23 2.3 2.00 0 0 0 0 0 0 0 0 0 0 1
24 2.8 2.21 0 0 0 0 0 0 0 0 0 0 0
25 2.8 2.25 1 0 0 0 0 0 0 0 0 0 0
26 2.8 2.25 0 1 0 0 0 0 0 0 0 0 0
27 2.2 2.45 0 0 1 0 0 0 0 0 0 0 0
28 2.6 2.50 0 0 0 1 0 0 0 0 0 0 0
29 2.8 2.50 0 0 0 0 1 0 0 0 0 0 0
30 2.5 2.64 0 0 0 0 0 1 0 0 0 0 0
31 2.4 2.75 0 0 0 0 0 0 1 0 0 0 0
32 2.3 2.93 0 0 0 0 0 0 0 1 0 0 0
33 1.9 3.00 0 0 0 0 0 0 0 0 1 0 0
34 1.7 3.17 0 0 0 0 0 0 0 0 0 1 0
35 2.0 3.25 0 0 0 0 0 0 0 0 0 0 1
36 2.1 3.39 0 0 0 0 0 0 0 0 0 0 0
37 1.7 3.50 1 0 0 0 0 0 0 0 0 0 0
38 1.8 3.50 0 1 0 0 0 0 0 0 0 0 0
39 1.8 3.65 0 0 1 0 0 0 0 0 0 0 0
40 1.8 3.75 0 0 0 1 0 0 0 0 0 0 0
41 1.3 3.75 0 0 0 0 1 0 0 0 0 0 0
42 1.3 3.90 0 0 0 0 0 1 0 0 0 0 0
43 1.3 4.00 0 0 0 0 0 0 1 0 0 0 0
44 1.2 4.00 0 0 0 0 0 0 0 1 0 0 0
45 1.4 4.00 0 0 0 0 0 0 0 0 1 0 0
46 2.2 4.00 0 0 0 0 0 0 0 0 0 1 0
47 2.9 4.00 0 0 0 0 0 0 0 0 0 0 1
48 3.1 4.00 0 0 0 0 0 0 0 0 0 0 0
49 3.5 4.00 1 0 0 0 0 0 0 0 0 0 0
50 3.6 4.00 0 1 0 0 0 0 0 0 0 0 0
51 4.4 4.00 0 0 1 0 0 0 0 0 0 0 0
52 4.1 4.00 0 0 0 1 0 0 0 0 0 0 0
53 5.1 4.00 0 0 0 0 1 0 0 0 0 0 0
54 5.8 4.00 0 0 0 0 0 1 0 0 0 0 0
55 5.9 4.18 0 0 0 0 0 0 1 0 0 0 0
56 5.4 4.25 0 0 0 0 0 0 0 1 0 0 0
57 5.5 4.25 0 0 0 0 0 0 0 0 1 0 0
58 4.8 3.97 0 0 0 0 0 0 0 0 0 1 0
59 3.2 3.42 0 0 0 0 0 0 0 0 0 0 1
60 2.7 2.75 0 0 0 0 0 0 0 0 0 0 0
61 2.1 2.31 1 0 0 0 0 0 0 0 0 0 0
62 1.9 2.00 0 1 0 0 0 0 0 0 0 0 0
63 0.6 1.66 0 0 1 0 0 0 0 0 0 0 0
64 0.7 1.31 0 0 0 1 0 0 0 0 0 0 0
65 -0.2 1.09 0 0 0 0 1 0 0 0 0 0 0
66 -1.0 1.00 0 0 0 0 0 1 0 0 0 0 0
67 -1.7 1.00 0 0 0 0 0 0 1 0 0 0 0
68 -0.7 1.00 0 0 0 0 0 0 0 1 0 0 0
69 -1.0 1.00 0 0 0 0 0 0 0 0 1 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
-0.11377 0.91769 -0.09258 -0.02850 -0.16336 -0.04944
M5 M6 M7 M8 M9 M10
0.05088 -0.04638 -0.20603 -0.17760 -0.27164 0.05500
M11
-0.03873
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.2979 -0.8300 0.2140 0.6275 2.3839
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.11377 0.71670 -0.159 0.874
X 0.91769 0.15556 5.899 2.22e-07 ***
M1 -0.09258 0.75973 -0.122 0.903
M2 -0.02850 0.76009 -0.037 0.970
M3 -0.16336 0.76008 -0.215 0.831
M4 -0.04944 0.76035 -0.065 0.948
M5 0.05088 0.76070 0.067 0.947
M6 -0.04638 0.76038 -0.061 0.952
M7 -0.20603 0.75987 -0.271 0.787
M8 -0.17760 0.75961 -0.234 0.816
M9 -0.27164 0.75955 -0.358 0.722
M10 0.05500 0.79327 0.069 0.945
M11 -0.03873 0.79295 -0.049 0.961
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.254 on 56 degrees of freedom
Multiple R-squared: 0.3926, Adjusted R-squared: 0.2624
F-statistic: 3.016 on 12 and 56 DF, p-value: 0.002552
> 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,] 2.716171e-01 5.432341e-01 0.7283829
[2,] 1.385418e-01 2.770835e-01 0.8614582
[3,] 7.892778e-02 1.578556e-01 0.9210722
[4,] 4.354197e-02 8.708394e-02 0.9564580
[5,] 2.961192e-02 5.922385e-02 0.9703881
[6,] 2.744161e-02 5.488322e-02 0.9725584
[7,] 1.463193e-02 2.926386e-02 0.9853681
[8,] 6.989829e-03 1.397966e-02 0.9930102
[9,] 3.305795e-03 6.611591e-03 0.9966942
[10,] 1.541738e-03 3.083476e-03 0.9984583
[11,] 6.918482e-04 1.383696e-03 0.9993082
[12,] 6.825030e-04 1.365006e-03 0.9993175
[13,] 3.332321e-04 6.664643e-04 0.9996668
[14,] 1.605965e-04 3.211929e-04 0.9998394
[15,] 9.202620e-05 1.840524e-04 0.9999080
[16,] 5.373668e-05 1.074734e-04 0.9999463
[17,] 2.975966e-05 5.951932e-05 0.9999702
[18,] 1.810128e-05 3.620255e-05 0.9999819
[19,] 1.158723e-05 2.317446e-05 0.9999884
[20,] 4.212093e-06 8.424186e-06 0.9999958
[21,] 1.505006e-06 3.010013e-06 0.9999985
[22,] 6.542421e-07 1.308484e-06 0.9999993
[23,] 2.806148e-07 5.612297e-07 0.9999997
[24,] 1.299796e-07 2.599593e-07 0.9999999
[25,] 7.608703e-08 1.521741e-07 0.9999999
[26,] 1.966748e-07 3.933496e-07 0.9999998
[27,] 6.427572e-07 1.285514e-06 0.9999994
[28,] 2.164526e-06 4.329052e-06 0.9999978
[29,] 2.654736e-05 5.309473e-05 0.9999735
[30,] 4.498909e-04 8.997818e-04 0.9995501
[31,] 2.542547e-03 5.085094e-03 0.9974575
[32,] 5.003794e-03 1.000759e-02 0.9949962
[33,] 1.642911e-02 3.285822e-02 0.9835709
[34,] 6.808938e-02 1.361788e-01 0.9319106
[35,] 2.763853e-01 5.527707e-01 0.7236147
[36,] 3.600779e-01 7.201557e-01 0.6399221
[37,] 7.618741e-01 4.762517e-01 0.2381259
[38,] 7.756739e-01 4.486522e-01 0.2243261
> postscript(file="/var/www/html/rcomp/tmp/1h9cs1258717427.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/2fol71258717427.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/328m81258717427.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/4yl8w1258717427.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/5sepq1258717427.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 = 69
Frequency = 1
1 2 3 4 5 6
-0.22902952 -0.49311019 -0.55824737 0.02782962 0.62751430 0.32477065
7 8 9 10 11 12
0.58442052 0.45599095 0.35003067 0.92338569 0.61712280 0.17839062
13 14 15 16 17 18
0.37097048 0.60688981 1.24175263 0.72782962 0.52751430 1.02477065
19 20 21 22 23 24
1.18442052 1.35599095 1.55003067 0.42338569 0.61712280 0.88567564
25 26 27 28 29 30
0.94154789 0.87746722 0.22879196 0.46898443 0.56866912 0.23744882
31 32 33 34 35 36
0.19615275 -0.09746108 -0.46765969 -1.15031203 -0.82999015 -0.89719899
37 38 39 40 41 42
-1.30556507 -1.26964573 -1.27243647 -1.47812852 -2.07844383 -2.11884104
43 44 45 46 47 48
-2.05096021 -2.17938977 -1.88535005 -1.41199503 -0.61825793 -0.45699011
49 50 51 52 53 54
0.03558975 0.07150908 1.00637190 0.59244889 1.49213358 2.28938992
55 56 57 58 59 60
2.38385553 1.79118764 1.98522736 1.21553568 0.21400248 0.29012284
61 62 63 64 65 66
0.18648647 0.20688981 -0.64623265 -0.33896403 -1.13738747 -1.75753899
67 68 69
-2.29788912 -1.32631868 -1.53227896
> postscript(file="/var/www/html/rcomp/tmp/6f3vy1258717427.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.22902952 NA
1 -0.49311019 -0.22902952
2 -0.55824737 -0.49311019
3 0.02782962 -0.55824737
4 0.62751430 0.02782962
5 0.32477065 0.62751430
6 0.58442052 0.32477065
7 0.45599095 0.58442052
8 0.35003067 0.45599095
9 0.92338569 0.35003067
10 0.61712280 0.92338569
11 0.17839062 0.61712280
12 0.37097048 0.17839062
13 0.60688981 0.37097048
14 1.24175263 0.60688981
15 0.72782962 1.24175263
16 0.52751430 0.72782962
17 1.02477065 0.52751430
18 1.18442052 1.02477065
19 1.35599095 1.18442052
20 1.55003067 1.35599095
21 0.42338569 1.55003067
22 0.61712280 0.42338569
23 0.88567564 0.61712280
24 0.94154789 0.88567564
25 0.87746722 0.94154789
26 0.22879196 0.87746722
27 0.46898443 0.22879196
28 0.56866912 0.46898443
29 0.23744882 0.56866912
30 0.19615275 0.23744882
31 -0.09746108 0.19615275
32 -0.46765969 -0.09746108
33 -1.15031203 -0.46765969
34 -0.82999015 -1.15031203
35 -0.89719899 -0.82999015
36 -1.30556507 -0.89719899
37 -1.26964573 -1.30556507
38 -1.27243647 -1.26964573
39 -1.47812852 -1.27243647
40 -2.07844383 -1.47812852
41 -2.11884104 -2.07844383
42 -2.05096021 -2.11884104
43 -2.17938977 -2.05096021
44 -1.88535005 -2.17938977
45 -1.41199503 -1.88535005
46 -0.61825793 -1.41199503
47 -0.45699011 -0.61825793
48 0.03558975 -0.45699011
49 0.07150908 0.03558975
50 1.00637190 0.07150908
51 0.59244889 1.00637190
52 1.49213358 0.59244889
53 2.28938992 1.49213358
54 2.38385553 2.28938992
55 1.79118764 2.38385553
56 1.98522736 1.79118764
57 1.21553568 1.98522736
58 0.21400248 1.21553568
59 0.29012284 0.21400248
60 0.18648647 0.29012284
61 0.20688981 0.18648647
62 -0.64623265 0.20688981
63 -0.33896403 -0.64623265
64 -1.13738747 -0.33896403
65 -1.75753899 -1.13738747
66 -2.29788912 -1.75753899
67 -1.32631868 -2.29788912
68 -1.53227896 -1.32631868
69 NA -1.53227896
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.49311019 -0.22902952
[2,] -0.55824737 -0.49311019
[3,] 0.02782962 -0.55824737
[4,] 0.62751430 0.02782962
[5,] 0.32477065 0.62751430
[6,] 0.58442052 0.32477065
[7,] 0.45599095 0.58442052
[8,] 0.35003067 0.45599095
[9,] 0.92338569 0.35003067
[10,] 0.61712280 0.92338569
[11,] 0.17839062 0.61712280
[12,] 0.37097048 0.17839062
[13,] 0.60688981 0.37097048
[14,] 1.24175263 0.60688981
[15,] 0.72782962 1.24175263
[16,] 0.52751430 0.72782962
[17,] 1.02477065 0.52751430
[18,] 1.18442052 1.02477065
[19,] 1.35599095 1.18442052
[20,] 1.55003067 1.35599095
[21,] 0.42338569 1.55003067
[22,] 0.61712280 0.42338569
[23,] 0.88567564 0.61712280
[24,] 0.94154789 0.88567564
[25,] 0.87746722 0.94154789
[26,] 0.22879196 0.87746722
[27,] 0.46898443 0.22879196
[28,] 0.56866912 0.46898443
[29,] 0.23744882 0.56866912
[30,] 0.19615275 0.23744882
[31,] -0.09746108 0.19615275
[32,] -0.46765969 -0.09746108
[33,] -1.15031203 -0.46765969
[34,] -0.82999015 -1.15031203
[35,] -0.89719899 -0.82999015
[36,] -1.30556507 -0.89719899
[37,] -1.26964573 -1.30556507
[38,] -1.27243647 -1.26964573
[39,] -1.47812852 -1.27243647
[40,] -2.07844383 -1.47812852
[41,] -2.11884104 -2.07844383
[42,] -2.05096021 -2.11884104
[43,] -2.17938977 -2.05096021
[44,] -1.88535005 -2.17938977
[45,] -1.41199503 -1.88535005
[46,] -0.61825793 -1.41199503
[47,] -0.45699011 -0.61825793
[48,] 0.03558975 -0.45699011
[49,] 0.07150908 0.03558975
[50,] 1.00637190 0.07150908
[51,] 0.59244889 1.00637190
[52,] 1.49213358 0.59244889
[53,] 2.28938992 1.49213358
[54,] 2.38385553 2.28938992
[55,] 1.79118764 2.38385553
[56,] 1.98522736 1.79118764
[57,] 1.21553568 1.98522736
[58,] 0.21400248 1.21553568
[59,] 0.29012284 0.21400248
[60,] 0.18648647 0.29012284
[61,] 0.20688981 0.18648647
[62,] -0.64623265 0.20688981
[63,] -0.33896403 -0.64623265
[64,] -1.13738747 -0.33896403
[65,] -1.75753899 -1.13738747
[66,] -2.29788912 -1.75753899
[67,] -1.32631868 -2.29788912
[68,] -1.53227896 -1.32631868
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.49311019 -0.22902952
2 -0.55824737 -0.49311019
3 0.02782962 -0.55824737
4 0.62751430 0.02782962
5 0.32477065 0.62751430
6 0.58442052 0.32477065
7 0.45599095 0.58442052
8 0.35003067 0.45599095
9 0.92338569 0.35003067
10 0.61712280 0.92338569
11 0.17839062 0.61712280
12 0.37097048 0.17839062
13 0.60688981 0.37097048
14 1.24175263 0.60688981
15 0.72782962 1.24175263
16 0.52751430 0.72782962
17 1.02477065 0.52751430
18 1.18442052 1.02477065
19 1.35599095 1.18442052
20 1.55003067 1.35599095
21 0.42338569 1.55003067
22 0.61712280 0.42338569
23 0.88567564 0.61712280
24 0.94154789 0.88567564
25 0.87746722 0.94154789
26 0.22879196 0.87746722
27 0.46898443 0.22879196
28 0.56866912 0.46898443
29 0.23744882 0.56866912
30 0.19615275 0.23744882
31 -0.09746108 0.19615275
32 -0.46765969 -0.09746108
33 -1.15031203 -0.46765969
34 -0.82999015 -1.15031203
35 -0.89719899 -0.82999015
36 -1.30556507 -0.89719899
37 -1.26964573 -1.30556507
38 -1.27243647 -1.26964573
39 -1.47812852 -1.27243647
40 -2.07844383 -1.47812852
41 -2.11884104 -2.07844383
42 -2.05096021 -2.11884104
43 -2.17938977 -2.05096021
44 -1.88535005 -2.17938977
45 -1.41199503 -1.88535005
46 -0.61825793 -1.41199503
47 -0.45699011 -0.61825793
48 0.03558975 -0.45699011
49 0.07150908 0.03558975
50 1.00637190 0.07150908
51 0.59244889 1.00637190
52 1.49213358 0.59244889
53 2.28938992 1.49213358
54 2.38385553 2.28938992
55 1.79118764 2.38385553
56 1.98522736 1.79118764
57 1.21553568 1.98522736
58 0.21400248 1.21553568
59 0.29012284 0.21400248
60 0.18648647 0.29012284
61 0.20688981 0.18648647
62 -0.64623265 0.20688981
63 -0.33896403 -0.64623265
64 -1.13738747 -0.33896403
65 -1.75753899 -1.13738747
66 -2.29788912 -1.75753899
67 -1.32631868 -2.29788912
68 -1.53227896 -1.32631868
> 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/71gn61258717427.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/8rjqv1258717427.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/9tw981258717427.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/10njb21258717427.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/11pz0w1258717427.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/129tey1258717427.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/13x2cx1258717427.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/14qxvz1258717427.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/15738w1258717427.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/16sunf1258717427.tab")
+ }
>
> system("convert tmp/1h9cs1258717427.ps tmp/1h9cs1258717427.png")
> system("convert tmp/2fol71258717427.ps tmp/2fol71258717427.png")
> system("convert tmp/328m81258717427.ps tmp/328m81258717427.png")
> system("convert tmp/4yl8w1258717427.ps tmp/4yl8w1258717427.png")
> system("convert tmp/5sepq1258717427.ps tmp/5sepq1258717427.png")
> system("convert tmp/6f3vy1258717427.ps tmp/6f3vy1258717427.png")
> system("convert tmp/71gn61258717427.ps tmp/71gn61258717427.png")
> system("convert tmp/8rjqv1258717427.ps tmp/8rjqv1258717427.png")
> system("convert tmp/9tw981258717427.ps tmp/9tw981258717427.png")
> system("convert tmp/10njb21258717427.ps tmp/10njb21258717427.png")
>
>
> proc.time()
user system elapsed
2.462 1.503 2.880