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(99.4,0,97.5,0,94.6,0,92.6,0,92.5,0,89.8,0,88.8,0,87.4,0,85.2,0,83.1,0,84.7,0,84.8,0,85.8,0,86.3,0,89,0,89,0,89.3,0,91.9,0,94.9,0,94.4,0,96.8,0,96.9,0,98,0,97.9,0,100.9,0,103.9,0,103.1,0,102.5,0,104.3,0,102.6,0,101.7,0,102.8,0,105.4,0,110.9,1,113.5,1,116.3,1,124,1,128.8,1,133.5,1,132.6,1,128.4,1,127.3,1,126.7,1,123.3,1,123.2,1,124.4,1,128.2,1,128.7,1,135.7,1,139,1,145.4,1,142.4,1,137.7,1,137,1,137.1,1,139.3,1,139.6,1,140.4,1,142.3,1,148.3,1),dim=c(2,60),dimnames=list(c('Grondstofprijzen','Wet'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Grondstofprijzen','Wet'),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
Grondstofprijzen Wet M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 99.4 0 1 0 0 0 0 0 0 0 0 0 0 1
2 97.5 0 0 1 0 0 0 0 0 0 0 0 0 2
3 94.6 0 0 0 1 0 0 0 0 0 0 0 0 3
4 92.6 0 0 0 0 1 0 0 0 0 0 0 0 4
5 92.5 0 0 0 0 0 1 0 0 0 0 0 0 5
6 89.8 0 0 0 0 0 0 1 0 0 0 0 0 6
7 88.8 0 0 0 0 0 0 0 1 0 0 0 0 7
8 87.4 0 0 0 0 0 0 0 0 1 0 0 0 8
9 85.2 0 0 0 0 0 0 0 0 0 1 0 0 9
10 83.1 0 0 0 0 0 0 0 0 0 0 1 0 10
11 84.7 0 0 0 0 0 0 0 0 0 0 0 1 11
12 84.8 0 0 0 0 0 0 0 0 0 0 0 0 12
13 85.8 0 1 0 0 0 0 0 0 0 0 0 0 13
14 86.3 0 0 1 0 0 0 0 0 0 0 0 0 14
15 89.0 0 0 0 1 0 0 0 0 0 0 0 0 15
16 89.0 0 0 0 0 1 0 0 0 0 0 0 0 16
17 89.3 0 0 0 0 0 1 0 0 0 0 0 0 17
18 91.9 0 0 0 0 0 0 1 0 0 0 0 0 18
19 94.9 0 0 0 0 0 0 0 1 0 0 0 0 19
20 94.4 0 0 0 0 0 0 0 0 1 0 0 0 20
21 96.8 0 0 0 0 0 0 0 0 0 1 0 0 21
22 96.9 0 0 0 0 0 0 0 0 0 0 1 0 22
23 98.0 0 0 0 0 0 0 0 0 0 0 0 1 23
24 97.9 0 0 0 0 0 0 0 0 0 0 0 0 24
25 100.9 0 1 0 0 0 0 0 0 0 0 0 0 25
26 103.9 0 0 1 0 0 0 0 0 0 0 0 0 26
27 103.1 0 0 0 1 0 0 0 0 0 0 0 0 27
28 102.5 0 0 0 0 1 0 0 0 0 0 0 0 28
29 104.3 0 0 0 0 0 1 0 0 0 0 0 0 29
30 102.6 0 0 0 0 0 0 1 0 0 0 0 0 30
31 101.7 0 0 0 0 0 0 0 1 0 0 0 0 31
32 102.8 0 0 0 0 0 0 0 0 1 0 0 0 32
33 105.4 0 0 0 0 0 0 0 0 0 1 0 0 33
34 110.9 1 0 0 0 0 0 0 0 0 0 1 0 34
35 113.5 1 0 0 0 0 0 0 0 0 0 0 1 35
36 116.3 1 0 0 0 0 0 0 0 0 0 0 0 36
37 124.0 1 1 0 0 0 0 0 0 0 0 0 0 37
38 128.8 1 0 1 0 0 0 0 0 0 0 0 0 38
39 133.5 1 0 0 1 0 0 0 0 0 0 0 0 39
40 132.6 1 0 0 0 1 0 0 0 0 0 0 0 40
41 128.4 1 0 0 0 0 1 0 0 0 0 0 0 41
42 127.3 1 0 0 0 0 0 1 0 0 0 0 0 42
43 126.7 1 0 0 0 0 0 0 1 0 0 0 0 43
44 123.3 1 0 0 0 0 0 0 0 1 0 0 0 44
45 123.2 1 0 0 0 0 0 0 0 0 1 0 0 45
46 124.4 1 0 0 0 0 0 0 0 0 0 1 0 46
47 128.2 1 0 0 0 0 0 0 0 0 0 0 1 47
48 128.7 1 0 0 0 0 0 0 0 0 0 0 0 48
49 135.7 1 1 0 0 0 0 0 0 0 0 0 0 49
50 139.0 1 0 1 0 0 0 0 0 0 0 0 0 50
51 145.4 1 0 0 1 0 0 0 0 0 0 0 0 51
52 142.4 1 0 0 0 1 0 0 0 0 0 0 0 52
53 137.7 1 0 0 0 0 1 0 0 0 0 0 0 53
54 137.0 1 0 0 0 0 0 1 0 0 0 0 0 54
55 137.1 1 0 0 0 0 0 0 1 0 0 0 0 55
56 139.3 1 0 0 0 0 0 0 0 1 0 0 0 56
57 139.6 1 0 0 0 0 0 0 0 0 1 0 0 57
58 140.4 1 0 0 0 0 0 0 0 0 0 1 0 58
59 142.3 1 0 0 0 0 0 0 0 0 0 0 1 59
60 148.3 1 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) Wet M1 M2 M3 M4
79.3417 16.0139 5.1836 6.3944 7.6853 5.6561
M5 M6 M7 M8 M9 M10
3.5469 2.0978 1.4886 0.3594 0.2303 -2.6017
M11 t
-1.1308 0.7292
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-9.644 -2.088 -0.265 2.527 14.146
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 79.34167 2.82528 28.083 < 2e-16 ***
Wet 16.01389 2.71862 5.890 4.23e-07 ***
M1 5.18361 3.29709 1.572 0.1228
M2 6.39444 3.28867 1.944 0.0580 .
M3 7.68528 3.28211 2.342 0.0236 *
M4 5.65611 3.27742 1.726 0.0911 .
M5 3.54694 3.27460 1.083 0.2844
M6 2.09778 3.27366 0.641 0.5248
M7 1.48861 3.27460 0.455 0.6515
M8 0.35944 3.27742 0.110 0.9131
M9 0.23028 3.28211 0.070 0.9444
M10 -2.60167 3.26612 -0.797 0.4298
M11 -1.13083 3.26329 -0.347 0.7305
t 0.72917 0.07848 9.291 4.00e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.158 on 46 degrees of freedom
Multiple R-squared: 0.9498, Adjusted R-squared: 0.9356
F-statistic: 66.95 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,] 0.7997383 4.005234e-01 2.002617e-01
[2,] 0.9075891 1.848217e-01 9.241086e-02
[3,] 0.9771846 4.563089e-02 2.281545e-02
[4,] 0.9910760 1.784793e-02 8.923967e-03
[5,] 0.9988034 2.393150e-03 1.196575e-03
[6,] 0.9999077 1.845413e-04 9.227067e-05
[7,] 0.9999853 2.946397e-05 1.473198e-05
[8,] 0.9999908 1.839717e-05 9.198586e-06
[9,] 0.9999773 4.539618e-05 2.269809e-05
[10,] 0.9999548 9.048264e-05 4.524132e-05
[11,] 0.9999793 4.144405e-05 2.072203e-05
[12,] 0.9999917 1.665567e-05 8.327834e-06
[13,] 0.9999786 4.274811e-05 2.137406e-05
[14,] 0.9999425 1.149212e-04 5.746060e-05
[15,] 0.9998753 2.493057e-04 1.246529e-04
[16,] 0.9997111 5.778027e-04 2.889013e-04
[17,] 0.9993098 1.380455e-03 6.902273e-04
[18,] 0.9982683 3.463354e-03 1.731677e-03
[19,] 0.9961112 7.777618e-03 3.888809e-03
[20,] 0.9929378 1.412430e-02 7.062152e-03
[21,] 0.9846591 3.068174e-02 1.534087e-02
[22,] 0.9756824 4.863524e-02 2.431762e-02
[23,] 0.9602881 7.942387e-02 3.971193e-02
[24,] 0.9470057 1.059887e-01 5.299434e-02
[25,] 0.9363247 1.273507e-01 6.367535e-02
[26,] 0.9357257 1.285486e-01 6.427429e-02
[27,] 0.9619134 7.617318e-02 3.808659e-02
> postscript(file="/var/www/html/rcomp/tmp/1gb5v1227471131.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/2xcdg1227471131.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/3p0ka1227471131.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/4bwzb1227471131.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/5mfww1227471131.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 7
14.1455556 10.3055556 5.3855556 4.6855556 5.9655556 3.9855556 2.8655556
8 9 10 11 12 13 14
1.8655556 -0.9344444 -0.9316667 -1.5316667 -3.2916667 -8.2044444 -9.6444444
15 16 17 18 19 20 21
-8.9644444 -7.6644444 -5.9844444 -2.6644444 0.2155556 0.1155556 1.9155556
22 23 24 25 26 27 28
4.1183333 3.0183333 1.0583333 -1.8544444 -0.7944444 -3.6144444 -2.9144444
29 30 31 32 33 34 35
0.2655556 -0.7144444 -1.7344444 -0.2344444 1.7655556 -6.6455556 -6.2455556
36 37 38 39 40 41 42
-5.3055556 -3.5183333 -0.6583333 2.0216667 2.4216667 -0.3983333 -0.7783333
43 44 45 46 47 48 49
-1.4983333 -4.4983333 -5.1983333 -1.8955556 -0.2955556 -1.6555556 -0.5683333
50 51 52 53 54 55 56
0.7916667 5.1716667 3.4716667 0.1516667 0.1716667 0.1516667 2.7516667
57 58 59 60
2.4516667 5.3544444 5.0544444 9.1944444
> postscript(file="/var/www/html/rcomp/tmp/6jog51227471131.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.1455556 NA
1 10.3055556 14.1455556
2 5.3855556 10.3055556
3 4.6855556 5.3855556
4 5.9655556 4.6855556
5 3.9855556 5.9655556
6 2.8655556 3.9855556
7 1.8655556 2.8655556
8 -0.9344444 1.8655556
9 -0.9316667 -0.9344444
10 -1.5316667 -0.9316667
11 -3.2916667 -1.5316667
12 -8.2044444 -3.2916667
13 -9.6444444 -8.2044444
14 -8.9644444 -9.6444444
15 -7.6644444 -8.9644444
16 -5.9844444 -7.6644444
17 -2.6644444 -5.9844444
18 0.2155556 -2.6644444
19 0.1155556 0.2155556
20 1.9155556 0.1155556
21 4.1183333 1.9155556
22 3.0183333 4.1183333
23 1.0583333 3.0183333
24 -1.8544444 1.0583333
25 -0.7944444 -1.8544444
26 -3.6144444 -0.7944444
27 -2.9144444 -3.6144444
28 0.2655556 -2.9144444
29 -0.7144444 0.2655556
30 -1.7344444 -0.7144444
31 -0.2344444 -1.7344444
32 1.7655556 -0.2344444
33 -6.6455556 1.7655556
34 -6.2455556 -6.6455556
35 -5.3055556 -6.2455556
36 -3.5183333 -5.3055556
37 -0.6583333 -3.5183333
38 2.0216667 -0.6583333
39 2.4216667 2.0216667
40 -0.3983333 2.4216667
41 -0.7783333 -0.3983333
42 -1.4983333 -0.7783333
43 -4.4983333 -1.4983333
44 -5.1983333 -4.4983333
45 -1.8955556 -5.1983333
46 -0.2955556 -1.8955556
47 -1.6555556 -0.2955556
48 -0.5683333 -1.6555556
49 0.7916667 -0.5683333
50 5.1716667 0.7916667
51 3.4716667 5.1716667
52 0.1516667 3.4716667
53 0.1716667 0.1516667
54 0.1516667 0.1716667
55 2.7516667 0.1516667
56 2.4516667 2.7516667
57 5.3544444 2.4516667
58 5.0544444 5.3544444
59 9.1944444 5.0544444
60 NA 9.1944444
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 10.3055556 14.1455556
[2,] 5.3855556 10.3055556
[3,] 4.6855556 5.3855556
[4,] 5.9655556 4.6855556
[5,] 3.9855556 5.9655556
[6,] 2.8655556 3.9855556
[7,] 1.8655556 2.8655556
[8,] -0.9344444 1.8655556
[9,] -0.9316667 -0.9344444
[10,] -1.5316667 -0.9316667
[11,] -3.2916667 -1.5316667
[12,] -8.2044444 -3.2916667
[13,] -9.6444444 -8.2044444
[14,] -8.9644444 -9.6444444
[15,] -7.6644444 -8.9644444
[16,] -5.9844444 -7.6644444
[17,] -2.6644444 -5.9844444
[18,] 0.2155556 -2.6644444
[19,] 0.1155556 0.2155556
[20,] 1.9155556 0.1155556
[21,] 4.1183333 1.9155556
[22,] 3.0183333 4.1183333
[23,] 1.0583333 3.0183333
[24,] -1.8544444 1.0583333
[25,] -0.7944444 -1.8544444
[26,] -3.6144444 -0.7944444
[27,] -2.9144444 -3.6144444
[28,] 0.2655556 -2.9144444
[29,] -0.7144444 0.2655556
[30,] -1.7344444 -0.7144444
[31,] -0.2344444 -1.7344444
[32,] 1.7655556 -0.2344444
[33,] -6.6455556 1.7655556
[34,] -6.2455556 -6.6455556
[35,] -5.3055556 -6.2455556
[36,] -3.5183333 -5.3055556
[37,] -0.6583333 -3.5183333
[38,] 2.0216667 -0.6583333
[39,] 2.4216667 2.0216667
[40,] -0.3983333 2.4216667
[41,] -0.7783333 -0.3983333
[42,] -1.4983333 -0.7783333
[43,] -4.4983333 -1.4983333
[44,] -5.1983333 -4.4983333
[45,] -1.8955556 -5.1983333
[46,] -0.2955556 -1.8955556
[47,] -1.6555556 -0.2955556
[48,] -0.5683333 -1.6555556
[49,] 0.7916667 -0.5683333
[50,] 5.1716667 0.7916667
[51,] 3.4716667 5.1716667
[52,] 0.1516667 3.4716667
[53,] 0.1716667 0.1516667
[54,] 0.1516667 0.1716667
[55,] 2.7516667 0.1516667
[56,] 2.4516667 2.7516667
[57,] 5.3544444 2.4516667
[58,] 5.0544444 5.3544444
[59,] 9.1944444 5.0544444
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 10.3055556 14.1455556
2 5.3855556 10.3055556
3 4.6855556 5.3855556
4 5.9655556 4.6855556
5 3.9855556 5.9655556
6 2.8655556 3.9855556
7 1.8655556 2.8655556
8 -0.9344444 1.8655556
9 -0.9316667 -0.9344444
10 -1.5316667 -0.9316667
11 -3.2916667 -1.5316667
12 -8.2044444 -3.2916667
13 -9.6444444 -8.2044444
14 -8.9644444 -9.6444444
15 -7.6644444 -8.9644444
16 -5.9844444 -7.6644444
17 -2.6644444 -5.9844444
18 0.2155556 -2.6644444
19 0.1155556 0.2155556
20 1.9155556 0.1155556
21 4.1183333 1.9155556
22 3.0183333 4.1183333
23 1.0583333 3.0183333
24 -1.8544444 1.0583333
25 -0.7944444 -1.8544444
26 -3.6144444 -0.7944444
27 -2.9144444 -3.6144444
28 0.2655556 -2.9144444
29 -0.7144444 0.2655556
30 -1.7344444 -0.7144444
31 -0.2344444 -1.7344444
32 1.7655556 -0.2344444
33 -6.6455556 1.7655556
34 -6.2455556 -6.6455556
35 -5.3055556 -6.2455556
36 -3.5183333 -5.3055556
37 -0.6583333 -3.5183333
38 2.0216667 -0.6583333
39 2.4216667 2.0216667
40 -0.3983333 2.4216667
41 -0.7783333 -0.3983333
42 -1.4983333 -0.7783333
43 -4.4983333 -1.4983333
44 -5.1983333 -4.4983333
45 -1.8955556 -5.1983333
46 -0.2955556 -1.8955556
47 -1.6555556 -0.2955556
48 -0.5683333 -1.6555556
49 0.7916667 -0.5683333
50 5.1716667 0.7916667
51 3.4716667 5.1716667
52 0.1516667 3.4716667
53 0.1716667 0.1516667
54 0.1516667 0.1716667
55 2.7516667 0.1516667
56 2.4516667 2.7516667
57 5.3544444 2.4516667
58 5.0544444 5.3544444
59 9.1944444 5.0544444
> 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/7bc2j1227471131.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/89s441227471131.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/9qio41227471131.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/1099xy1227471131.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/11xh8o1227471131.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/12suu11227471131.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/13bra91227471131.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/149kh11227471131.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/155kgc1227471131.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/16kk2x1227471131.tab")
+ }
>
> system("convert tmp/1gb5v1227471131.ps tmp/1gb5v1227471131.png")
> system("convert tmp/2xcdg1227471131.ps tmp/2xcdg1227471131.png")
> system("convert tmp/3p0ka1227471131.ps tmp/3p0ka1227471131.png")
> system("convert tmp/4bwzb1227471131.ps tmp/4bwzb1227471131.png")
> system("convert tmp/5mfww1227471131.ps tmp/5mfww1227471131.png")
> system("convert tmp/6jog51227471131.ps tmp/6jog51227471131.png")
> system("convert tmp/7bc2j1227471131.ps tmp/7bc2j1227471131.png")
> system("convert tmp/89s441227471131.ps tmp/89s441227471131.png")
> system("convert tmp/9qio41227471131.ps tmp/9qio41227471131.png")
> system("convert tmp/1099xy1227471131.ps tmp/1099xy1227471131.png")
>
>
> proc.time()
user system elapsed
2.392 1.531 2.766