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(89.1,72.7,82.6,79.7,102.7,115.8,91.8,87.8,94.1,99.2,103.1,111.4,93.2,102.3,91,94.4,94.3,118.5,99.4,112.1,115.7,136.5,116.8,139.8,99.8,104.5,96,123.3,115.9,156.6,109.1,136.2,117.3,147.5,109.8,143.8,112.8,135.8,110.7,121.6,100,128,113.3,129.7,122.4,136.2,112.5,130.5,104.2,99.2,92.5,110.4,117.2,151.6,109.3,129.6,106.1,123.6,118.8,142.7,105.3,119,106,118.1,102,120,112.9,124.3,116.5,123.3,114.8,122.4,100.5,90.5,85.4,91,114.6,137,109.9,127.7,100.7,105.1,115.5,135.6,100.7,112.4,99,102.5,102.3,112.6,108.8,110.8,105.9,103.4,113.2,117.6,95.7,87.5,80.9,87,113.9,130,98.1,102.9,102.8,111.1,104.7,128.9,95.9,106.3,94.6,99,101.6,109.9,103.9,104,110.3,112.9,114.1,113.6),dim=c(2,60),dimnames=list(c('TotaleIndustrieleProductie','Investeringsgoederen'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('TotaleIndustrieleProductie','Investeringsgoederen'),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
TotaleIndustrieleProductie Investeringsgoederen M1 M2 M3 M4 M5 M6 M7 M8 M9
1 89.1 72.7 1 0 0 0 0 0 0 0 0
2 82.6 79.7 0 1 0 0 0 0 0 0 0
3 102.7 115.8 0 0 1 0 0 0 0 0 0
4 91.8 87.8 0 0 0 1 0 0 0 0 0
5 94.1 99.2 0 0 0 0 1 0 0 0 0
6 103.1 111.4 0 0 0 0 0 1 0 0 0
7 93.2 102.3 0 0 0 0 0 0 1 0 0
8 91.0 94.4 0 0 0 0 0 0 0 1 0
9 94.3 118.5 0 0 0 0 0 0 0 0 1
10 99.4 112.1 0 0 0 0 0 0 0 0 0
11 115.7 136.5 0 0 0 0 0 0 0 0 0
12 116.8 139.8 0 0 0 0 0 0 0 0 0
13 99.8 104.5 1 0 0 0 0 0 0 0 0
14 96.0 123.3 0 1 0 0 0 0 0 0 0
15 115.9 156.6 0 0 1 0 0 0 0 0 0
16 109.1 136.2 0 0 0 1 0 0 0 0 0
17 117.3 147.5 0 0 0 0 1 0 0 0 0
18 109.8 143.8 0 0 0 0 0 1 0 0 0
19 112.8 135.8 0 0 0 0 0 0 1 0 0
20 110.7 121.6 0 0 0 0 0 0 0 1 0
21 100.0 128.0 0 0 0 0 0 0 0 0 1
22 113.3 129.7 0 0 0 0 0 0 0 0 0
23 122.4 136.2 0 0 0 0 0 0 0 0 0
24 112.5 130.5 0 0 0 0 0 0 0 0 0
25 104.2 99.2 1 0 0 0 0 0 0 0 0
26 92.5 110.4 0 1 0 0 0 0 0 0 0
27 117.2 151.6 0 0 1 0 0 0 0 0 0
28 109.3 129.6 0 0 0 1 0 0 0 0 0
29 106.1 123.6 0 0 0 0 1 0 0 0 0
30 118.8 142.7 0 0 0 0 0 1 0 0 0
31 105.3 119.0 0 0 0 0 0 0 1 0 0
32 106.0 118.1 0 0 0 0 0 0 0 1 0
33 102.0 120.0 0 0 0 0 0 0 0 0 1
34 112.9 124.3 0 0 0 0 0 0 0 0 0
35 116.5 123.3 0 0 0 0 0 0 0 0 0
36 114.8 122.4 0 0 0 0 0 0 0 0 0
37 100.5 90.5 1 0 0 0 0 0 0 0 0
38 85.4 91.0 0 1 0 0 0 0 0 0 0
39 114.6 137.0 0 0 1 0 0 0 0 0 0
40 109.9 127.7 0 0 0 1 0 0 0 0 0
41 100.7 105.1 0 0 0 0 1 0 0 0 0
42 115.5 135.6 0 0 0 0 0 1 0 0 0
43 100.7 112.4 0 0 0 0 0 0 1 0 0
44 99.0 102.5 0 0 0 0 0 0 0 1 0
45 102.3 112.6 0 0 0 0 0 0 0 0 1
46 108.8 110.8 0 0 0 0 0 0 0 0 0
47 105.9 103.4 0 0 0 0 0 0 0 0 0
48 113.2 117.6 0 0 0 0 0 0 0 0 0
49 95.7 87.5 1 0 0 0 0 0 0 0 0
50 80.9 87.0 0 1 0 0 0 0 0 0 0
51 113.9 130.0 0 0 1 0 0 0 0 0 0
52 98.1 102.9 0 0 0 1 0 0 0 0 0
53 102.8 111.1 0 0 0 0 1 0 0 0 0
54 104.7 128.9 0 0 0 0 0 1 0 0 0
55 95.9 106.3 0 0 0 0 0 0 1 0 0
56 94.6 99.0 0 0 0 0 0 0 0 1 0
57 101.6 109.9 0 0 0 0 0 0 0 0 1
58 103.9 104.0 0 0 0 0 0 0 0 0 0
59 110.3 112.9 0 0 0 0 0 0 0 0 0
60 114.1 113.6 0 0 0 0 0 0 0 0 0
M10 M11 t
1 0 0 1
2 0 0 2
3 0 0 3
4 0 0 4
5 0 0 5
6 0 0 6
7 0 0 7
8 0 0 8
9 0 0 9
10 1 0 10
11 0 1 11
12 0 0 12
13 0 0 13
14 0 0 14
15 0 0 15
16 0 0 16
17 0 0 17
18 0 0 18
19 0 0 19
20 0 0 20
21 0 0 21
22 1 0 22
23 0 1 23
24 0 0 24
25 0 0 25
26 0 0 26
27 0 0 27
28 0 0 28
29 0 0 29
30 0 0 30
31 0 0 31
32 0 0 32
33 0 0 33
34 1 0 34
35 0 1 35
36 0 0 36
37 0 0 37
38 0 0 38
39 0 0 39
40 0 0 40
41 0 0 41
42 0 0 42
43 0 0 43
44 0 0 44
45 0 0 45
46 1 0 46
47 0 1 47
48 0 0 48
49 0 0 49
50 0 0 50
51 0 0 51
52 0 0 52
53 0 0 53
54 0 0 54
55 0 0 55
56 0 0 56
57 0 0 57
58 1 0 58
59 0 1 59
60 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Investeringsgoederen M1
60.57233 0.40941 -1.73997
M2 M3 M4
-15.22241 -6.25875 -6.80666
M5 M6 M7
-6.50782 -6.61544 -8.39735
M8 M9 M10
-6.49856 -11.16385 -2.95344
M11 t
0.90265 0.07283
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.9623 -1.7409 0.2445 1.9402 5.3858
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 60.57233 4.25941 14.221 < 2e-16 ***
Investeringsgoederen 0.40941 0.03049 13.429 < 2e-16 ***
M1 -1.73997 2.19911 -0.791 0.432877
M2 -15.22241 2.09314 -7.273 3.56e-09 ***
M3 -6.25875 1.94557 -3.217 0.002374 **
M4 -6.80666 1.92751 -3.531 0.000953 ***
M5 -6.50782 1.92294 -3.384 0.001468 **
M6 -6.61544 1.91523 -3.454 0.001197 **
M7 -8.39735 1.92811 -4.355 7.36e-05 ***
M8 -6.49856 1.98052 -3.281 0.001976 **
M9 -11.16385 1.91342 -5.834 5.12e-07 ***
M10 -2.95344 1.91859 -1.539 0.130562
M11 0.90265 1.90046 0.475 0.637058
t 0.07283 0.02324 3.135 0.002996 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.002 on 46 degrees of freedom
Multiple R-squared: 0.9241, Adjusted R-squared: 0.9027
F-statistic: 43.11 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.4123068 0.82461358 0.58769321
[2,] 0.3385704 0.67714084 0.66142958
[3,] 0.7756423 0.44871548 0.22435774
[4,] 0.9120712 0.17585768 0.08792884
[5,] 0.9345969 0.13080628 0.06540314
[6,] 0.9203536 0.15929289 0.07964645
[7,] 0.8769554 0.24608911 0.12304455
[8,] 0.9675506 0.06489875 0.03244938
[9,] 0.9458107 0.10837867 0.05418934
[10,] 0.9400039 0.11999211 0.05999605
[11,] 0.9532825 0.09343500 0.04671750
[12,] 0.9291540 0.14169209 0.07084605
[13,] 0.9300406 0.13991889 0.06995944
[14,] 0.9339293 0.13214141 0.06607070
[15,] 0.9084260 0.18314795 0.09157398
[16,] 0.8653167 0.26936660 0.13468330
[17,] 0.9028572 0.19428564 0.09714282
[18,] 0.8778167 0.24436655 0.12218328
[19,] 0.8170512 0.36589761 0.18294881
[20,] 0.8217968 0.35640648 0.17820324
[21,] 0.7699964 0.46000719 0.23000360
[22,] 0.7309058 0.53818842 0.26909421
[23,] 0.7222830 0.55543399 0.27771700
[24,] 0.6799394 0.64012114 0.32006057
[25,] 0.5577152 0.88456950 0.44228475
[26,] 0.8306218 0.33875643 0.16937822
[27,] 0.7308661 0.53826788 0.26913394
> postscript(file="/var/www/html/rcomp/tmp/1wyjx1258721972.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/2ob6k1258721972.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/3j7c41258721972.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/4cl571258721972.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/5qglz1258721972.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 60
Frequency = 1
1 2 3 4 5 6
0.43096483 4.47472743 0.75866022 1.79711933 -0.94178769 3.09824366
7 8 9 10 11 12
-1.36707770 -2.30439461 -4.27862972 -4.84166669 -2.46011232 -1.88133068
13 14 15 16 17 18
-2.76213900 -0.84937300 -3.61910203 -1.59213226 1.60990138 -4.34050407
19 20 21 22 23 24
3.64382743 5.38577142 -3.34196880 0.97880165 3.48873222 -3.24782772
25 26 27 28 29 30
2.93373799 0.05799333 -1.14604698 0.43597317 -0.67926089 4.23586566
31 32 33 34 35 36
2.14787907 1.24471673 1.05930572 1.91561929 1.99609855 1.49438745
37 38 39 40 41 42
1.92159706 0.02650185 1.35731038 0.93986810 0.62078178 2.96867434
43 44 45 46 47 48
-0.62401549 -0.24251942 3.51493634 2.46862951 -1.33068969 0.98556120
49 50 51 52 53 54
-2.52416088 -3.70984960 2.64917841 -1.58082833 -0.60963458 -5.96227958
55 56 57 58 59 60
-3.80061331 -4.08357411 3.04635646 -0.52138376 -1.69402877 2.64920975
> postscript(file="/var/www/html/rcomp/tmp/6lan61258721972.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 0.43096483 NA
1 4.47472743 0.43096483
2 0.75866022 4.47472743
3 1.79711933 0.75866022
4 -0.94178769 1.79711933
5 3.09824366 -0.94178769
6 -1.36707770 3.09824366
7 -2.30439461 -1.36707770
8 -4.27862972 -2.30439461
9 -4.84166669 -4.27862972
10 -2.46011232 -4.84166669
11 -1.88133068 -2.46011232
12 -2.76213900 -1.88133068
13 -0.84937300 -2.76213900
14 -3.61910203 -0.84937300
15 -1.59213226 -3.61910203
16 1.60990138 -1.59213226
17 -4.34050407 1.60990138
18 3.64382743 -4.34050407
19 5.38577142 3.64382743
20 -3.34196880 5.38577142
21 0.97880165 -3.34196880
22 3.48873222 0.97880165
23 -3.24782772 3.48873222
24 2.93373799 -3.24782772
25 0.05799333 2.93373799
26 -1.14604698 0.05799333
27 0.43597317 -1.14604698
28 -0.67926089 0.43597317
29 4.23586566 -0.67926089
30 2.14787907 4.23586566
31 1.24471673 2.14787907
32 1.05930572 1.24471673
33 1.91561929 1.05930572
34 1.99609855 1.91561929
35 1.49438745 1.99609855
36 1.92159706 1.49438745
37 0.02650185 1.92159706
38 1.35731038 0.02650185
39 0.93986810 1.35731038
40 0.62078178 0.93986810
41 2.96867434 0.62078178
42 -0.62401549 2.96867434
43 -0.24251942 -0.62401549
44 3.51493634 -0.24251942
45 2.46862951 3.51493634
46 -1.33068969 2.46862951
47 0.98556120 -1.33068969
48 -2.52416088 0.98556120
49 -3.70984960 -2.52416088
50 2.64917841 -3.70984960
51 -1.58082833 2.64917841
52 -0.60963458 -1.58082833
53 -5.96227958 -0.60963458
54 -3.80061331 -5.96227958
55 -4.08357411 -3.80061331
56 3.04635646 -4.08357411
57 -0.52138376 3.04635646
58 -1.69402877 -0.52138376
59 2.64920975 -1.69402877
60 NA 2.64920975
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4.47472743 0.43096483
[2,] 0.75866022 4.47472743
[3,] 1.79711933 0.75866022
[4,] -0.94178769 1.79711933
[5,] 3.09824366 -0.94178769
[6,] -1.36707770 3.09824366
[7,] -2.30439461 -1.36707770
[8,] -4.27862972 -2.30439461
[9,] -4.84166669 -4.27862972
[10,] -2.46011232 -4.84166669
[11,] -1.88133068 -2.46011232
[12,] -2.76213900 -1.88133068
[13,] -0.84937300 -2.76213900
[14,] -3.61910203 -0.84937300
[15,] -1.59213226 -3.61910203
[16,] 1.60990138 -1.59213226
[17,] -4.34050407 1.60990138
[18,] 3.64382743 -4.34050407
[19,] 5.38577142 3.64382743
[20,] -3.34196880 5.38577142
[21,] 0.97880165 -3.34196880
[22,] 3.48873222 0.97880165
[23,] -3.24782772 3.48873222
[24,] 2.93373799 -3.24782772
[25,] 0.05799333 2.93373799
[26,] -1.14604698 0.05799333
[27,] 0.43597317 -1.14604698
[28,] -0.67926089 0.43597317
[29,] 4.23586566 -0.67926089
[30,] 2.14787907 4.23586566
[31,] 1.24471673 2.14787907
[32,] 1.05930572 1.24471673
[33,] 1.91561929 1.05930572
[34,] 1.99609855 1.91561929
[35,] 1.49438745 1.99609855
[36,] 1.92159706 1.49438745
[37,] 0.02650185 1.92159706
[38,] 1.35731038 0.02650185
[39,] 0.93986810 1.35731038
[40,] 0.62078178 0.93986810
[41,] 2.96867434 0.62078178
[42,] -0.62401549 2.96867434
[43,] -0.24251942 -0.62401549
[44,] 3.51493634 -0.24251942
[45,] 2.46862951 3.51493634
[46,] -1.33068969 2.46862951
[47,] 0.98556120 -1.33068969
[48,] -2.52416088 0.98556120
[49,] -3.70984960 -2.52416088
[50,] 2.64917841 -3.70984960
[51,] -1.58082833 2.64917841
[52,] -0.60963458 -1.58082833
[53,] -5.96227958 -0.60963458
[54,] -3.80061331 -5.96227958
[55,] -4.08357411 -3.80061331
[56,] 3.04635646 -4.08357411
[57,] -0.52138376 3.04635646
[58,] -1.69402877 -0.52138376
[59,] 2.64920975 -1.69402877
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4.47472743 0.43096483
2 0.75866022 4.47472743
3 1.79711933 0.75866022
4 -0.94178769 1.79711933
5 3.09824366 -0.94178769
6 -1.36707770 3.09824366
7 -2.30439461 -1.36707770
8 -4.27862972 -2.30439461
9 -4.84166669 -4.27862972
10 -2.46011232 -4.84166669
11 -1.88133068 -2.46011232
12 -2.76213900 -1.88133068
13 -0.84937300 -2.76213900
14 -3.61910203 -0.84937300
15 -1.59213226 -3.61910203
16 1.60990138 -1.59213226
17 -4.34050407 1.60990138
18 3.64382743 -4.34050407
19 5.38577142 3.64382743
20 -3.34196880 5.38577142
21 0.97880165 -3.34196880
22 3.48873222 0.97880165
23 -3.24782772 3.48873222
24 2.93373799 -3.24782772
25 0.05799333 2.93373799
26 -1.14604698 0.05799333
27 0.43597317 -1.14604698
28 -0.67926089 0.43597317
29 4.23586566 -0.67926089
30 2.14787907 4.23586566
31 1.24471673 2.14787907
32 1.05930572 1.24471673
33 1.91561929 1.05930572
34 1.99609855 1.91561929
35 1.49438745 1.99609855
36 1.92159706 1.49438745
37 0.02650185 1.92159706
38 1.35731038 0.02650185
39 0.93986810 1.35731038
40 0.62078178 0.93986810
41 2.96867434 0.62078178
42 -0.62401549 2.96867434
43 -0.24251942 -0.62401549
44 3.51493634 -0.24251942
45 2.46862951 3.51493634
46 -1.33068969 2.46862951
47 0.98556120 -1.33068969
48 -2.52416088 0.98556120
49 -3.70984960 -2.52416088
50 2.64917841 -3.70984960
51 -1.58082833 2.64917841
52 -0.60963458 -1.58082833
53 -5.96227958 -0.60963458
54 -3.80061331 -5.96227958
55 -4.08357411 -3.80061331
56 3.04635646 -4.08357411
57 -0.52138376 3.04635646
58 -1.69402877 -0.52138376
59 2.64920975 -1.69402877
> 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/7pd8v1258721972.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/8t88m1258721972.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/9kn7e1258721972.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/10t5wu1258721972.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/11inz01258721972.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/123ouw1258721972.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/132e2e1258721972.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/14rlli1258721972.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/15ij5h1258721972.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/16qric1258721972.tab")
+ }
>
> system("convert tmp/1wyjx1258721972.ps tmp/1wyjx1258721972.png")
> system("convert tmp/2ob6k1258721972.ps tmp/2ob6k1258721972.png")
> system("convert tmp/3j7c41258721972.ps tmp/3j7c41258721972.png")
> system("convert tmp/4cl571258721972.ps tmp/4cl571258721972.png")
> system("convert tmp/5qglz1258721972.ps tmp/5qglz1258721972.png")
> system("convert tmp/6lan61258721972.ps tmp/6lan61258721972.png")
> system("convert tmp/7pd8v1258721972.ps tmp/7pd8v1258721972.png")
> system("convert tmp/8t88m1258721972.ps tmp/8t88m1258721972.png")
> system("convert tmp/9kn7e1258721972.ps tmp/9kn7e1258721972.png")
> system("convert tmp/10t5wu1258721972.ps tmp/10t5wu1258721972.png")
>
>
> proc.time()
user system elapsed
2.369 1.568 2.767