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(7.3,0,7.1,0,6.9,0,6.8,0,7.5,0,7.6,0,7.8,0,8,0,8.1,0,8.2,0,8.3,0,8.2,0,8,0,7.9,0,7.6,0,7.6,0,8.2,0,8.3,0,8.4,0,8.4,0,8.4,0,8.6,0,8.9,0,8.8,0,8.3,0,7.5,0,7.2,0,7.5,0,8.8,0,9.3,0,9.3,0,8.7,0,8.2,0,8.3,0,8.5,0,8.6,0,8.6,0,8.2,0,8.1,0,8,0,8.6,0,8.7,0,8.8,0,8.5,0,8.4,0,8.5,0,8.7,0,8.7,0,8.6,0,8.5,0,8.3,0,8.1,0,8.2,0,8.1,0,8.1,0,7.9,0,7.9,0,7.9,0,8,0,8,0,7.9,0,8,0,7.7,1,7.2,1,7.5,1,7.3,1,7,1,7,1,7,1,7.2,1,7.3,1,7.1,1,6.8,1,6.6,1,6.2,1,6.2,1,6.8,1,6.9,1,6.8,1),dim=c(2,79),dimnames=list(c('y','x'),1:79))
> y <- array(NA,dim=c(2,79),dimnames=list(c('y','x'),1:79))
> 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 = 'Do not include Seasonal 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
1 7.3 0
2 7.1 0
3 6.9 0
4 6.8 0
5 7.5 0
6 7.6 0
7 7.8 0
8 8.0 0
9 8.1 0
10 8.2 0
11 8.3 0
12 8.2 0
13 8.0 0
14 7.9 0
15 7.6 0
16 7.6 0
17 8.2 0
18 8.3 0
19 8.4 0
20 8.4 0
21 8.4 0
22 8.6 0
23 8.9 0
24 8.8 0
25 8.3 0
26 7.5 0
27 7.2 0
28 7.5 0
29 8.8 0
30 9.3 0
31 9.3 0
32 8.7 0
33 8.2 0
34 8.3 0
35 8.5 0
36 8.6 0
37 8.6 0
38 8.2 0
39 8.1 0
40 8.0 0
41 8.6 0
42 8.7 0
43 8.8 0
44 8.5 0
45 8.4 0
46 8.5 0
47 8.7 0
48 8.7 0
49 8.6 0
50 8.5 0
51 8.3 0
52 8.1 0
53 8.2 0
54 8.1 0
55 8.1 0
56 7.9 0
57 7.9 0
58 7.9 0
59 8.0 0
60 8.0 0
61 7.9 0
62 8.0 0
63 7.7 1
64 7.2 1
65 7.5 1
66 7.3 1
67 7.0 1
68 7.0 1
69 7.0 1
70 7.2 1
71 7.3 1
72 7.1 1
73 6.8 1
74 6.6 1
75 6.2 1
76 6.2 1
77 6.8 1
78 6.9 1
79 6.8 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x
8.168 -1.191
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.36774 -0.22211 0.03226 0.33226 1.13226
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.16774 0.06306 129.517 < 2e-16 ***
x -1.19127 0.13595 -8.763 3.44e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4966 on 77 degrees of freedom
Multiple R-squared: 0.4993, Adjusted R-squared: 0.4928
F-statistic: 76.79 on 1 and 77 DF, p-value: 3.437e-13
> 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.3959551 7.919102e-01 6.040449e-01
[2,] 0.4411752 8.823505e-01 5.588248e-01
[3,] 0.5548246 8.903508e-01 4.451754e-01
[4,] 0.6973109 6.053781e-01 3.026891e-01
[5,] 0.7889814 4.220373e-01 2.110186e-01
[6,] 0.8518098 2.963805e-01 1.481902e-01
[7,] 0.8966650 2.066701e-01 1.033350e-01
[8,] 0.9004817 1.990365e-01 9.951826e-02
[9,] 0.8765193 2.469614e-01 1.234807e-01
[10,] 0.8422274 3.155452e-01 1.577726e-01
[11,] 0.8168657 3.662686e-01 1.831343e-01
[12,] 0.7940654 4.118692e-01 2.059346e-01
[13,] 0.7907499 4.185003e-01 2.092501e-01
[14,] 0.7992666 4.014668e-01 2.007334e-01
[15,] 0.8187897 3.624207e-01 1.812103e-01
[16,] 0.8273570 3.452860e-01 1.726430e-01
[17,] 0.8283105 3.433789e-01 1.716895e-01
[18,] 0.8595971 2.808059e-01 1.404029e-01
[19,] 0.9294261 1.411478e-01 7.057392e-02
[20,] 0.9535971 9.280579e-02 4.640290e-02
[21,] 0.9396522 1.206956e-01 6.034779e-02
[22,] 0.9522295 9.554093e-02 4.777047e-02
[23,] 0.9844697 3.106054e-02 1.553027e-02
[24,] 0.9906482 1.870369e-02 9.351846e-03
[25,] 0.9942523 1.149548e-02 5.747738e-03
[26,] 0.9995392 9.215359e-04 4.607679e-04
[27,] 0.9999768 4.632905e-05 2.316453e-05
[28,] 0.9999791 4.174601e-05 2.087301e-05
[29,] 0.9999594 8.110065e-05 4.055032e-05
[30,] 0.9999239 1.522678e-04 7.613388e-05
[31,] 0.9998861 2.278625e-04 1.139312e-04
[32,] 0.9998609 2.782997e-04 1.391499e-04
[33,] 0.9998309 3.382634e-04 1.691317e-04
[34,] 0.9996878 6.243435e-04 3.121718e-04
[35,] 0.9994634 1.073129e-03 5.365645e-04
[36,] 0.9991916 1.616878e-03 8.084389e-04
[37,] 0.9990132 1.973505e-03 9.867523e-04
[38,] 0.9990601 1.879857e-03 9.399287e-04
[39,] 0.9993720 1.256076e-03 6.280381e-04
[40,] 0.9991240 1.751918e-03 8.759591e-04
[41,] 0.9986147 2.770591e-03 1.385296e-03
[42,] 0.9981561 3.687842e-03 1.843921e-03
[43,] 0.9985678 2.864490e-03 1.432245e-03
[44,] 0.9990308 1.938388e-03 9.691941e-04
[45,] 0.9992105 1.578998e-03 7.894988e-04
[46,] 0.9992238 1.552327e-03 7.761635e-04
[47,] 0.9988536 2.292775e-03 1.146388e-03
[48,] 0.9979515 4.096965e-03 2.048483e-03
[49,] 0.9967400 6.519991e-03 3.259995e-03
[50,] 0.9945279 1.094426e-02 5.472130e-03
[51,] 0.9911200 1.775998e-02 8.879989e-03
[52,] 0.9856568 2.868634e-02 1.434317e-02
[53,] 0.9773856 4.522871e-02 2.261436e-02
[54,] 0.9652545 6.949093e-02 3.474546e-02
[55,] 0.9468158 1.063684e-01 5.318422e-02
[56,] 0.9209599 1.580802e-01 7.904012e-02
[57,] 0.8875426 2.249147e-01 1.124574e-01
[58,] 0.8415871 3.168258e-01 1.584129e-01
[59,] 0.9026940 1.946120e-01 9.730600e-02
[60,] 0.8769538 2.460923e-01 1.230462e-01
[61,] 0.9083810 1.832380e-01 9.161900e-02
[62,] 0.9080766 1.838468e-01 9.192338e-02
[63,] 0.8683355 2.633290e-01 1.316645e-01
[64,] 0.8157493 3.685014e-01 1.842507e-01
[65,] 0.7498267 5.003467e-01 2.501733e-01
[66,] 0.7351204 5.297592e-01 2.648796e-01
[67,] 0.8056505 3.886989e-01 1.943495e-01
[68,] 0.8205527 3.588946e-01 1.794473e-01
[69,] 0.7300172 5.399655e-01 2.699828e-01
[70,] 0.5709599 8.580802e-01 4.290401e-01
> postscript(file="/var/www/html/freestat/rcomp/tmp/1spjt1227548950.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/freestat/rcomp/tmp/2spxv1227548950.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/freestat/rcomp/tmp/32lki1227548951.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/freestat/rcomp/tmp/4m6ta1227548951.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/freestat/rcomp/tmp/5arsw1227548951.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 = 79
Frequency = 1
1 2 3 4 5 6
-0.86774194 -1.06774194 -1.26774194 -1.36774194 -0.66774194 -0.56774194
7 8 9 10 11 12
-0.36774194 -0.16774194 -0.06774194 0.03225806 0.13225806 0.03225806
13 14 15 16 17 18
-0.16774194 -0.26774194 -0.56774194 -0.56774194 0.03225806 0.13225806
19 20 21 22 23 24
0.23225806 0.23225806 0.23225806 0.43225806 0.73225806 0.63225806
25 26 27 28 29 30
0.13225806 -0.66774194 -0.96774194 -0.66774194 0.63225806 1.13225806
31 32 33 34 35 36
1.13225806 0.53225806 0.03225806 0.13225806 0.33225806 0.43225806
37 38 39 40 41 42
0.43225806 0.03225806 -0.06774194 -0.16774194 0.43225806 0.53225806
43 44 45 46 47 48
0.63225806 0.33225806 0.23225806 0.33225806 0.53225806 0.53225806
49 50 51 52 53 54
0.43225806 0.33225806 0.13225806 -0.06774194 0.03225806 -0.06774194
55 56 57 58 59 60
-0.06774194 -0.26774194 -0.26774194 -0.26774194 -0.16774194 -0.16774194
61 62 63 64 65 66
-0.26774194 -0.16774194 0.72352941 0.22352941 0.52352941 0.32352941
67 68 69 70 71 72
0.02352941 0.02352941 0.02352941 0.22352941 0.32352941 0.12352941
73 74 75 76 77 78
-0.17647059 -0.37647059 -0.77647059 -0.77647059 -0.17647059 -0.07647059
79
-0.17647059
> postscript(file="/var/www/html/freestat/rcomp/tmp/6yczw1227548951.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 = 79
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.86774194 NA
1 -1.06774194 -0.86774194
2 -1.26774194 -1.06774194
3 -1.36774194 -1.26774194
4 -0.66774194 -1.36774194
5 -0.56774194 -0.66774194
6 -0.36774194 -0.56774194
7 -0.16774194 -0.36774194
8 -0.06774194 -0.16774194
9 0.03225806 -0.06774194
10 0.13225806 0.03225806
11 0.03225806 0.13225806
12 -0.16774194 0.03225806
13 -0.26774194 -0.16774194
14 -0.56774194 -0.26774194
15 -0.56774194 -0.56774194
16 0.03225806 -0.56774194
17 0.13225806 0.03225806
18 0.23225806 0.13225806
19 0.23225806 0.23225806
20 0.23225806 0.23225806
21 0.43225806 0.23225806
22 0.73225806 0.43225806
23 0.63225806 0.73225806
24 0.13225806 0.63225806
25 -0.66774194 0.13225806
26 -0.96774194 -0.66774194
27 -0.66774194 -0.96774194
28 0.63225806 -0.66774194
29 1.13225806 0.63225806
30 1.13225806 1.13225806
31 0.53225806 1.13225806
32 0.03225806 0.53225806
33 0.13225806 0.03225806
34 0.33225806 0.13225806
35 0.43225806 0.33225806
36 0.43225806 0.43225806
37 0.03225806 0.43225806
38 -0.06774194 0.03225806
39 -0.16774194 -0.06774194
40 0.43225806 -0.16774194
41 0.53225806 0.43225806
42 0.63225806 0.53225806
43 0.33225806 0.63225806
44 0.23225806 0.33225806
45 0.33225806 0.23225806
46 0.53225806 0.33225806
47 0.53225806 0.53225806
48 0.43225806 0.53225806
49 0.33225806 0.43225806
50 0.13225806 0.33225806
51 -0.06774194 0.13225806
52 0.03225806 -0.06774194
53 -0.06774194 0.03225806
54 -0.06774194 -0.06774194
55 -0.26774194 -0.06774194
56 -0.26774194 -0.26774194
57 -0.26774194 -0.26774194
58 -0.16774194 -0.26774194
59 -0.16774194 -0.16774194
60 -0.26774194 -0.16774194
61 -0.16774194 -0.26774194
62 0.72352941 -0.16774194
63 0.22352941 0.72352941
64 0.52352941 0.22352941
65 0.32352941 0.52352941
66 0.02352941 0.32352941
67 0.02352941 0.02352941
68 0.02352941 0.02352941
69 0.22352941 0.02352941
70 0.32352941 0.22352941
71 0.12352941 0.32352941
72 -0.17647059 0.12352941
73 -0.37647059 -0.17647059
74 -0.77647059 -0.37647059
75 -0.77647059 -0.77647059
76 -0.17647059 -0.77647059
77 -0.07647059 -0.17647059
78 -0.17647059 -0.07647059
79 NA -0.17647059
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.06774194 -0.86774194
[2,] -1.26774194 -1.06774194
[3,] -1.36774194 -1.26774194
[4,] -0.66774194 -1.36774194
[5,] -0.56774194 -0.66774194
[6,] -0.36774194 -0.56774194
[7,] -0.16774194 -0.36774194
[8,] -0.06774194 -0.16774194
[9,] 0.03225806 -0.06774194
[10,] 0.13225806 0.03225806
[11,] 0.03225806 0.13225806
[12,] -0.16774194 0.03225806
[13,] -0.26774194 -0.16774194
[14,] -0.56774194 -0.26774194
[15,] -0.56774194 -0.56774194
[16,] 0.03225806 -0.56774194
[17,] 0.13225806 0.03225806
[18,] 0.23225806 0.13225806
[19,] 0.23225806 0.23225806
[20,] 0.23225806 0.23225806
[21,] 0.43225806 0.23225806
[22,] 0.73225806 0.43225806
[23,] 0.63225806 0.73225806
[24,] 0.13225806 0.63225806
[25,] -0.66774194 0.13225806
[26,] -0.96774194 -0.66774194
[27,] -0.66774194 -0.96774194
[28,] 0.63225806 -0.66774194
[29,] 1.13225806 0.63225806
[30,] 1.13225806 1.13225806
[31,] 0.53225806 1.13225806
[32,] 0.03225806 0.53225806
[33,] 0.13225806 0.03225806
[34,] 0.33225806 0.13225806
[35,] 0.43225806 0.33225806
[36,] 0.43225806 0.43225806
[37,] 0.03225806 0.43225806
[38,] -0.06774194 0.03225806
[39,] -0.16774194 -0.06774194
[40,] 0.43225806 -0.16774194
[41,] 0.53225806 0.43225806
[42,] 0.63225806 0.53225806
[43,] 0.33225806 0.63225806
[44,] 0.23225806 0.33225806
[45,] 0.33225806 0.23225806
[46,] 0.53225806 0.33225806
[47,] 0.53225806 0.53225806
[48,] 0.43225806 0.53225806
[49,] 0.33225806 0.43225806
[50,] 0.13225806 0.33225806
[51,] -0.06774194 0.13225806
[52,] 0.03225806 -0.06774194
[53,] -0.06774194 0.03225806
[54,] -0.06774194 -0.06774194
[55,] -0.26774194 -0.06774194
[56,] -0.26774194 -0.26774194
[57,] -0.26774194 -0.26774194
[58,] -0.16774194 -0.26774194
[59,] -0.16774194 -0.16774194
[60,] -0.26774194 -0.16774194
[61,] -0.16774194 -0.26774194
[62,] 0.72352941 -0.16774194
[63,] 0.22352941 0.72352941
[64,] 0.52352941 0.22352941
[65,] 0.32352941 0.52352941
[66,] 0.02352941 0.32352941
[67,] 0.02352941 0.02352941
[68,] 0.02352941 0.02352941
[69,] 0.22352941 0.02352941
[70,] 0.32352941 0.22352941
[71,] 0.12352941 0.32352941
[72,] -0.17647059 0.12352941
[73,] -0.37647059 -0.17647059
[74,] -0.77647059 -0.37647059
[75,] -0.77647059 -0.77647059
[76,] -0.17647059 -0.77647059
[77,] -0.07647059 -0.17647059
[78,] -0.17647059 -0.07647059
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.06774194 -0.86774194
2 -1.26774194 -1.06774194
3 -1.36774194 -1.26774194
4 -0.66774194 -1.36774194
5 -0.56774194 -0.66774194
6 -0.36774194 -0.56774194
7 -0.16774194 -0.36774194
8 -0.06774194 -0.16774194
9 0.03225806 -0.06774194
10 0.13225806 0.03225806
11 0.03225806 0.13225806
12 -0.16774194 0.03225806
13 -0.26774194 -0.16774194
14 -0.56774194 -0.26774194
15 -0.56774194 -0.56774194
16 0.03225806 -0.56774194
17 0.13225806 0.03225806
18 0.23225806 0.13225806
19 0.23225806 0.23225806
20 0.23225806 0.23225806
21 0.43225806 0.23225806
22 0.73225806 0.43225806
23 0.63225806 0.73225806
24 0.13225806 0.63225806
25 -0.66774194 0.13225806
26 -0.96774194 -0.66774194
27 -0.66774194 -0.96774194
28 0.63225806 -0.66774194
29 1.13225806 0.63225806
30 1.13225806 1.13225806
31 0.53225806 1.13225806
32 0.03225806 0.53225806
33 0.13225806 0.03225806
34 0.33225806 0.13225806
35 0.43225806 0.33225806
36 0.43225806 0.43225806
37 0.03225806 0.43225806
38 -0.06774194 0.03225806
39 -0.16774194 -0.06774194
40 0.43225806 -0.16774194
41 0.53225806 0.43225806
42 0.63225806 0.53225806
43 0.33225806 0.63225806
44 0.23225806 0.33225806
45 0.33225806 0.23225806
46 0.53225806 0.33225806
47 0.53225806 0.53225806
48 0.43225806 0.53225806
49 0.33225806 0.43225806
50 0.13225806 0.33225806
51 -0.06774194 0.13225806
52 0.03225806 -0.06774194
53 -0.06774194 0.03225806
54 -0.06774194 -0.06774194
55 -0.26774194 -0.06774194
56 -0.26774194 -0.26774194
57 -0.26774194 -0.26774194
58 -0.16774194 -0.26774194
59 -0.16774194 -0.16774194
60 -0.26774194 -0.16774194
61 -0.16774194 -0.26774194
62 0.72352941 -0.16774194
63 0.22352941 0.72352941
64 0.52352941 0.22352941
65 0.32352941 0.52352941
66 0.02352941 0.32352941
67 0.02352941 0.02352941
68 0.02352941 0.02352941
69 0.22352941 0.02352941
70 0.32352941 0.22352941
71 0.12352941 0.32352941
72 -0.17647059 0.12352941
73 -0.37647059 -0.17647059
74 -0.77647059 -0.37647059
75 -0.77647059 -0.77647059
76 -0.17647059 -0.77647059
77 -0.07647059 -0.17647059
78 -0.17647059 -0.07647059
> 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/freestat/rcomp/tmp/765i71227548951.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/freestat/rcomp/tmp/82y651227548951.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/freestat/rcomp/tmp/9l7qi1227548951.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/freestat/rcomp/tmp/10mfan1227548951.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11dvtv1227548951.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/freestat/rcomp/tmp/12ef9j1227548951.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/freestat/rcomp/tmp/13lb311227548951.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/freestat/rcomp/tmp/14q0s21227548951.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/freestat/rcomp/tmp/15wumr1227548951.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/freestat/rcomp/tmp/16gy671227548951.tab")
+ }
>
> system("convert tmp/1spjt1227548950.ps tmp/1spjt1227548950.png")
> system("convert tmp/2spxv1227548950.ps tmp/2spxv1227548950.png")
> system("convert tmp/32lki1227548951.ps tmp/32lki1227548951.png")
> system("convert tmp/4m6ta1227548951.ps tmp/4m6ta1227548951.png")
> system("convert tmp/5arsw1227548951.ps tmp/5arsw1227548951.png")
> system("convert tmp/6yczw1227548951.ps tmp/6yczw1227548951.png")
> system("convert tmp/765i71227548951.ps tmp/765i71227548951.png")
> system("convert tmp/82y651227548951.ps tmp/82y651227548951.png")
> system("convert tmp/9l7qi1227548951.ps tmp/9l7qi1227548951.png")
> system("convert tmp/10mfan1227548951.ps tmp/10mfan1227548951.png")
>
>
> proc.time()
user system elapsed
3.921 2.513 4.442