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(127.96,0,127.47,0,126.47,0,125.75,0,125.42,0,125.14,0,125.15,0,125.51,0,125.63,0,126.22,0,126.88,0,127.96,0,128.74,0,129.6,0,131.2,0,132.72,0,134.67,0,135.94,0,136.39,0,136.74,0,137.2,0,137.36,0,138.63,0,141.07,0,143.32,0,147.91,0,152.56,0,151.61,0,156.56,0,157.45,0,158.13,0,159.18,0,159.47,0,159.79,0,161.65,0,162.77,0,163.48,0,166.16,0,163.86,0,162.12,0,149.08,0,145.32,0,141.21,0,134.68,0,133.65,0,139.17,0,138.61,0,144.96,1,157.99,1,167.18,1,174.48,1,182.77,1,190.00,1,189.70,1,188.90,1,198.28,1,201.18,1,204.14,1,221.02,1,221.12,1,220.68,1),dim=c(2,61),dimnames=list(c('Gasindex','dumivariable'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Gasindex','dumivariable'),1:61))
> 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
Gasindex dumivariable M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 127.96 0 1 0 0 0 0 0 0 0 0 0 0 1
2 127.47 0 0 1 0 0 0 0 0 0 0 0 0 2
3 126.47 0 0 0 1 0 0 0 0 0 0 0 0 3
4 125.75 0 0 0 0 1 0 0 0 0 0 0 0 4
5 125.42 0 0 0 0 0 1 0 0 0 0 0 0 5
6 125.14 0 0 0 0 0 0 1 0 0 0 0 0 6
7 125.15 0 0 0 0 0 0 0 1 0 0 0 0 7
8 125.51 0 0 0 0 0 0 0 0 1 0 0 0 8
9 125.63 0 0 0 0 0 0 0 0 0 1 0 0 9
10 126.22 0 0 0 0 0 0 0 0 0 0 1 0 10
11 126.88 0 0 0 0 0 0 0 0 0 0 0 1 11
12 127.96 0 0 0 0 0 0 0 0 0 0 0 0 12
13 128.74 0 1 0 0 0 0 0 0 0 0 0 0 13
14 129.60 0 0 1 0 0 0 0 0 0 0 0 0 14
15 131.20 0 0 0 1 0 0 0 0 0 0 0 0 15
16 132.72 0 0 0 0 1 0 0 0 0 0 0 0 16
17 134.67 0 0 0 0 0 1 0 0 0 0 0 0 17
18 135.94 0 0 0 0 0 0 1 0 0 0 0 0 18
19 136.39 0 0 0 0 0 0 0 1 0 0 0 0 19
20 136.74 0 0 0 0 0 0 0 0 1 0 0 0 20
21 137.20 0 0 0 0 0 0 0 0 0 1 0 0 21
22 137.36 0 0 0 0 0 0 0 0 0 0 1 0 22
23 138.63 0 0 0 0 0 0 0 0 0 0 0 1 23
24 141.07 0 0 0 0 0 0 0 0 0 0 0 0 24
25 143.32 0 1 0 0 0 0 0 0 0 0 0 0 25
26 147.91 0 0 1 0 0 0 0 0 0 0 0 0 26
27 152.56 0 0 0 1 0 0 0 0 0 0 0 0 27
28 151.61 0 0 0 0 1 0 0 0 0 0 0 0 28
29 156.56 0 0 0 0 0 1 0 0 0 0 0 0 29
30 157.45 0 0 0 0 0 0 1 0 0 0 0 0 30
31 158.13 0 0 0 0 0 0 0 1 0 0 0 0 31
32 159.18 0 0 0 0 0 0 0 0 1 0 0 0 32
33 159.47 0 0 0 0 0 0 0 0 0 1 0 0 33
34 159.79 0 0 0 0 0 0 0 0 0 0 1 0 34
35 161.65 0 0 0 0 0 0 0 0 0 0 0 1 35
36 162.77 0 0 0 0 0 0 0 0 0 0 0 0 36
37 163.48 0 1 0 0 0 0 0 0 0 0 0 0 37
38 166.16 0 0 1 0 0 0 0 0 0 0 0 0 38
39 163.86 0 0 0 1 0 0 0 0 0 0 0 0 39
40 162.12 0 0 0 0 1 0 0 0 0 0 0 0 40
41 149.08 0 0 0 0 0 1 0 0 0 0 0 0 41
42 145.32 0 0 0 0 0 0 1 0 0 0 0 0 42
43 141.21 0 0 0 0 0 0 0 1 0 0 0 0 43
44 134.68 0 0 0 0 0 0 0 0 1 0 0 0 44
45 133.65 0 0 0 0 0 0 0 0 0 1 0 0 45
46 139.17 0 0 0 0 0 0 0 0 0 0 1 0 46
47 138.61 0 0 0 0 0 0 0 0 0 0 0 1 47
48 144.96 1 0 0 0 0 0 0 0 0 0 0 0 48
49 157.99 1 1 0 0 0 0 0 0 0 0 0 0 49
50 167.18 1 0 1 0 0 0 0 0 0 0 0 0 50
51 174.48 1 0 0 1 0 0 0 0 0 0 0 0 51
52 182.77 1 0 0 0 1 0 0 0 0 0 0 0 52
53 190.00 1 0 0 0 0 1 0 0 0 0 0 0 53
54 189.70 1 0 0 0 0 0 1 0 0 0 0 0 54
55 188.90 1 0 0 0 0 0 0 1 0 0 0 0 55
56 198.28 1 0 0 0 0 0 0 0 1 0 0 0 56
57 201.18 1 0 0 0 0 0 0 0 0 1 0 0 57
58 204.14 1 0 0 0 0 0 0 0 0 0 1 0 58
59 221.02 1 0 0 0 0 0 0 0 0 0 0 1 59
60 221.12 1 0 0 0 0 0 0 0 0 0 0 0 60
61 220.68 1 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) dumivariable M1 M2 M3
119.90184 22.92431 3.21735 1.14632 2.34897
M4 M5 M6 M7 M8
2.78163 2.08628 0.80293 -0.79841 -0.72376
M9 M10 M11 t
-1.02310 0.03955 3.21421 0.84735
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-38.5387 -3.7610 0.1783 8.7263 27.4531
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 119.90184 7.47380 16.043 < 2e-16 ***
dumivariable 22.92431 6.41243 3.575 0.000823 ***
M1 3.21735 8.53786 0.377 0.707994
M2 1.14632 8.96149 0.128 0.898761
M3 2.34897 8.95267 0.262 0.794178
M4 2.78163 8.94648 0.311 0.757237
M5 2.08628 8.94292 0.233 0.816551
M6 0.80293 8.94201 0.090 0.928833
M7 -0.79841 8.94372 -0.089 0.929246
M8 -0.72376 8.94808 -0.081 0.935878
M9 -1.02310 8.95507 -0.114 0.909528
M10 0.03955 8.96468 0.004 0.996498
M11 3.21421 8.97692 0.358 0.721907
t 0.84735 0.15359 5.517 1.44e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.07 on 47 degrees of freedom
Multiple R-squared: 0.7749, Adjusted R-squared: 0.7127
F-statistic: 12.45 on 13 and 47 DF, p-value: 3.858e-11
> 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,] 5.883841e-03 1.176768e-02 0.9941162
[2,] 1.980316e-03 3.960631e-03 0.9980197
[3,] 5.770351e-04 1.154070e-03 0.9994230
[4,] 1.438188e-04 2.876377e-04 0.9998562
[5,] 3.422279e-05 6.844558e-05 0.9999658
[6,] 6.930906e-06 1.386181e-05 0.9999931
[7,] 1.512722e-06 3.025444e-06 0.9999985
[8,] 3.841677e-07 7.683353e-07 0.9999996
[9,] 6.130417e-08 1.226083e-07 0.9999999
[10,] 2.417798e-08 4.835595e-08 1.0000000
[11,] 3.180561e-08 6.361122e-08 1.0000000
[12,] 1.290670e-08 2.581339e-08 1.0000000
[13,] 1.450198e-08 2.900396e-08 1.0000000
[14,] 1.131982e-08 2.263964e-08 1.0000000
[15,] 8.056275e-09 1.611255e-08 1.0000000
[16,] 5.962611e-09 1.192522e-08 1.0000000
[17,] 4.528807e-09 9.057613e-09 1.0000000
[18,] 3.890215e-09 7.780431e-09 1.0000000
[19,] 5.953434e-09 1.190687e-08 1.0000000
[20,] 1.675227e-08 3.350454e-08 1.0000000
[21,] 1.298298e-07 2.596596e-07 0.9999999
[22,] 1.271366e-06 2.542731e-06 0.9999987
[23,] 2.654062e-05 5.308125e-05 0.9999735
[24,] 1.204915e-03 2.409831e-03 0.9987951
[25,] 2.151245e-02 4.302489e-02 0.9784876
[26,] 1.573276e-01 3.146552e-01 0.8426724
[27,] 5.749125e-01 8.501750e-01 0.4250875
[28,] 6.150249e-01 7.699501e-01 0.3849751
> postscript(file="/var/www/html/rcomp/tmp/1o6o91229945793.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/2vfyf1229945793.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/3s3w51229945793.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/403l31229945793.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/5l0wp1229945793.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 = 61
Frequency = 1
1 2 3 4 5 6
3.99346906 4.72715472 1.67715472 -0.32284528 -0.80484528 -0.64884528
7 8 9 10 11 12
0.11515472 -0.44684528 -0.87484528 -2.19484528 -5.55684528 -2.10998371
13 14 15 16 17 18
-5.39467752 -3.31099186 -3.76099186 -3.52099186 -1.72299186 -0.01699186
19 20 21 22 23 24
1.18700814 0.61500814 0.52700814 -1.22299186 -3.97499186 0.83186971
25 26 27 28 29 30
-0.98282410 4.83086156 7.43086156 5.20086156 9.99886156 11.32486156
31 32 33 34 35 36
12.75886156 12.88686156 12.62886156 11.03886156 8.87686156 12.36372313
37 38 39 40 41 42
9.00902932 12.91271498 8.56271498 5.54271498 -7.64928502 -10.97328502
43 44 45 46 47 48
-14.32928502 -21.78128502 -23.35928502 -19.74928502 -24.33128502 -38.53873127
49 50 51 52 53 54
-29.57342508 -19.15973941 -13.90973941 -6.89973941 0.17826059 0.31426059
55 56 57 58 59 60
0.26826059 8.72626059 11.07826059 12.12826059 24.98626059 27.45312215
61
22.94842834
> postscript(file="/var/www/html/rcomp/tmp/6hzu31229945793.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 3.99346906 NA
1 4.72715472 3.99346906
2 1.67715472 4.72715472
3 -0.32284528 1.67715472
4 -0.80484528 -0.32284528
5 -0.64884528 -0.80484528
6 0.11515472 -0.64884528
7 -0.44684528 0.11515472
8 -0.87484528 -0.44684528
9 -2.19484528 -0.87484528
10 -5.55684528 -2.19484528
11 -2.10998371 -5.55684528
12 -5.39467752 -2.10998371
13 -3.31099186 -5.39467752
14 -3.76099186 -3.31099186
15 -3.52099186 -3.76099186
16 -1.72299186 -3.52099186
17 -0.01699186 -1.72299186
18 1.18700814 -0.01699186
19 0.61500814 1.18700814
20 0.52700814 0.61500814
21 -1.22299186 0.52700814
22 -3.97499186 -1.22299186
23 0.83186971 -3.97499186
24 -0.98282410 0.83186971
25 4.83086156 -0.98282410
26 7.43086156 4.83086156
27 5.20086156 7.43086156
28 9.99886156 5.20086156
29 11.32486156 9.99886156
30 12.75886156 11.32486156
31 12.88686156 12.75886156
32 12.62886156 12.88686156
33 11.03886156 12.62886156
34 8.87686156 11.03886156
35 12.36372313 8.87686156
36 9.00902932 12.36372313
37 12.91271498 9.00902932
38 8.56271498 12.91271498
39 5.54271498 8.56271498
40 -7.64928502 5.54271498
41 -10.97328502 -7.64928502
42 -14.32928502 -10.97328502
43 -21.78128502 -14.32928502
44 -23.35928502 -21.78128502
45 -19.74928502 -23.35928502
46 -24.33128502 -19.74928502
47 -38.53873127 -24.33128502
48 -29.57342508 -38.53873127
49 -19.15973941 -29.57342508
50 -13.90973941 -19.15973941
51 -6.89973941 -13.90973941
52 0.17826059 -6.89973941
53 0.31426059 0.17826059
54 0.26826059 0.31426059
55 8.72626059 0.26826059
56 11.07826059 8.72626059
57 12.12826059 11.07826059
58 24.98626059 12.12826059
59 27.45312215 24.98626059
60 22.94842834 27.45312215
61 NA 22.94842834
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4.72715472 3.99346906
[2,] 1.67715472 4.72715472
[3,] -0.32284528 1.67715472
[4,] -0.80484528 -0.32284528
[5,] -0.64884528 -0.80484528
[6,] 0.11515472 -0.64884528
[7,] -0.44684528 0.11515472
[8,] -0.87484528 -0.44684528
[9,] -2.19484528 -0.87484528
[10,] -5.55684528 -2.19484528
[11,] -2.10998371 -5.55684528
[12,] -5.39467752 -2.10998371
[13,] -3.31099186 -5.39467752
[14,] -3.76099186 -3.31099186
[15,] -3.52099186 -3.76099186
[16,] -1.72299186 -3.52099186
[17,] -0.01699186 -1.72299186
[18,] 1.18700814 -0.01699186
[19,] 0.61500814 1.18700814
[20,] 0.52700814 0.61500814
[21,] -1.22299186 0.52700814
[22,] -3.97499186 -1.22299186
[23,] 0.83186971 -3.97499186
[24,] -0.98282410 0.83186971
[25,] 4.83086156 -0.98282410
[26,] 7.43086156 4.83086156
[27,] 5.20086156 7.43086156
[28,] 9.99886156 5.20086156
[29,] 11.32486156 9.99886156
[30,] 12.75886156 11.32486156
[31,] 12.88686156 12.75886156
[32,] 12.62886156 12.88686156
[33,] 11.03886156 12.62886156
[34,] 8.87686156 11.03886156
[35,] 12.36372313 8.87686156
[36,] 9.00902932 12.36372313
[37,] 12.91271498 9.00902932
[38,] 8.56271498 12.91271498
[39,] 5.54271498 8.56271498
[40,] -7.64928502 5.54271498
[41,] -10.97328502 -7.64928502
[42,] -14.32928502 -10.97328502
[43,] -21.78128502 -14.32928502
[44,] -23.35928502 -21.78128502
[45,] -19.74928502 -23.35928502
[46,] -24.33128502 -19.74928502
[47,] -38.53873127 -24.33128502
[48,] -29.57342508 -38.53873127
[49,] -19.15973941 -29.57342508
[50,] -13.90973941 -19.15973941
[51,] -6.89973941 -13.90973941
[52,] 0.17826059 -6.89973941
[53,] 0.31426059 0.17826059
[54,] 0.26826059 0.31426059
[55,] 8.72626059 0.26826059
[56,] 11.07826059 8.72626059
[57,] 12.12826059 11.07826059
[58,] 24.98626059 12.12826059
[59,] 27.45312215 24.98626059
[60,] 22.94842834 27.45312215
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4.72715472 3.99346906
2 1.67715472 4.72715472
3 -0.32284528 1.67715472
4 -0.80484528 -0.32284528
5 -0.64884528 -0.80484528
6 0.11515472 -0.64884528
7 -0.44684528 0.11515472
8 -0.87484528 -0.44684528
9 -2.19484528 -0.87484528
10 -5.55684528 -2.19484528
11 -2.10998371 -5.55684528
12 -5.39467752 -2.10998371
13 -3.31099186 -5.39467752
14 -3.76099186 -3.31099186
15 -3.52099186 -3.76099186
16 -1.72299186 -3.52099186
17 -0.01699186 -1.72299186
18 1.18700814 -0.01699186
19 0.61500814 1.18700814
20 0.52700814 0.61500814
21 -1.22299186 0.52700814
22 -3.97499186 -1.22299186
23 0.83186971 -3.97499186
24 -0.98282410 0.83186971
25 4.83086156 -0.98282410
26 7.43086156 4.83086156
27 5.20086156 7.43086156
28 9.99886156 5.20086156
29 11.32486156 9.99886156
30 12.75886156 11.32486156
31 12.88686156 12.75886156
32 12.62886156 12.88686156
33 11.03886156 12.62886156
34 8.87686156 11.03886156
35 12.36372313 8.87686156
36 9.00902932 12.36372313
37 12.91271498 9.00902932
38 8.56271498 12.91271498
39 5.54271498 8.56271498
40 -7.64928502 5.54271498
41 -10.97328502 -7.64928502
42 -14.32928502 -10.97328502
43 -21.78128502 -14.32928502
44 -23.35928502 -21.78128502
45 -19.74928502 -23.35928502
46 -24.33128502 -19.74928502
47 -38.53873127 -24.33128502
48 -29.57342508 -38.53873127
49 -19.15973941 -29.57342508
50 -13.90973941 -19.15973941
51 -6.89973941 -13.90973941
52 0.17826059 -6.89973941
53 0.31426059 0.17826059
54 0.26826059 0.31426059
55 8.72626059 0.26826059
56 11.07826059 8.72626059
57 12.12826059 11.07826059
58 24.98626059 12.12826059
59 27.45312215 24.98626059
60 22.94842834 27.45312215
> 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/7rlkr1229945793.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/8bvpr1229945793.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/9ie3h1229945794.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/10gia91229945794.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/11aibj1229945794.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/12jjb31229945794.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/131w8t1229945794.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/14ln8f1229945794.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/15hgb41229945794.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/16vs1x1229945794.tab")
+ }
>
> system("convert tmp/1o6o91229945793.ps tmp/1o6o91229945793.png")
> system("convert tmp/2vfyf1229945793.ps tmp/2vfyf1229945793.png")
> system("convert tmp/3s3w51229945793.ps tmp/3s3w51229945793.png")
> system("convert tmp/403l31229945793.ps tmp/403l31229945793.png")
> system("convert tmp/5l0wp1229945793.ps tmp/5l0wp1229945793.png")
> system("convert tmp/6hzu31229945793.ps tmp/6hzu31229945793.png")
> system("convert tmp/7rlkr1229945793.ps tmp/7rlkr1229945793.png")
> system("convert tmp/8bvpr1229945793.ps tmp/8bvpr1229945793.png")
> system("convert tmp/9ie3h1229945794.ps tmp/9ie3h1229945794.png")
> system("convert tmp/10gia91229945794.ps tmp/10gia91229945794.png")
>
>
> proc.time()
user system elapsed
2.436 1.588 3.086