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(8.8,8.1,0,8.5,9.9,0,8.6,11.5,0,8.7,23.4,0,9.1,25.4,0,8.8,27.9,0,6.3,26.1,0,2.5,18.8,0,-2.7,14.1,0,-4.5,11.5,0,-7,15.8,0,-9.3,12.4,0,-12.2,4.5,0,-13.2,-2.2,1,-13.7,-4.2,1,-15,-9.4,1,-16.9,-14.5,1,-16.3,-17.9,1,-16.7,-15.1,1,-16,-15.2,1,-14.5,-15.7,1,-12.2,-18,1,-7.5,-18.1,1,-4.4,-13.5,1,-1.1,-9.9,1,1.3,-4.8,1,-0.1,-1.7,0,0.4,-0.1,0,2.4,2.2,0,1,10.2,0,3.3,7.6,0,1.8,10.8,0,3.2,3.8,0,1.3,11,0,1.5,10.8,0,1.3,20.1,0,2,14.9,0,3,13,0,4.4,10.9,0,3.1,9.6,0,2.6,4,0,2.7,-1.1,0,4,-7.7,0,4.1,-8.9,0,3,-8,0,2.7,-7.1,0,4,-5.3,0,4.8,-2.5,0,6,-2.4,0,4.6,-2.9,0,4.4,-4.8,0,6.6,-7.2,0,4.7,1.7,0,7.6,2.2,0,5.3,13.4,0,6.6,12.3,0,4,13.7,0,3.8,4.4,0,1.2,-2.5,0),dim=c(3,59),dimnames=list(c('Industriële_productie','registratie_personenwagens','crisis'),1:59))
> y <- array(NA,dim=c(3,59),dimnames=list(c('Industriële_productie','registratie_personenwagens','crisis'),1:59))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> 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
Industri\303\253le_productie registratie_personenwagens crisis M1 M2 M3 M4
1 8.8 8.1 0 1 0 0 0
2 8.5 9.9 0 0 1 0 0
3 8.6 11.5 0 0 0 1 0
4 8.7 23.4 0 0 0 0 1
5 9.1 25.4 0 0 0 0 0
6 8.8 27.9 0 0 0 0 0
7 6.3 26.1 0 0 0 0 0
8 2.5 18.8 0 0 0 0 0
9 -2.7 14.1 0 0 0 0 0
10 -4.5 11.5 0 0 0 0 0
11 -7.0 15.8 0 0 0 0 0
12 -9.3 12.4 0 0 0 0 0
13 -12.2 4.5 0 1 0 0 0
14 -13.2 -2.2 1 0 1 0 0
15 -13.7 -4.2 1 0 0 1 0
16 -15.0 -9.4 1 0 0 0 1
17 -16.9 -14.5 1 0 0 0 0
18 -16.3 -17.9 1 0 0 0 0
19 -16.7 -15.1 1 0 0 0 0
20 -16.0 -15.2 1 0 0 0 0
21 -14.5 -15.7 1 0 0 0 0
22 -12.2 -18.0 1 0 0 0 0
23 -7.5 -18.1 1 0 0 0 0
24 -4.4 -13.5 1 0 0 0 0
25 -1.1 -9.9 1 1 0 0 0
26 1.3 -4.8 1 0 1 0 0
27 -0.1 -1.7 0 0 0 1 0
28 0.4 -0.1 0 0 0 0 1
29 2.4 2.2 0 0 0 0 0
30 1.0 10.2 0 0 0 0 0
31 3.3 7.6 0 0 0 0 0
32 1.8 10.8 0 0 0 0 0
33 3.2 3.8 0 0 0 0 0
34 1.3 11.0 0 0 0 0 0
35 1.5 10.8 0 0 0 0 0
36 1.3 20.1 0 0 0 0 0
37 2.0 14.9 0 1 0 0 0
38 3.0 13.0 0 0 1 0 0
39 4.4 10.9 0 0 0 1 0
40 3.1 9.6 0 0 0 0 1
41 2.6 4.0 0 0 0 0 0
42 2.7 -1.1 0 0 0 0 0
43 4.0 -7.7 0 0 0 0 0
44 4.1 -8.9 0 0 0 0 0
45 3.0 -8.0 0 0 0 0 0
46 2.7 -7.1 0 0 0 0 0
47 4.0 -5.3 0 0 0 0 0
48 4.8 -2.5 0 0 0 0 0
49 6.0 -2.4 0 1 0 0 0
50 4.6 -2.9 0 0 1 0 0
51 4.4 -4.8 0 0 0 1 0
52 6.6 -7.2 0 0 0 0 1
53 4.7 1.7 0 0 0 0 0
54 7.6 2.2 0 0 0 0 0
55 5.3 13.4 0 0 0 0 0
56 6.6 12.3 0 0 0 0 0
57 4.0 13.7 0 0 0 0 0
58 3.8 4.4 0 0 0 0 0
59 1.2 -2.5 0 0 0 0 0
M5 M6 M7 M8 M9 M10 M11
1 0 0 0 0 0 0 0
2 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0
4 0 0 0 0 0 0 0
5 1 0 0 0 0 0 0
6 0 1 0 0 0 0 0
7 0 0 1 0 0 0 0
8 0 0 0 1 0 0 0
9 0 0 0 0 1 0 0
10 0 0 0 0 0 1 0
11 0 0 0 0 0 0 1
12 0 0 0 0 0 0 0
13 0 0 0 0 0 0 0
14 0 0 0 0 0 0 0
15 0 0 0 0 0 0 0
16 0 0 0 0 0 0 0
17 1 0 0 0 0 0 0
18 0 1 0 0 0 0 0
19 0 0 1 0 0 0 0
20 0 0 0 1 0 0 0
21 0 0 0 0 1 0 0
22 0 0 0 0 0 1 0
23 0 0 0 0 0 0 1
24 0 0 0 0 0 0 0
25 0 0 0 0 0 0 0
26 0 0 0 0 0 0 0
27 0 0 0 0 0 0 0
28 0 0 0 0 0 0 0
29 1 0 0 0 0 0 0
30 0 1 0 0 0 0 0
31 0 0 1 0 0 0 0
32 0 0 0 1 0 0 0
33 0 0 0 0 1 0 0
34 0 0 0 0 0 1 0
35 0 0 0 0 0 0 1
36 0 0 0 0 0 0 0
37 0 0 0 0 0 0 0
38 0 0 0 0 0 0 0
39 0 0 0 0 0 0 0
40 0 0 0 0 0 0 0
41 1 0 0 0 0 0 0
42 0 1 0 0 0 0 0
43 0 0 1 0 0 0 0
44 0 0 0 1 0 0 0
45 0 0 0 0 1 0 0
46 0 0 0 0 0 1 0
47 0 0 0 0 0 0 1
48 0 0 0 0 0 0 0
49 0 0 0 0 0 0 0
50 0 0 0 0 0 0 0
51 0 0 0 0 0 0 0
52 0 0 0 0 0 0 0
53 1 0 0 0 0 0 0
54 0 1 0 0 0 0 0
55 0 0 1 0 0 0 0
56 0 0 0 1 0 0 0
57 0 0 0 0 1 0 0
58 0 0 0 0 0 1 0
59 0 0 0 0 0 0 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) registratie_personenwagens
1.52136 0.02403
crisis M1
-14.08189 1.92198
M2 M3
4.88893 1.95880
M4 M5
1.97669 1.58468
M6 M7
1.95266 1.61824
M8 M9
1.00948 -0.14294
M10 M11
-0.49363 -0.26834
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-15.75146 -2.85822 0.00796 2.64044 9.77644
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.52136 2.74460 0.554 0.582
registratie_personenwagens 0.02403 0.07858 0.306 0.761
crisis -14.08189 2.26455 -6.218 1.48e-07 ***
M1 1.92198 3.51810 0.546 0.588
M2 4.88893 3.52394 1.387 0.172
M3 1.95880 3.52106 0.556 0.581
M4 1.97669 3.51734 0.562 0.577
M5 1.58468 3.51594 0.451 0.654
M6 1.95266 3.51498 0.556 0.581
M7 1.61824 3.51441 0.460 0.647
M8 1.00948 3.51645 0.287 0.775
M9 -0.14294 3.52525 -0.041 0.968
M10 -0.49363 3.53408 -0.140 0.890
M11 -0.26834 3.53594 -0.076 0.940
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.237 on 45 degrees of freedom
Multiple R-squared: 0.6381, Adjusted R-squared: 0.5336
F-statistic: 6.104 on 13 and 45 DF, p-value: 2.237e-06
> 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.9934970 1.300601e-02 6.503004e-03
[2,] 0.9899827 2.003468e-02 1.001734e-02
[3,] 0.9907819 1.843622e-02 9.218111e-03
[4,] 0.9961455 7.708968e-03 3.854484e-03
[5,] 0.9996154 7.692550e-04 3.846275e-04
[6,] 0.9999854 2.920854e-05 1.460427e-05
[7,] 0.9999994 1.240251e-06 6.201256e-07
[8,] 0.9999999 1.281429e-07 6.407144e-08
[9,] 0.9999999 1.801671e-07 9.008357e-08
[10,] 0.9999997 5.165339e-07 2.582670e-07
[11,] 0.9999998 4.375751e-07 2.187876e-07
[12,] 0.9999998 3.526232e-07 1.763116e-07
[13,] 0.9999994 1.194387e-06 5.971934e-07
[14,] 0.9999991 1.835057e-06 9.175286e-07
[15,] 0.9999968 6.342566e-06 3.171283e-06
[16,] 0.9999943 1.135963e-05 5.679815e-06
[17,] 0.9999810 3.801387e-05 1.900694e-05
[18,] 0.9999377 1.245522e-04 6.227609e-05
[19,] 0.9997769 4.461947e-04 2.230973e-04
[20,] 0.9994935 1.013058e-03 5.065290e-04
[21,] 0.9992691 1.461706e-03 7.308531e-04
[22,] 0.9979589 4.082142e-03 2.041071e-03
[23,] 0.9931974 1.360511e-02 6.802553e-03
[24,] 0.9918414 1.631716e-02 8.158580e-03
[25,] 0.9787181 4.256381e-02 2.128191e-02
[26,] 0.9858075 2.838507e-02 1.419253e-02
> postscript(file="/var/www/html/rcomp/tmp/1yxsp1292950901.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/www/html/rcomp/tmp/2yxsp1292950901.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/www/html/rcomp/tmp/3r69s1292950901.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/www/html/rcomp/tmp/4r69s1292950901.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/www/html/rcomp/tmp/5r69s1292950901.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 = 59
Frequency = 1
1 2 3 4 5
5.162038000 1.851835979 4.843521720 4.639690317 5.383647776
6 7 8 9 10
4.655591054 2.533259120 -0.482570714 -4.417213572 -5.804054434
11 12 13 14 15
-8.632662626 -11.119307803 -15.751459900 -5.475528388 -2.997340546
16 17 18 19 20
-4.190286977 -5.575728156 -5.262017548 -5.394879942 -4.083713975
21 22 23 24 25
-1.419275950 1.286674680 5.763791277 8.484919212 9.776440964
26 27 28 29 30
9.086945351 -3.539303915 -3.095643200 -0.758894250 -2.719106957
31 32 33 34 35
-0.022216202 -0.990343827 1.730278545 0.007959746 -0.012520821
36 37 38 39 40
-0.704326183 -1.801354854 -3.722651940 0.657938737 -0.628718302
41 42 43 44 45
-0.602145300 -0.747586479 1.045417720 1.783014884 1.813813205
46 47 48 49 50
1.842873079 2.874335790 3.338714775 2.614335790 -1.740601001
51 52 53 54 55
1.035184004 3.274958162 1.553119930 4.073119930 1.838419304
56 57 58 59
3.773613632 2.292397772 2.666546928 0.007056380
> postscript(file="/var/www/html/rcomp/tmp/6r69s1292950901.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 = 59
Frequency = 1
lag(myerror, k = 1) myerror
0 5.162038000 NA
1 1.851835979 5.162038000
2 4.843521720 1.851835979
3 4.639690317 4.843521720
4 5.383647776 4.639690317
5 4.655591054 5.383647776
6 2.533259120 4.655591054
7 -0.482570714 2.533259120
8 -4.417213572 -0.482570714
9 -5.804054434 -4.417213572
10 -8.632662626 -5.804054434
11 -11.119307803 -8.632662626
12 -15.751459900 -11.119307803
13 -5.475528388 -15.751459900
14 -2.997340546 -5.475528388
15 -4.190286977 -2.997340546
16 -5.575728156 -4.190286977
17 -5.262017548 -5.575728156
18 -5.394879942 -5.262017548
19 -4.083713975 -5.394879942
20 -1.419275950 -4.083713975
21 1.286674680 -1.419275950
22 5.763791277 1.286674680
23 8.484919212 5.763791277
24 9.776440964 8.484919212
25 9.086945351 9.776440964
26 -3.539303915 9.086945351
27 -3.095643200 -3.539303915
28 -0.758894250 -3.095643200
29 -2.719106957 -0.758894250
30 -0.022216202 -2.719106957
31 -0.990343827 -0.022216202
32 1.730278545 -0.990343827
33 0.007959746 1.730278545
34 -0.012520821 0.007959746
35 -0.704326183 -0.012520821
36 -1.801354854 -0.704326183
37 -3.722651940 -1.801354854
38 0.657938737 -3.722651940
39 -0.628718302 0.657938737
40 -0.602145300 -0.628718302
41 -0.747586479 -0.602145300
42 1.045417720 -0.747586479
43 1.783014884 1.045417720
44 1.813813205 1.783014884
45 1.842873079 1.813813205
46 2.874335790 1.842873079
47 3.338714775 2.874335790
48 2.614335790 3.338714775
49 -1.740601001 2.614335790
50 1.035184004 -1.740601001
51 3.274958162 1.035184004
52 1.553119930 3.274958162
53 4.073119930 1.553119930
54 1.838419304 4.073119930
55 3.773613632 1.838419304
56 2.292397772 3.773613632
57 2.666546928 2.292397772
58 0.007056380 2.666546928
59 NA 0.007056380
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1.851835979 5.162038000
[2,] 4.843521720 1.851835979
[3,] 4.639690317 4.843521720
[4,] 5.383647776 4.639690317
[5,] 4.655591054 5.383647776
[6,] 2.533259120 4.655591054
[7,] -0.482570714 2.533259120
[8,] -4.417213572 -0.482570714
[9,] -5.804054434 -4.417213572
[10,] -8.632662626 -5.804054434
[11,] -11.119307803 -8.632662626
[12,] -15.751459900 -11.119307803
[13,] -5.475528388 -15.751459900
[14,] -2.997340546 -5.475528388
[15,] -4.190286977 -2.997340546
[16,] -5.575728156 -4.190286977
[17,] -5.262017548 -5.575728156
[18,] -5.394879942 -5.262017548
[19,] -4.083713975 -5.394879942
[20,] -1.419275950 -4.083713975
[21,] 1.286674680 -1.419275950
[22,] 5.763791277 1.286674680
[23,] 8.484919212 5.763791277
[24,] 9.776440964 8.484919212
[25,] 9.086945351 9.776440964
[26,] -3.539303915 9.086945351
[27,] -3.095643200 -3.539303915
[28,] -0.758894250 -3.095643200
[29,] -2.719106957 -0.758894250
[30,] -0.022216202 -2.719106957
[31,] -0.990343827 -0.022216202
[32,] 1.730278545 -0.990343827
[33,] 0.007959746 1.730278545
[34,] -0.012520821 0.007959746
[35,] -0.704326183 -0.012520821
[36,] -1.801354854 -0.704326183
[37,] -3.722651940 -1.801354854
[38,] 0.657938737 -3.722651940
[39,] -0.628718302 0.657938737
[40,] -0.602145300 -0.628718302
[41,] -0.747586479 -0.602145300
[42,] 1.045417720 -0.747586479
[43,] 1.783014884 1.045417720
[44,] 1.813813205 1.783014884
[45,] 1.842873079 1.813813205
[46,] 2.874335790 1.842873079
[47,] 3.338714775 2.874335790
[48,] 2.614335790 3.338714775
[49,] -1.740601001 2.614335790
[50,] 1.035184004 -1.740601001
[51,] 3.274958162 1.035184004
[52,] 1.553119930 3.274958162
[53,] 4.073119930 1.553119930
[54,] 1.838419304 4.073119930
[55,] 3.773613632 1.838419304
[56,] 2.292397772 3.773613632
[57,] 2.666546928 2.292397772
[58,] 0.007056380 2.666546928
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1.851835979 5.162038000
2 4.843521720 1.851835979
3 4.639690317 4.843521720
4 5.383647776 4.639690317
5 4.655591054 5.383647776
6 2.533259120 4.655591054
7 -0.482570714 2.533259120
8 -4.417213572 -0.482570714
9 -5.804054434 -4.417213572
10 -8.632662626 -5.804054434
11 -11.119307803 -8.632662626
12 -15.751459900 -11.119307803
13 -5.475528388 -15.751459900
14 -2.997340546 -5.475528388
15 -4.190286977 -2.997340546
16 -5.575728156 -4.190286977
17 -5.262017548 -5.575728156
18 -5.394879942 -5.262017548
19 -4.083713975 -5.394879942
20 -1.419275950 -4.083713975
21 1.286674680 -1.419275950
22 5.763791277 1.286674680
23 8.484919212 5.763791277
24 9.776440964 8.484919212
25 9.086945351 9.776440964
26 -3.539303915 9.086945351
27 -3.095643200 -3.539303915
28 -0.758894250 -3.095643200
29 -2.719106957 -0.758894250
30 -0.022216202 -2.719106957
31 -0.990343827 -0.022216202
32 1.730278545 -0.990343827
33 0.007959746 1.730278545
34 -0.012520821 0.007959746
35 -0.704326183 -0.012520821
36 -1.801354854 -0.704326183
37 -3.722651940 -1.801354854
38 0.657938737 -3.722651940
39 -0.628718302 0.657938737
40 -0.602145300 -0.628718302
41 -0.747586479 -0.602145300
42 1.045417720 -0.747586479
43 1.783014884 1.045417720
44 1.813813205 1.783014884
45 1.842873079 1.813813205
46 2.874335790 1.842873079
47 3.338714775 2.874335790
48 2.614335790 3.338714775
49 -1.740601001 2.614335790
50 1.035184004 -1.740601001
51 3.274958162 1.035184004
52 1.553119930 3.274958162
53 4.073119930 1.553119930
54 1.838419304 4.073119930
55 3.773613632 1.838419304
56 2.292397772 3.773613632
57 2.666546928 2.292397772
58 0.007056380 2.666546928
> 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/71yrv1292950901.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/www/html/rcomp/tmp/8c7qg1292950901.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/www/html/rcomp/tmp/9c7qg1292950901.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/www/html/rcomp/tmp/10c7qg1292950901.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/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/11qz671292950901.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/121q5a1292950901.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/13qr2l1292950901.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/1400j61292950901.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/15mj0u1292950901.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/16isfl1292950901.tab")
+ }
>
> try(system("convert tmp/1yxsp1292950901.ps tmp/1yxsp1292950901.png",intern=TRUE))
character(0)
> try(system("convert tmp/2yxsp1292950901.ps tmp/2yxsp1292950901.png",intern=TRUE))
character(0)
> try(system("convert tmp/3r69s1292950901.ps tmp/3r69s1292950901.png",intern=TRUE))
character(0)
> try(system("convert tmp/4r69s1292950901.ps tmp/4r69s1292950901.png",intern=TRUE))
character(0)
> try(system("convert tmp/5r69s1292950901.ps tmp/5r69s1292950901.png",intern=TRUE))
character(0)
> try(system("convert tmp/6r69s1292950901.ps tmp/6r69s1292950901.png",intern=TRUE))
character(0)
> try(system("convert tmp/71yrv1292950901.ps tmp/71yrv1292950901.png",intern=TRUE))
character(0)
> try(system("convert tmp/8c7qg1292950901.ps tmp/8c7qg1292950901.png",intern=TRUE))
character(0)
> try(system("convert tmp/9c7qg1292950901.ps tmp/9c7qg1292950901.png",intern=TRUE))
character(0)
> try(system("convert tmp/10c7qg1292950901.ps tmp/10c7qg1292950901.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.456 1.655 5.982