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(98.71,153.4,98.54,145,98.2,137.7,96.92,148.3,99.06,152.2,99.65,169.4,99.82,168.6,99.99,161.1,100.33,174.1,99.31,179,101.1,190.6,101.1,190,100.93,181.6,100.85,174.8,100.93,180.5,99.6,196.8,101.88,193.8,101.81,197,102.38,216.3,102.74,221.4,102.82,217.9,101.72,229.7,103.47,227.4,102.98,204.2,102.68,196.6,102.9,198.8,103.03,207.5,101.29,190.7,103.69,201.6,103.68,210.5,104.2,223.5,104.08,223.8,104.16,231.2,103.05,244,104.66,234.7,104.46,250.2,104.95,265.7,105.85,287.6,106.23,283.3,104.86,295.4,107.44,312.3,108.23,333.8,108.45,347.7,109.39,383.2,110.15,407.1,109.13,413.6,110.28,362.7,110.17,321.9,109.99,239.4,109.26,191,109.11,159.7,107.06,163.4,109.53,157.6,108.92,166.2,109.24,176.7,109.12,198.3,109,226.2,107.23,216.2,109.49,235.9,109.04,226.9),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),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 = '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 t
1 98.71 153.4 1 0 0 0 0 0 0 0 0 0 0 1
2 98.54 145.0 0 1 0 0 0 0 0 0 0 0 0 2
3 98.20 137.7 0 0 1 0 0 0 0 0 0 0 0 3
4 96.92 148.3 0 0 0 1 0 0 0 0 0 0 0 4
5 99.06 152.2 0 0 0 0 1 0 0 0 0 0 0 5
6 99.65 169.4 0 0 0 0 0 1 0 0 0 0 0 6
7 99.82 168.6 0 0 0 0 0 0 1 0 0 0 0 7
8 99.99 161.1 0 0 0 0 0 0 0 1 0 0 0 8
9 100.33 174.1 0 0 0 0 0 0 0 0 1 0 0 9
10 99.31 179.0 0 0 0 0 0 0 0 0 0 1 0 10
11 101.10 190.6 0 0 0 0 0 0 0 0 0 0 1 11
12 101.10 190.0 0 0 0 0 0 0 0 0 0 0 0 12
13 100.93 181.6 1 0 0 0 0 0 0 0 0 0 0 13
14 100.85 174.8 0 1 0 0 0 0 0 0 0 0 0 14
15 100.93 180.5 0 0 1 0 0 0 0 0 0 0 0 15
16 99.60 196.8 0 0 0 1 0 0 0 0 0 0 0 16
17 101.88 193.8 0 0 0 0 1 0 0 0 0 0 0 17
18 101.81 197.0 0 0 0 0 0 1 0 0 0 0 0 18
19 102.38 216.3 0 0 0 0 0 0 1 0 0 0 0 19
20 102.74 221.4 0 0 0 0 0 0 0 1 0 0 0 20
21 102.82 217.9 0 0 0 0 0 0 0 0 1 0 0 21
22 101.72 229.7 0 0 0 0 0 0 0 0 0 1 0 22
23 103.47 227.4 0 0 0 0 0 0 0 0 0 0 1 23
24 102.98 204.2 0 0 0 0 0 0 0 0 0 0 0 24
25 102.68 196.6 1 0 0 0 0 0 0 0 0 0 0 25
26 102.90 198.8 0 1 0 0 0 0 0 0 0 0 0 26
27 103.03 207.5 0 0 1 0 0 0 0 0 0 0 0 27
28 101.29 190.7 0 0 0 1 0 0 0 0 0 0 0 28
29 103.69 201.6 0 0 0 0 1 0 0 0 0 0 0 29
30 103.68 210.5 0 0 0 0 0 1 0 0 0 0 0 30
31 104.20 223.5 0 0 0 0 0 0 1 0 0 0 0 31
32 104.08 223.8 0 0 0 0 0 0 0 1 0 0 0 32
33 104.16 231.2 0 0 0 0 0 0 0 0 1 0 0 33
34 103.05 244.0 0 0 0 0 0 0 0 0 0 1 0 34
35 104.66 234.7 0 0 0 0 0 0 0 0 0 0 1 35
36 104.46 250.2 0 0 0 0 0 0 0 0 0 0 0 36
37 104.95 265.7 1 0 0 0 0 0 0 0 0 0 0 37
38 105.85 287.6 0 1 0 0 0 0 0 0 0 0 0 38
39 106.23 283.3 0 0 1 0 0 0 0 0 0 0 0 39
40 104.86 295.4 0 0 0 1 0 0 0 0 0 0 0 40
41 107.44 312.3 0 0 0 0 1 0 0 0 0 0 0 41
42 108.23 333.8 0 0 0 0 0 1 0 0 0 0 0 42
43 108.45 347.7 0 0 0 0 0 0 1 0 0 0 0 43
44 109.39 383.2 0 0 0 0 0 0 0 1 0 0 0 44
45 110.15 407.1 0 0 0 0 0 0 0 0 1 0 0 45
46 109.13 413.6 0 0 0 0 0 0 0 0 0 1 0 46
47 110.28 362.7 0 0 0 0 0 0 0 0 0 0 1 47
48 110.17 321.9 0 0 0 0 0 0 0 0 0 0 0 48
49 109.99 239.4 1 0 0 0 0 0 0 0 0 0 0 49
50 109.26 191.0 0 1 0 0 0 0 0 0 0 0 0 50
51 109.11 159.7 0 0 1 0 0 0 0 0 0 0 0 51
52 107.06 163.4 0 0 0 1 0 0 0 0 0 0 0 52
53 109.53 157.6 0 0 0 0 1 0 0 0 0 0 0 53
54 108.92 166.2 0 0 0 0 0 1 0 0 0 0 0 54
55 109.24 176.7 0 0 0 0 0 0 1 0 0 0 0 55
56 109.12 198.3 0 0 0 0 0 0 0 1 0 0 0 56
57 109.00 226.2 0 0 0 0 0 0 0 0 1 0 0 57
58 107.23 216.2 0 0 0 0 0 0 0 0 0 1 0 58
59 109.49 235.9 0 0 0 0 0 0 0 0 0 0 1 59
60 109.04 226.9 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
95.98827 0.01125 0.35562 0.28140 0.17445 -1.62886
M5 M6 M7 M8 M9 M10
0.50257 0.31592 0.35914 0.29038 0.17280 -1.28073
M11 t
0.31038 0.19106
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.45099 -0.58647 -0.02623 0.40541 1.59208
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 95.988271 0.540531 177.581 < 2e-16 ***
X 0.011246 0.001895 5.935 3.63e-07 ***
M1 0.355621 0.510941 0.696 0.48992
M2 0.281404 0.511393 0.550 0.58480
M3 0.174447 0.512072 0.341 0.73490
M4 -1.628863 0.510748 -3.189 0.00257 **
M5 0.502574 0.509681 0.986 0.32926
M6 0.315919 0.507863 0.622 0.53698
M7 0.359135 0.506874 0.709 0.48219
M8 0.290376 0.506595 0.573 0.56931
M9 0.172803 0.507326 0.341 0.73494
M10 -1.280732 0.507632 -2.523 0.01516 *
M11 0.310383 0.506564 0.613 0.54308
t 0.191058 0.006825 27.993 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7999 on 46 degrees of freedom
Multiple R-squared: 0.9664, Adjusted R-squared: 0.9569
F-statistic: 101.8 on 13 and 46 DF, p-value: < 2.2e-16
> 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,] 1.856889e-03 3.713777e-03 0.9981431
[2,] 2.027335e-04 4.054671e-04 0.9997973
[3,] 6.057359e-05 1.211472e-04 0.9999394
[4,] 1.579291e-05 3.158582e-05 0.9999842
[5,] 2.272127e-06 4.544253e-06 0.9999977
[6,] 1.211360e-06 2.422721e-06 0.9999988
[7,] 2.132665e-07 4.265331e-07 0.9999998
[8,] 5.890061e-08 1.178012e-07 0.9999999
[9,] 1.675160e-08 3.350319e-08 1.0000000
[10,] 1.991934e-09 3.983868e-09 1.0000000
[11,] 2.347649e-10 4.695297e-10 1.0000000
[12,] 8.408964e-11 1.681793e-10 1.0000000
[13,] 1.492335e-11 2.984670e-11 1.0000000
[14,] 2.135575e-12 4.271149e-12 1.0000000
[15,] 2.240781e-13 4.481562e-13 1.0000000
[16,] 1.088558e-12 2.177116e-12 1.0000000
[17,] 1.214370e-11 2.428739e-11 1.0000000
[18,] 1.191288e-10 2.382576e-10 1.0000000
[19,] 1.534099e-10 3.068198e-10 1.0000000
[20,] 1.272314e-08 2.544627e-08 1.0000000
[21,] 1.158424e-08 2.316848e-08 1.0000000
[22,] 1.030339e-08 2.060678e-08 1.0000000
[23,] 1.000640e-07 2.001280e-07 0.9999999
[24,] 7.917842e-07 1.583568e-06 0.9999992
[25,] 1.022863e-04 2.045727e-04 0.9998977
[26,] 6.621152e-04 1.324230e-03 0.9993379
[27,] 4.781810e-02 9.563620e-02 0.9521819
> postscript(file="/var/www/html/rcomp/tmp/1irqg1258720136.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/2mqt01258720136.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/3mh9x1258720136.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/4d6s61258720136.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/50wwc1258720136.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
0.44997076 0.25759357 -0.08441343 0.12863527 -0.09771770 0.29445556
7 8 9 10 11 12
0.23917818 0.37122229 0.49154405 0.67891774 0.55629575 0.68236893
13 14 15 16 17 18
0.06015326 -0.06021693 -0.12841709 -0.02946847 -0.03822661 -0.14861456
19 20 21 22 23 24
-0.02992906 0.15042014 0.19629476 0.22607362 0.21976585 0.10999022
25 26 27 28 29 30
-0.65122196 -0.57280280 -0.62473984 -0.56356094 -0.60863331 -0.72312134
31 32 33 34 35 36
-0.58358838 -0.82926016 -0.90596289 -0.89742966 -0.96501804 -1.21999945
37 38 39 40 41 42
-1.45098563 -0.91410533 -0.56984921 -0.46366895 -0.39621508 0.14760198
43 44 45 46 47 48
-0.02298613 0.39549600 0.81324042 0.98262111 0.92285083 1.39099825
49 50 51 52 53 54
1.59208358 1.28953149 1.40741956 0.92806309 1.14079270 0.42967836
55 56 57 58 59 60
0.39732539 -0.08787827 -0.59511635 -0.99018281 -0.73389439 -0.96335794
> postscript(file="/var/www/html/rcomp/tmp/6qodt1258720136.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 0.44997076 NA
1 0.25759357 0.44997076
2 -0.08441343 0.25759357
3 0.12863527 -0.08441343
4 -0.09771770 0.12863527
5 0.29445556 -0.09771770
6 0.23917818 0.29445556
7 0.37122229 0.23917818
8 0.49154405 0.37122229
9 0.67891774 0.49154405
10 0.55629575 0.67891774
11 0.68236893 0.55629575
12 0.06015326 0.68236893
13 -0.06021693 0.06015326
14 -0.12841709 -0.06021693
15 -0.02946847 -0.12841709
16 -0.03822661 -0.02946847
17 -0.14861456 -0.03822661
18 -0.02992906 -0.14861456
19 0.15042014 -0.02992906
20 0.19629476 0.15042014
21 0.22607362 0.19629476
22 0.21976585 0.22607362
23 0.10999022 0.21976585
24 -0.65122196 0.10999022
25 -0.57280280 -0.65122196
26 -0.62473984 -0.57280280
27 -0.56356094 -0.62473984
28 -0.60863331 -0.56356094
29 -0.72312134 -0.60863331
30 -0.58358838 -0.72312134
31 -0.82926016 -0.58358838
32 -0.90596289 -0.82926016
33 -0.89742966 -0.90596289
34 -0.96501804 -0.89742966
35 -1.21999945 -0.96501804
36 -1.45098563 -1.21999945
37 -0.91410533 -1.45098563
38 -0.56984921 -0.91410533
39 -0.46366895 -0.56984921
40 -0.39621508 -0.46366895
41 0.14760198 -0.39621508
42 -0.02298613 0.14760198
43 0.39549600 -0.02298613
44 0.81324042 0.39549600
45 0.98262111 0.81324042
46 0.92285083 0.98262111
47 1.39099825 0.92285083
48 1.59208358 1.39099825
49 1.28953149 1.59208358
50 1.40741956 1.28953149
51 0.92806309 1.40741956
52 1.14079270 0.92806309
53 0.42967836 1.14079270
54 0.39732539 0.42967836
55 -0.08787827 0.39732539
56 -0.59511635 -0.08787827
57 -0.99018281 -0.59511635
58 -0.73389439 -0.99018281
59 -0.96335794 -0.73389439
60 NA -0.96335794
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.25759357 0.44997076
[2,] -0.08441343 0.25759357
[3,] 0.12863527 -0.08441343
[4,] -0.09771770 0.12863527
[5,] 0.29445556 -0.09771770
[6,] 0.23917818 0.29445556
[7,] 0.37122229 0.23917818
[8,] 0.49154405 0.37122229
[9,] 0.67891774 0.49154405
[10,] 0.55629575 0.67891774
[11,] 0.68236893 0.55629575
[12,] 0.06015326 0.68236893
[13,] -0.06021693 0.06015326
[14,] -0.12841709 -0.06021693
[15,] -0.02946847 -0.12841709
[16,] -0.03822661 -0.02946847
[17,] -0.14861456 -0.03822661
[18,] -0.02992906 -0.14861456
[19,] 0.15042014 -0.02992906
[20,] 0.19629476 0.15042014
[21,] 0.22607362 0.19629476
[22,] 0.21976585 0.22607362
[23,] 0.10999022 0.21976585
[24,] -0.65122196 0.10999022
[25,] -0.57280280 -0.65122196
[26,] -0.62473984 -0.57280280
[27,] -0.56356094 -0.62473984
[28,] -0.60863331 -0.56356094
[29,] -0.72312134 -0.60863331
[30,] -0.58358838 -0.72312134
[31,] -0.82926016 -0.58358838
[32,] -0.90596289 -0.82926016
[33,] -0.89742966 -0.90596289
[34,] -0.96501804 -0.89742966
[35,] -1.21999945 -0.96501804
[36,] -1.45098563 -1.21999945
[37,] -0.91410533 -1.45098563
[38,] -0.56984921 -0.91410533
[39,] -0.46366895 -0.56984921
[40,] -0.39621508 -0.46366895
[41,] 0.14760198 -0.39621508
[42,] -0.02298613 0.14760198
[43,] 0.39549600 -0.02298613
[44,] 0.81324042 0.39549600
[45,] 0.98262111 0.81324042
[46,] 0.92285083 0.98262111
[47,] 1.39099825 0.92285083
[48,] 1.59208358 1.39099825
[49,] 1.28953149 1.59208358
[50,] 1.40741956 1.28953149
[51,] 0.92806309 1.40741956
[52,] 1.14079270 0.92806309
[53,] 0.42967836 1.14079270
[54,] 0.39732539 0.42967836
[55,] -0.08787827 0.39732539
[56,] -0.59511635 -0.08787827
[57,] -0.99018281 -0.59511635
[58,] -0.73389439 -0.99018281
[59,] -0.96335794 -0.73389439
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.25759357 0.44997076
2 -0.08441343 0.25759357
3 0.12863527 -0.08441343
4 -0.09771770 0.12863527
5 0.29445556 -0.09771770
6 0.23917818 0.29445556
7 0.37122229 0.23917818
8 0.49154405 0.37122229
9 0.67891774 0.49154405
10 0.55629575 0.67891774
11 0.68236893 0.55629575
12 0.06015326 0.68236893
13 -0.06021693 0.06015326
14 -0.12841709 -0.06021693
15 -0.02946847 -0.12841709
16 -0.03822661 -0.02946847
17 -0.14861456 -0.03822661
18 -0.02992906 -0.14861456
19 0.15042014 -0.02992906
20 0.19629476 0.15042014
21 0.22607362 0.19629476
22 0.21976585 0.22607362
23 0.10999022 0.21976585
24 -0.65122196 0.10999022
25 -0.57280280 -0.65122196
26 -0.62473984 -0.57280280
27 -0.56356094 -0.62473984
28 -0.60863331 -0.56356094
29 -0.72312134 -0.60863331
30 -0.58358838 -0.72312134
31 -0.82926016 -0.58358838
32 -0.90596289 -0.82926016
33 -0.89742966 -0.90596289
34 -0.96501804 -0.89742966
35 -1.21999945 -0.96501804
36 -1.45098563 -1.21999945
37 -0.91410533 -1.45098563
38 -0.56984921 -0.91410533
39 -0.46366895 -0.56984921
40 -0.39621508 -0.46366895
41 0.14760198 -0.39621508
42 -0.02298613 0.14760198
43 0.39549600 -0.02298613
44 0.81324042 0.39549600
45 0.98262111 0.81324042
46 0.92285083 0.98262111
47 1.39099825 0.92285083
48 1.59208358 1.39099825
49 1.28953149 1.59208358
50 1.40741956 1.28953149
51 0.92806309 1.40741956
52 1.14079270 0.92806309
53 0.42967836 1.14079270
54 0.39732539 0.42967836
55 -0.08787827 0.39732539
56 -0.59511635 -0.08787827
57 -0.99018281 -0.59511635
58 -0.73389439 -0.99018281
59 -0.96335794 -0.73389439
> 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/76fo61258720136.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/8x2k91258720136.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/961uz1258720136.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/103eq01258720136.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/11kfb11258720136.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/12ccl71258720137.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/13gcfl1258720137.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/143ln61258720137.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/15fn6h1258720137.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/16bcqw1258720137.tab")
+ }
>
> system("convert tmp/1irqg1258720136.ps tmp/1irqg1258720136.png")
> system("convert tmp/2mqt01258720136.ps tmp/2mqt01258720136.png")
> system("convert tmp/3mh9x1258720136.ps tmp/3mh9x1258720136.png")
> system("convert tmp/4d6s61258720136.ps tmp/4d6s61258720136.png")
> system("convert tmp/50wwc1258720136.ps tmp/50wwc1258720136.png")
> system("convert tmp/6qodt1258720136.ps tmp/6qodt1258720136.png")
> system("convert tmp/76fo61258720136.ps tmp/76fo61258720136.png")
> system("convert tmp/8x2k91258720136.ps tmp/8x2k91258720136.png")
> system("convert tmp/961uz1258720136.ps tmp/961uz1258720136.png")
> system("convert tmp/103eq01258720136.ps tmp/103eq01258720136.png")
>
>
> proc.time()
user system elapsed
2.379 1.572 2.753