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(121.6,0,118.8,0,114.0,1,111.5,1,97.2,1,102.5,1,113.4,1,109.8,1,104.9,1,126.1,1,80.0,1,96.8,1,117.2,1,112.3,1,117.3,1,111.1,0,102.2,0,104.3,0,122.9,0,107.6,0,121.3,0,131.5,0,89.0,0,104.4,0,128.9,0,135.9,0,133.3,0,121.3,0,120.5,0,120.4,0,137.9,0,126.1,0,133.2,0,151.1,0,105.0,0,119.0,0,140.4,0,156.6,0,137.1,0,122.7,0,125.8,0,139.3,0,134.9,0,149.2,1,132.3,0,149.0,1,117.2,1,119.6,1,152.0,1,149.4,1,127.3,1,114.1,1,102.1,1,107.7,1,104.4,1,102.1,1,96.0,1,109.3,1,90.0,1,83.9,1),dim=c(2,60),dimnames=list(c('Promet','Dummy'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Promet','Dummy'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = '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
Promet Dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 121.6 0 1 0 0 0 0 0 0 0 0 0 0
2 118.8 0 0 1 0 0 0 0 0 0 0 0 0
3 114.0 1 0 0 1 0 0 0 0 0 0 0 0
4 111.5 1 0 0 0 1 0 0 0 0 0 0 0
5 97.2 1 0 0 0 0 1 0 0 0 0 0 0
6 102.5 1 0 0 0 0 0 1 0 0 0 0 0
7 113.4 1 0 0 0 0 0 0 1 0 0 0 0
8 109.8 1 0 0 0 0 0 0 0 1 0 0 0
9 104.9 1 0 0 0 0 0 0 0 0 1 0 0
10 126.1 1 0 0 0 0 0 0 0 0 0 1 0
11 80.0 1 0 0 0 0 0 0 0 0 0 0 1
12 96.8 1 0 0 0 0 0 0 0 0 0 0 0
13 117.2 1 1 0 0 0 0 0 0 0 0 0 0
14 112.3 1 0 1 0 0 0 0 0 0 0 0 0
15 117.3 1 0 0 1 0 0 0 0 0 0 0 0
16 111.1 0 0 0 0 1 0 0 0 0 0 0 0
17 102.2 0 0 0 0 0 1 0 0 0 0 0 0
18 104.3 0 0 0 0 0 0 1 0 0 0 0 0
19 122.9 0 0 0 0 0 0 0 1 0 0 0 0
20 107.6 0 0 0 0 0 0 0 0 1 0 0 0
21 121.3 0 0 0 0 0 0 0 0 0 1 0 0
22 131.5 0 0 0 0 0 0 0 0 0 0 1 0
23 89.0 0 0 0 0 0 0 0 0 0 0 0 1
24 104.4 0 0 0 0 0 0 0 0 0 0 0 0
25 128.9 0 1 0 0 0 0 0 0 0 0 0 0
26 135.9 0 0 1 0 0 0 0 0 0 0 0 0
27 133.3 0 0 0 1 0 0 0 0 0 0 0 0
28 121.3 0 0 0 0 1 0 0 0 0 0 0 0
29 120.5 0 0 0 0 0 1 0 0 0 0 0 0
30 120.4 0 0 0 0 0 0 1 0 0 0 0 0
31 137.9 0 0 0 0 0 0 0 1 0 0 0 0
32 126.1 0 0 0 0 0 0 0 0 1 0 0 0
33 133.2 0 0 0 0 0 0 0 0 0 1 0 0
34 151.1 0 0 0 0 0 0 0 0 0 0 1 0
35 105.0 0 0 0 0 0 0 0 0 0 0 0 1
36 119.0 0 0 0 0 0 0 0 0 0 0 0 0
37 140.4 0 1 0 0 0 0 0 0 0 0 0 0
38 156.6 0 0 1 0 0 0 0 0 0 0 0 0
39 137.1 0 0 0 1 0 0 0 0 0 0 0 0
40 122.7 0 0 0 0 1 0 0 0 0 0 0 0
41 125.8 0 0 0 0 0 1 0 0 0 0 0 0
42 139.3 0 0 0 0 0 0 1 0 0 0 0 0
43 134.9 0 0 0 0 0 0 0 1 0 0 0 0
44 149.2 1 0 0 0 0 0 0 0 1 0 0 0
45 132.3 0 0 0 0 0 0 0 0 0 1 0 0
46 149.0 1 0 0 0 0 0 0 0 0 0 1 0
47 117.2 1 0 0 0 0 0 0 0 0 0 0 1
48 119.6 1 0 0 0 0 0 0 0 0 0 0 0
49 152.0 1 1 0 0 0 0 0 0 0 0 0 0
50 149.4 1 0 1 0 0 0 0 0 0 0 0 0
51 127.3 1 0 0 1 0 0 0 0 0 0 0 0
52 114.1 1 0 0 0 1 0 0 0 0 0 0 0
53 102.1 1 0 0 0 0 1 0 0 0 0 0 0
54 107.7 1 0 0 0 0 0 1 0 0 0 0 0
55 104.4 1 0 0 0 0 0 0 1 0 0 0 0
56 102.1 1 0 0 0 0 0 0 0 1 0 0 0
57 96.0 1 0 0 0 0 0 0 0 0 1 0 0
58 109.3 1 0 0 0 0 0 0 0 0 0 1 0
59 90.0 1 0 0 0 0 0 0 0 0 0 0 1
60 83.9 1 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dummy M1 M2 M3 M4
111.237 -10.828 25.114 27.694 21.060 9.234
M5 M6 M7 M8 M9 M10
2.654 7.934 15.794 14.220 10.634 28.660
M11
-8.500
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-20.131 -8.342 -1.436 6.897 34.571
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 111.237 6.575 16.918 < 2e-16 ***
Dummy -10.828 3.653 -2.964 0.00475 **
M1 25.114 8.797 2.855 0.00639 **
M2 27.694 8.797 3.148 0.00285 **
M3 21.060 8.767 2.402 0.02030 *
M4 9.234 8.797 1.050 0.29923
M5 2.654 8.797 0.302 0.76418
M6 7.934 8.797 0.902 0.37170
M7 15.794 8.797 1.795 0.07902 .
M8 14.220 8.767 1.622 0.11149
M9 10.634 8.797 1.209 0.23277
M10 28.660 8.767 3.269 0.00202 **
M11 -8.500 8.767 -0.970 0.33723
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.86 on 47 degrees of freedom
Multiple R-squared: 0.5068, Adjusted R-squared: 0.3808
F-statistic: 4.024 on 12 and 47 DF, p-value: 0.0002726
> 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,] 6.277305e-03 1.255461e-02 0.9937227
[2,] 9.575927e-04 1.915185e-03 0.9990424
[3,] 1.635744e-04 3.271489e-04 0.9998364
[4,] 9.956596e-05 1.991319e-04 0.9999004
[5,] 7.775686e-05 1.555137e-04 0.9999222
[6,] 3.247049e-04 6.494097e-04 0.9996753
[7,] 8.307090e-05 1.661418e-04 0.9999169
[8,] 3.243819e-05 6.487639e-05 0.9999676
[9,] 8.472809e-06 1.694562e-05 0.9999915
[10,] 7.863579e-06 1.572716e-05 0.9999921
[11,] 2.523734e-04 5.047469e-04 0.9997476
[12,] 2.285926e-04 4.571851e-04 0.9997714
[13,] 1.075289e-04 2.150577e-04 0.9998925
[14,] 2.980544e-04 5.961088e-04 0.9997019
[15,] 3.134717e-04 6.269435e-04 0.9996865
[16,] 4.293781e-04 8.587561e-04 0.9995706
[17,] 4.289086e-04 8.578172e-04 0.9995711
[18,] 4.787640e-04 9.575280e-04 0.9995212
[19,] 5.545455e-04 1.109091e-03 0.9994455
[20,] 5.724657e-04 1.144931e-03 0.9994275
[21,] 3.585305e-04 7.170609e-04 0.9996415
[22,] 6.504689e-04 1.300938e-03 0.9993495
[23,] 3.094100e-03 6.188201e-03 0.9969059
[24,] 1.783958e-03 3.567916e-03 0.9982160
[25,] 1.124273e-03 2.248546e-03 0.9988757
[26,] 7.001568e-04 1.400314e-03 0.9992998
[27,] 8.088506e-04 1.617701e-03 0.9991911
[28,] 2.828246e-04 5.656493e-04 0.9997172
[29,] 1.369424e-02 2.738848e-02 0.9863058
> postscript(file="/var/www/html/rcomp/tmp/12k2e1258620338.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/24cal1258620338.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/3yz5c1258620338.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/4158i1258620338.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/5yjq21258620338.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 60
Frequency = 1
1 2 3 4 5 6
-14.7511111 -20.1311111 -7.4688889 1.8566667 -5.8633333 -5.8433333
7 8 9 10 11 12
-2.8033333 -4.8288889 -6.1433333 -2.9688889 -11.9088889 -3.6088889
13 14 15 16 17 18
-8.3233333 -15.8033333 -4.1688889 -9.3711111 -11.6911111 -14.8711111
19 20 21 22 23 24
-4.1311111 -17.8566667 -0.5711111 -8.3966667 -13.7366667 -6.8366667
25 26 27 28 29 30
-7.4511111 -3.0311111 1.0033333 0.8288889 6.6088889 1.2288889
31 32 33 34 35 36
10.8688889 0.6433333 11.3288889 11.2033333 2.2633333 7.7633333
37 38 39 40 41 42
4.0488889 17.6688889 4.8033333 2.2288889 11.9088889 20.1288889
43 44 45 46 47 48
7.8688889 34.5711111 10.4288889 19.9311111 25.2911111 19.1911111
49 50 51 52 53 54
26.4766667 21.2966667 5.8311111 4.4566667 -0.9633333 -0.6433333
55 56 57 58 59 60
-11.8033333 -12.5288889 -15.0433333 -19.7688889 -1.9088889 -16.5088889
> postscript(file="/var/www/html/rcomp/tmp/62avc1258620338.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -14.7511111 NA
1 -20.1311111 -14.7511111
2 -7.4688889 -20.1311111
3 1.8566667 -7.4688889
4 -5.8633333 1.8566667
5 -5.8433333 -5.8633333
6 -2.8033333 -5.8433333
7 -4.8288889 -2.8033333
8 -6.1433333 -4.8288889
9 -2.9688889 -6.1433333
10 -11.9088889 -2.9688889
11 -3.6088889 -11.9088889
12 -8.3233333 -3.6088889
13 -15.8033333 -8.3233333
14 -4.1688889 -15.8033333
15 -9.3711111 -4.1688889
16 -11.6911111 -9.3711111
17 -14.8711111 -11.6911111
18 -4.1311111 -14.8711111
19 -17.8566667 -4.1311111
20 -0.5711111 -17.8566667
21 -8.3966667 -0.5711111
22 -13.7366667 -8.3966667
23 -6.8366667 -13.7366667
24 -7.4511111 -6.8366667
25 -3.0311111 -7.4511111
26 1.0033333 -3.0311111
27 0.8288889 1.0033333
28 6.6088889 0.8288889
29 1.2288889 6.6088889
30 10.8688889 1.2288889
31 0.6433333 10.8688889
32 11.3288889 0.6433333
33 11.2033333 11.3288889
34 2.2633333 11.2033333
35 7.7633333 2.2633333
36 4.0488889 7.7633333
37 17.6688889 4.0488889
38 4.8033333 17.6688889
39 2.2288889 4.8033333
40 11.9088889 2.2288889
41 20.1288889 11.9088889
42 7.8688889 20.1288889
43 34.5711111 7.8688889
44 10.4288889 34.5711111
45 19.9311111 10.4288889
46 25.2911111 19.9311111
47 19.1911111 25.2911111
48 26.4766667 19.1911111
49 21.2966667 26.4766667
50 5.8311111 21.2966667
51 4.4566667 5.8311111
52 -0.9633333 4.4566667
53 -0.6433333 -0.9633333
54 -11.8033333 -0.6433333
55 -12.5288889 -11.8033333
56 -15.0433333 -12.5288889
57 -19.7688889 -15.0433333
58 -1.9088889 -19.7688889
59 -16.5088889 -1.9088889
60 NA -16.5088889
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -20.1311111 -14.7511111
[2,] -7.4688889 -20.1311111
[3,] 1.8566667 -7.4688889
[4,] -5.8633333 1.8566667
[5,] -5.8433333 -5.8633333
[6,] -2.8033333 -5.8433333
[7,] -4.8288889 -2.8033333
[8,] -6.1433333 -4.8288889
[9,] -2.9688889 -6.1433333
[10,] -11.9088889 -2.9688889
[11,] -3.6088889 -11.9088889
[12,] -8.3233333 -3.6088889
[13,] -15.8033333 -8.3233333
[14,] -4.1688889 -15.8033333
[15,] -9.3711111 -4.1688889
[16,] -11.6911111 -9.3711111
[17,] -14.8711111 -11.6911111
[18,] -4.1311111 -14.8711111
[19,] -17.8566667 -4.1311111
[20,] -0.5711111 -17.8566667
[21,] -8.3966667 -0.5711111
[22,] -13.7366667 -8.3966667
[23,] -6.8366667 -13.7366667
[24,] -7.4511111 -6.8366667
[25,] -3.0311111 -7.4511111
[26,] 1.0033333 -3.0311111
[27,] 0.8288889 1.0033333
[28,] 6.6088889 0.8288889
[29,] 1.2288889 6.6088889
[30,] 10.8688889 1.2288889
[31,] 0.6433333 10.8688889
[32,] 11.3288889 0.6433333
[33,] 11.2033333 11.3288889
[34,] 2.2633333 11.2033333
[35,] 7.7633333 2.2633333
[36,] 4.0488889 7.7633333
[37,] 17.6688889 4.0488889
[38,] 4.8033333 17.6688889
[39,] 2.2288889 4.8033333
[40,] 11.9088889 2.2288889
[41,] 20.1288889 11.9088889
[42,] 7.8688889 20.1288889
[43,] 34.5711111 7.8688889
[44,] 10.4288889 34.5711111
[45,] 19.9311111 10.4288889
[46,] 25.2911111 19.9311111
[47,] 19.1911111 25.2911111
[48,] 26.4766667 19.1911111
[49,] 21.2966667 26.4766667
[50,] 5.8311111 21.2966667
[51,] 4.4566667 5.8311111
[52,] -0.9633333 4.4566667
[53,] -0.6433333 -0.9633333
[54,] -11.8033333 -0.6433333
[55,] -12.5288889 -11.8033333
[56,] -15.0433333 -12.5288889
[57,] -19.7688889 -15.0433333
[58,] -1.9088889 -19.7688889
[59,] -16.5088889 -1.9088889
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -20.1311111 -14.7511111
2 -7.4688889 -20.1311111
3 1.8566667 -7.4688889
4 -5.8633333 1.8566667
5 -5.8433333 -5.8633333
6 -2.8033333 -5.8433333
7 -4.8288889 -2.8033333
8 -6.1433333 -4.8288889
9 -2.9688889 -6.1433333
10 -11.9088889 -2.9688889
11 -3.6088889 -11.9088889
12 -8.3233333 -3.6088889
13 -15.8033333 -8.3233333
14 -4.1688889 -15.8033333
15 -9.3711111 -4.1688889
16 -11.6911111 -9.3711111
17 -14.8711111 -11.6911111
18 -4.1311111 -14.8711111
19 -17.8566667 -4.1311111
20 -0.5711111 -17.8566667
21 -8.3966667 -0.5711111
22 -13.7366667 -8.3966667
23 -6.8366667 -13.7366667
24 -7.4511111 -6.8366667
25 -3.0311111 -7.4511111
26 1.0033333 -3.0311111
27 0.8288889 1.0033333
28 6.6088889 0.8288889
29 1.2288889 6.6088889
30 10.8688889 1.2288889
31 0.6433333 10.8688889
32 11.3288889 0.6433333
33 11.2033333 11.3288889
34 2.2633333 11.2033333
35 7.7633333 2.2633333
36 4.0488889 7.7633333
37 17.6688889 4.0488889
38 4.8033333 17.6688889
39 2.2288889 4.8033333
40 11.9088889 2.2288889
41 20.1288889 11.9088889
42 7.8688889 20.1288889
43 34.5711111 7.8688889
44 10.4288889 34.5711111
45 19.9311111 10.4288889
46 25.2911111 19.9311111
47 19.1911111 25.2911111
48 26.4766667 19.1911111
49 21.2966667 26.4766667
50 5.8311111 21.2966667
51 4.4566667 5.8311111
52 -0.9633333 4.4566667
53 -0.6433333 -0.9633333
54 -11.8033333 -0.6433333
55 -12.5288889 -11.8033333
56 -15.0433333 -12.5288889
57 -19.7688889 -15.0433333
58 -1.9088889 -19.7688889
59 -16.5088889 -1.9088889
> 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/7ihcb1258620338.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/8idsy1258620338.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/9vo2p1258620338.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/101l7a1258620338.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/11j13s1258620338.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/12ljj41258620338.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/13ipjy1258620338.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/1412du1258620338.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/15qc5u1258620338.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/16kv8o1258620338.tab")
+ }
>
> system("convert tmp/12k2e1258620338.ps tmp/12k2e1258620338.png")
> system("convert tmp/24cal1258620338.ps tmp/24cal1258620338.png")
> system("convert tmp/3yz5c1258620338.ps tmp/3yz5c1258620338.png")
> system("convert tmp/4158i1258620338.ps tmp/4158i1258620338.png")
> system("convert tmp/5yjq21258620338.ps tmp/5yjq21258620338.png")
> system("convert tmp/62avc1258620338.ps tmp/62avc1258620338.png")
> system("convert tmp/7ihcb1258620338.ps tmp/7ihcb1258620338.png")
> system("convert tmp/8idsy1258620338.ps tmp/8idsy1258620338.png")
> system("convert tmp/9vo2p1258620338.ps tmp/9vo2p1258620338.png")
> system("convert tmp/101l7a1258620338.ps tmp/101l7a1258620338.png")
>
>
> proc.time()
user system elapsed
2.383 1.543 3.546