R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(1.464
+ ,1.487
+ ,1.232
+ ,0.6683
+ ,0.659
+ ,1.464
+ ,1.487
+ ,1.232
+ ,0.6683
+ ,0.659
+ ,1.464
+ ,1.487
+ ,1.268
+ ,0.7041
+ ,0.659
+ ,1.464
+ ,1.487
+ ,1.268
+ ,0.7041
+ ,0.659
+ ,1.464
+ ,1.487
+ ,1.268
+ ,0.7041
+ ,0.659
+ ,1.464
+ ,1.487
+ ,1.268
+ ,0.7041
+ ,0.659
+ ,1.498
+ ,1.523
+ ,1.268
+ ,0.7041
+ ,0.659
+ ,1.498
+ ,1.523
+ ,1.268
+ ,0.7041
+ ,0.659
+ ,1.498
+ ,1.523
+ ,1.268
+ ,0.7041
+ ,0.659
+ ,1.498
+ ,1.523
+ ,1.268
+ ,0.7041
+ ,0.732
+ ,1.498
+ ,1.523
+ ,1.268
+ ,0.7041
+ ,0.732
+ ,1.498
+ ,1.523
+ ,1.268
+ ,0.7041
+ ,0.732
+ ,1.498
+ ,1.523
+ ,1.298
+ ,0.7319
+ ,0.732
+ ,1.538
+ ,1.562
+ ,1.298
+ ,0.7319
+ ,0.732
+ ,1.538
+ ,1.562
+ ,1.298
+ ,0.7319
+ ,0.732
+ ,1.538
+ ,1.562
+ ,1.298
+ ,0.7319
+ ,0.732
+ ,1.538
+ ,1.562
+ ,1.298
+ ,0.7319
+ ,0.732
+ ,1.538
+ ,1.562
+ ,1.298
+ ,0.7319
+ ,0.732
+ ,1.538
+ ,1.562
+ ,1.298
+ ,0.7319
+ ,0.732
+ ,1.538
+ ,1.562
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.538
+ ,1.562
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.511
+ ,1.535
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.511
+ ,1.535
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.511
+ ,1.535
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.511
+ ,1.535
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.511
+ ,1.535
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.511
+ ,1.535
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.511
+ ,1.535
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.511
+ ,1.535
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.511
+ ,1.535
+ ,1.298
+ ,0.7319
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.324
+ ,0.7568
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.324
+ ,0.7568
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.324
+ ,0.7568
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.324
+ ,0.7568
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.324
+ ,0.7568
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.324
+ ,0.7568
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.324
+ ,0.7568
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.324
+ ,0.7568
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.324
+ ,0.7568
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.324
+ ,0.7568
+ ,0.725
+ ,1.544
+ ,1.569
+ ,1.31
+ ,0.7286
+ ,0.725
+ ,1.524
+ ,1.548
+ ,1.31
+ ,0.7286
+ ,0.725
+ ,1.524
+ ,1.548
+ ,1.31
+ ,0.7286
+ ,0.725
+ ,1.524
+ ,1.548
+ ,1.31
+ ,0.7286
+ ,0.725
+ ,1.524
+ ,1.548
+ ,1.31
+ ,0.7286
+ ,0.725
+ ,1.524
+ ,1.548
+ ,1.31
+ ,0.7286
+ ,0.725
+ ,1.524
+ ,1.548
+ ,1.31
+ ,0.7286
+ ,0.725
+ ,1.524
+ ,1.548
+ ,1.338
+ ,0.753
+ ,0.725
+ ,1.558
+ ,1.584
+ ,1.338
+ ,0.753
+ ,0.725)
+ ,dim=c(5
+ ,49)
+ ,dimnames=list(c('Super95'
+ ,'Super98'
+ ,'Diesel'
+ ,'Gasolie'
+ ,'LPG')
+ ,1:49))
> y <- array(NA,dim=c(5,49),dimnames=list(c('Super95','Super98','Diesel','Gasolie','LPG'),1:49))
> 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
> 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
Super95 Super98 Diesel Gasolie LPG
1 1.464 1.487 1.232 0.6683 0.659
2 1.464 1.487 1.232 0.6683 0.659
3 1.464 1.487 1.268 0.7041 0.659
4 1.464 1.487 1.268 0.7041 0.659
5 1.464 1.487 1.268 0.7041 0.659
6 1.464 1.487 1.268 0.7041 0.659
7 1.498 1.523 1.268 0.7041 0.659
8 1.498 1.523 1.268 0.7041 0.659
9 1.498 1.523 1.268 0.7041 0.659
10 1.498 1.523 1.268 0.7041 0.732
11 1.498 1.523 1.268 0.7041 0.732
12 1.498 1.523 1.268 0.7041 0.732
13 1.498 1.523 1.298 0.7319 0.732
14 1.538 1.562 1.298 0.7319 0.732
15 1.538 1.562 1.298 0.7319 0.732
16 1.538 1.562 1.298 0.7319 0.732
17 1.538 1.562 1.298 0.7319 0.732
18 1.538 1.562 1.298 0.7319 0.732
19 1.538 1.562 1.298 0.7319 0.732
20 1.538 1.562 1.298 0.7319 0.725
21 1.538 1.562 1.298 0.7319 0.725
22 1.511 1.535 1.298 0.7319 0.725
23 1.511 1.535 1.298 0.7319 0.725
24 1.511 1.535 1.298 0.7319 0.725
25 1.511 1.535 1.298 0.7319 0.725
26 1.511 1.535 1.298 0.7319 0.725
27 1.511 1.535 1.298 0.7319 0.725
28 1.511 1.535 1.298 0.7319 0.725
29 1.511 1.535 1.298 0.7319 0.725
30 1.511 1.535 1.298 0.7319 0.725
31 1.544 1.569 1.324 0.7568 0.725
32 1.544 1.569 1.324 0.7568 0.725
33 1.544 1.569 1.324 0.7568 0.725
34 1.544 1.569 1.324 0.7568 0.725
35 1.544 1.569 1.324 0.7568 0.725
36 1.544 1.569 1.324 0.7568 0.725
37 1.544 1.569 1.324 0.7568 0.725
38 1.544 1.569 1.324 0.7568 0.725
39 1.544 1.569 1.324 0.7568 0.725
40 1.544 1.569 1.324 0.7568 0.725
41 1.544 1.569 1.310 0.7286 0.725
42 1.524 1.548 1.310 0.7286 0.725
43 1.524 1.548 1.310 0.7286 0.725
44 1.524 1.548 1.310 0.7286 0.725
45 1.524 1.548 1.310 0.7286 0.725
46 1.524 1.548 1.310 0.7286 0.725
47 1.524 1.548 1.310 0.7286 0.725
48 1.524 1.548 1.338 0.7530 0.725
49 1.558 1.584 1.338 0.7530 0.725
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Super98 Diesel Gasolie LPG
-0.002484 0.979475 0.013001 -0.013389 0.003840
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.0011432 -0.0001768 0.0001300 0.0004048 0.0006842
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.002484 0.007972 -0.312 0.757
Super98 0.979475 0.006689 146.425 <2e-16 ***
Diesel 0.013001 0.013281 0.979 0.333
Gasolie -0.013389 0.014757 -0.907 0.369
LPG 0.003840 0.004725 0.813 0.421
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0005615 on 44 degrees of freedom
Multiple R-squared: 0.9996, Adjusted R-squared: 0.9995
F-statistic: 2.634e+04 on 4 and 44 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,] 7.983658e-40 1.596732e-39 1.000000e+00
[2,] 8.960172e-50 1.792034e-49 1.000000e+00
[3,] 2.168803e-72 4.337606e-72 1.000000e+00
[4,] 5.489477e-75 1.097895e-74 1.000000e+00
[5,] 1.329175e-85 2.658351e-85 1.000000e+00
[6,] 4.347466e-111 8.694931e-111 1.000000e+00
[7,] 9.999923e-01 1.547697e-05 7.738484e-06
[8,] 9.999968e-01 6.494354e-06 3.247177e-06
[9,] 9.999954e-01 9.297443e-06 4.648722e-06
[10,] 9.999907e-01 1.864408e-05 9.322041e-06
[11,] 9.999790e-01 4.195386e-05 2.097693e-05
[12,] 9.999511e-01 9.775126e-05 4.887563e-05
[13,] 9.999799e-01 4.021757e-05 2.010878e-05
[14,] 9.999998e-01 4.194103e-07 2.097051e-07
[15,] 9.999995e-01 9.421667e-07 4.710833e-07
[16,] 9.999989e-01 2.267283e-06 1.133642e-06
[17,] 9.999973e-01 5.460277e-06 2.730139e-06
[18,] 9.999937e-01 1.259293e-05 6.296465e-06
[19,] 9.999867e-01 2.660007e-05 1.330003e-05
[20,] 9.999761e-01 4.786858e-05 2.393429e-05
[21,] 9.999692e-01 6.150817e-05 3.075408e-05
[22,] 9.999859e-01 2.818283e-05 1.409142e-05
[23,] 1.000000e+00 3.773509e-14 1.886755e-14
[24,] 1.000000e+00 2.791619e-13 1.395809e-13
[25,] 1.000000e+00 2.917411e-12 1.458706e-12
[26,] 1.000000e+00 3.496463e-11 1.748232e-11
[27,] 1.000000e+00 4.402373e-10 2.201186e-10
[28,] 1.000000e+00 5.568225e-09 2.784113e-09
[29,] 1.000000e+00 6.888119e-08 3.444059e-08
[30,] 9.999996e-01 8.180814e-07 4.090407e-07
[31,] 9.999954e-01 9.190285e-06 4.595142e-06
[32,] 9.999518e-01 9.630039e-05 4.815019e-05
[33,] 9.995364e-01 9.271171e-04 4.635586e-04
[34,] 1.000000e+00 2.674714e-37 1.337357e-37
> postscript(file="/var/wessaorg/rcomp/tmp/1l4cz1322054586.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/27udy1322054586.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/3zupp1322054586.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/4vg4a1322054586.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/53do61322054586.ps",horizontal=F,onefile=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 = 49
Frequency = 1
1 2 3 4 5
0.0004048121 0.0004048121 0.0004160890 0.0004160890 0.0004160890
6 7 8 9 10
0.0004160890 -0.0008449959 -0.0008449959 -0.0008449959 -0.0011253381
11 12 13 14 15
-0.0011253381 -0.0011253381 -0.0011431639 0.0006573274 0.0006573274
16 17 18 19 20
0.0006573274 0.0006573274 0.0006573274 0.0006573274 0.0006842096
21 22 23 24 25
0.0006842096 0.0001300232 0.0001300232 0.0001300232 0.0001300232
26 27 28 29 30
0.0001300232 0.0001300232 0.0001300232 0.0001300232 0.0001300232
31 32 33 34 35
-0.0001767615 -0.0001767615 -0.0001767615 -0.0001767615 -0.0001767615
36 37 38 39 40
-0.0001767615 -0.0001767615 -0.0001767615 -0.0001767615 -0.0001767615
41 42 43 44 45
-0.0003723052 0.0001966610 0.0001966610 0.0001966610 0.0001966610
46 47 48 49
0.0001966610 0.0001966610 0.0001593160 -0.0011017689
> postscript(file="/var/wessaorg/rcomp/tmp/66evg1322054586.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 49
Frequency = 1
lag(myerror, k = 1) myerror
0 0.0004048121 NA
1 0.0004048121 0.0004048121
2 0.0004160890 0.0004048121
3 0.0004160890 0.0004160890
4 0.0004160890 0.0004160890
5 0.0004160890 0.0004160890
6 -0.0008449959 0.0004160890
7 -0.0008449959 -0.0008449959
8 -0.0008449959 -0.0008449959
9 -0.0011253381 -0.0008449959
10 -0.0011253381 -0.0011253381
11 -0.0011253381 -0.0011253381
12 -0.0011431639 -0.0011253381
13 0.0006573274 -0.0011431639
14 0.0006573274 0.0006573274
15 0.0006573274 0.0006573274
16 0.0006573274 0.0006573274
17 0.0006573274 0.0006573274
18 0.0006573274 0.0006573274
19 0.0006842096 0.0006573274
20 0.0006842096 0.0006842096
21 0.0001300232 0.0006842096
22 0.0001300232 0.0001300232
23 0.0001300232 0.0001300232
24 0.0001300232 0.0001300232
25 0.0001300232 0.0001300232
26 0.0001300232 0.0001300232
27 0.0001300232 0.0001300232
28 0.0001300232 0.0001300232
29 0.0001300232 0.0001300232
30 -0.0001767615 0.0001300232
31 -0.0001767615 -0.0001767615
32 -0.0001767615 -0.0001767615
33 -0.0001767615 -0.0001767615
34 -0.0001767615 -0.0001767615
35 -0.0001767615 -0.0001767615
36 -0.0001767615 -0.0001767615
37 -0.0001767615 -0.0001767615
38 -0.0001767615 -0.0001767615
39 -0.0001767615 -0.0001767615
40 -0.0003723052 -0.0001767615
41 0.0001966610 -0.0003723052
42 0.0001966610 0.0001966610
43 0.0001966610 0.0001966610
44 0.0001966610 0.0001966610
45 0.0001966610 0.0001966610
46 0.0001966610 0.0001966610
47 0.0001593160 0.0001966610
48 -0.0011017689 0.0001593160
49 NA -0.0011017689
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.0004048121 0.0004048121
[2,] 0.0004160890 0.0004048121
[3,] 0.0004160890 0.0004160890
[4,] 0.0004160890 0.0004160890
[5,] 0.0004160890 0.0004160890
[6,] -0.0008449959 0.0004160890
[7,] -0.0008449959 -0.0008449959
[8,] -0.0008449959 -0.0008449959
[9,] -0.0011253381 -0.0008449959
[10,] -0.0011253381 -0.0011253381
[11,] -0.0011253381 -0.0011253381
[12,] -0.0011431639 -0.0011253381
[13,] 0.0006573274 -0.0011431639
[14,] 0.0006573274 0.0006573274
[15,] 0.0006573274 0.0006573274
[16,] 0.0006573274 0.0006573274
[17,] 0.0006573274 0.0006573274
[18,] 0.0006573274 0.0006573274
[19,] 0.0006842096 0.0006573274
[20,] 0.0006842096 0.0006842096
[21,] 0.0001300232 0.0006842096
[22,] 0.0001300232 0.0001300232
[23,] 0.0001300232 0.0001300232
[24,] 0.0001300232 0.0001300232
[25,] 0.0001300232 0.0001300232
[26,] 0.0001300232 0.0001300232
[27,] 0.0001300232 0.0001300232
[28,] 0.0001300232 0.0001300232
[29,] 0.0001300232 0.0001300232
[30,] -0.0001767615 0.0001300232
[31,] -0.0001767615 -0.0001767615
[32,] -0.0001767615 -0.0001767615
[33,] -0.0001767615 -0.0001767615
[34,] -0.0001767615 -0.0001767615
[35,] -0.0001767615 -0.0001767615
[36,] -0.0001767615 -0.0001767615
[37,] -0.0001767615 -0.0001767615
[38,] -0.0001767615 -0.0001767615
[39,] -0.0001767615 -0.0001767615
[40,] -0.0003723052 -0.0001767615
[41,] 0.0001966610 -0.0003723052
[42,] 0.0001966610 0.0001966610
[43,] 0.0001966610 0.0001966610
[44,] 0.0001966610 0.0001966610
[45,] 0.0001966610 0.0001966610
[46,] 0.0001966610 0.0001966610
[47,] 0.0001593160 0.0001966610
[48,] -0.0011017689 0.0001593160
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.0004048121 0.0004048121
2 0.0004160890 0.0004048121
3 0.0004160890 0.0004160890
4 0.0004160890 0.0004160890
5 0.0004160890 0.0004160890
6 -0.0008449959 0.0004160890
7 -0.0008449959 -0.0008449959
8 -0.0008449959 -0.0008449959
9 -0.0011253381 -0.0008449959
10 -0.0011253381 -0.0011253381
11 -0.0011253381 -0.0011253381
12 -0.0011431639 -0.0011253381
13 0.0006573274 -0.0011431639
14 0.0006573274 0.0006573274
15 0.0006573274 0.0006573274
16 0.0006573274 0.0006573274
17 0.0006573274 0.0006573274
18 0.0006573274 0.0006573274
19 0.0006842096 0.0006573274
20 0.0006842096 0.0006842096
21 0.0001300232 0.0006842096
22 0.0001300232 0.0001300232
23 0.0001300232 0.0001300232
24 0.0001300232 0.0001300232
25 0.0001300232 0.0001300232
26 0.0001300232 0.0001300232
27 0.0001300232 0.0001300232
28 0.0001300232 0.0001300232
29 0.0001300232 0.0001300232
30 -0.0001767615 0.0001300232
31 -0.0001767615 -0.0001767615
32 -0.0001767615 -0.0001767615
33 -0.0001767615 -0.0001767615
34 -0.0001767615 -0.0001767615
35 -0.0001767615 -0.0001767615
36 -0.0001767615 -0.0001767615
37 -0.0001767615 -0.0001767615
38 -0.0001767615 -0.0001767615
39 -0.0001767615 -0.0001767615
40 -0.0003723052 -0.0001767615
41 0.0001966610 -0.0003723052
42 0.0001966610 0.0001966610
43 0.0001966610 0.0001966610
44 0.0001966610 0.0001966610
45 0.0001966610 0.0001966610
46 0.0001966610 0.0001966610
47 0.0001593160 0.0001966610
48 -0.0011017689 0.0001593160
> 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/wessaorg/rcomp/tmp/7dtke1322054586.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/8y0zj1322054586.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/95cic1322054586.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/106wb51322054586.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11q0um1322054586.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/wessaorg/rcomp/tmp/12gzxk1322054586.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/wessaorg/rcomp/tmp/13dlqj1322054586.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/wessaorg/rcomp/tmp/14zm0d1322054586.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/wessaorg/rcomp/tmp/15fkse1322054586.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/wessaorg/rcomp/tmp/16lrxq1322054586.tab")
+ }
>
> try(system("convert tmp/1l4cz1322054586.ps tmp/1l4cz1322054586.png",intern=TRUE))
character(0)
> try(system("convert tmp/27udy1322054586.ps tmp/27udy1322054586.png",intern=TRUE))
character(0)
> try(system("convert tmp/3zupp1322054586.ps tmp/3zupp1322054586.png",intern=TRUE))
character(0)
> try(system("convert tmp/4vg4a1322054586.ps tmp/4vg4a1322054586.png",intern=TRUE))
character(0)
> try(system("convert tmp/53do61322054586.ps tmp/53do61322054586.png",intern=TRUE))
character(0)
> try(system("convert tmp/66evg1322054586.ps tmp/66evg1322054586.png",intern=TRUE))
character(0)
> try(system("convert tmp/7dtke1322054586.ps tmp/7dtke1322054586.png",intern=TRUE))
character(0)
> try(system("convert tmp/8y0zj1322054586.ps tmp/8y0zj1322054586.png",intern=TRUE))
character(0)
> try(system("convert tmp/95cic1322054586.ps tmp/95cic1322054586.png",intern=TRUE))
character(0)
> try(system("convert tmp/106wb51322054586.ps tmp/106wb51322054586.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.036 0.480 3.546