R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(117.6
+ ,112.3
+ ,103.9
+ ,100.4
+ ,118.6
+ ,125.7
+ ,110.9
+ ,108.6
+ ,110.1
+ ,117.5
+ ,106.6
+ ,109.9
+ ,113.1
+ ,105.8
+ ,103.2
+ ,109.5
+ ,91.2
+ ,115.1
+ ,108.7
+ ,109.2
+ ,106.2
+ ,119.3
+ ,104.6
+ ,100.9
+ ,115.5
+ ,119.0
+ ,108.7
+ ,98.4
+ ,106.2
+ ,126.8
+ ,109.5
+ ,94.2
+ ,95.9
+ ,116.3
+ ,102.6
+ ,94.7
+ ,113.2
+ ,127.1
+ ,109.5
+ ,95.2
+ ,78.3
+ ,106.5
+ ,108.1
+ ,100.3
+ ,79.8
+ ,95.5
+ ,96.1
+ ,100.9
+ ,121.2
+ ,121.3
+ ,118.5
+ ,97.9
+ ,125.6
+ ,118.3
+ ,116.3
+ ,106.9
+ ,97.2
+ ,95.6
+ ,105.9
+ ,100.8
+ ,102.8
+ ,90.1
+ ,120.7
+ ,106.6
+ ,88.8
+ ,82.6
+ ,97.0
+ ,108.2
+ ,95.3
+ ,86.9
+ ,99.6
+ ,98.4
+ ,107.6
+ ,95.9
+ ,114.2
+ ,102.0
+ ,95.0
+ ,88.8
+ ,103.3
+ ,95.7
+ ,87.5
+ ,93.2
+ ,99.6
+ ,100.8
+ ,106.7
+ ,98.9
+ ,110.1
+ ,98.8
+ ,75.8
+ ,79.6
+ ,105.7
+ ,99.6
+ ,80.0
+ ,80.7
+ ,97.8
+ ,106.1
+ ,117.2
+ ,102.9
+ ,111.1
+ ,106.3
+ ,106.6
+ ,101.7
+ ,112.0
+ ,105.7
+ ,104.7
+ ,95.9
+ ,107.2
+ ,103.7
+ ,95.2
+ ,87.1
+ ,112.7
+ ,111.2
+ ,94.0
+ ,87.5
+ ,102.5
+ ,114.8
+ ,95.7
+ ,93.6
+ ,114.9
+ ,103.6
+ ,112.6
+ ,109.6
+ ,126.4
+ ,107.0
+ ,99.1
+ ,103.2
+ ,107.3
+ ,104.8
+ ,91.6
+ ,98.9
+ ,103.8
+ ,104.7
+ ,111.5
+ ,113.7
+ ,124.5
+ ,102.0
+ ,76.6
+ ,87.9
+ ,110.1
+ ,103.4
+ ,83.4
+ ,89.6
+ ,107.1
+ ,107.0
+ ,113.5
+ ,114.4
+ ,117.3
+ ,104.0
+ ,106.4
+ ,108.8
+ ,113.8
+ ,105.4
+ ,104.1
+ ,104.3
+ ,119.0
+ ,107.9
+ ,108.4
+ ,97.0
+ ,119.1
+ ,110.1
+ ,91.0
+ ,100.4
+ ,106.9
+ ,111.0
+ ,108.3
+ ,105.3
+ ,115.0
+ ,98.5
+ ,121.0
+ ,122.3
+ ,141.7
+ ,101.9
+ ,95.4
+ ,105.6
+ ,115.2
+ ,103.4
+ ,109.9
+ ,116.9
+ ,117.9
+ ,102.9
+ ,101.4
+ ,110.5
+ ,116.5
+ ,101.0
+ ,86.0
+ ,88.6
+ ,107.4
+ ,103.4
+ ,96.5
+ ,94.5
+ ,122.2
+ ,107.2
+ ,124.6
+ ,115.7
+ ,126.9
+ ,104.5
+ ,109.3
+ ,107.8
+ ,117.7
+ ,104.7
+ ,104.5
+ ,106.6
+ ,114.4
+ ,107.0
+ ,101.8
+ ,100.7
+ ,113.6
+ ,110.3
+ ,101.5
+ ,100.4
+ ,107.0
+ ,107.9
+ ,103.4
+ ,101.7
+ ,107.5
+ ,97.1
+ ,125.9
+ ,115.2
+ ,121.4
+ ,98.6
+ ,96.8
+ ,100.9
+ ,101.3
+ ,95.3
+ ,104.4
+ ,105.3
+ ,102.9
+ ,101.7
+ ,121.1
+ ,109.8
+ ,109.5
+ ,96.3
+ ,83.7
+ ,92.1
+ ,110.2
+ ,99.0
+ ,91.5
+ ,92.5
+ ,118.2
+ ,104.0)
+ ,dim=c(4
+ ,60)
+ ,dimnames=list(c('durable_consumer_goods'
+ ,'intermediate_and_capital_goods'
+ ,'non-durable_consumer_goods'
+ ,'energy')
+ ,1:60))
> y <- array(NA,dim=c(4,60),dimnames=list(c('durable_consumer_goods','intermediate_and_capital_goods','non-durable_consumer_goods','energy'),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 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '4'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
energy durable_consumer_goods intermediate_and_capital_goods
1 100.4 117.6 112.3
2 108.6 118.6 125.7
3 109.9 110.1 117.5
4 109.5 113.1 105.8
5 109.2 91.2 115.1
6 100.9 106.2 119.3
7 98.4 115.5 119.0
8 94.2 106.2 126.8
9 94.7 95.9 116.3
10 95.2 113.2 127.1
11 100.3 78.3 106.5
12 100.9 79.8 95.5
13 97.9 121.2 121.3
14 106.9 125.6 118.3
15 100.8 97.2 95.6
16 106.6 102.8 90.1
17 108.2 88.8 82.6
18 98.4 95.3 86.9
19 102.0 107.6 95.9
20 95.7 95.0 88.8
21 100.8 87.5 93.2
22 98.8 106.7 98.9
23 99.6 75.8 79.6
24 106.1 80.0 80.7
25 106.3 117.2 102.9
26 105.7 106.6 101.7
27 103.7 104.7 95.9
28 111.2 95.2 87.1
29 114.8 94.0 87.5
30 103.6 95.7 93.6
31 107.0 112.6 109.6
32 104.8 99.1 103.2
33 104.7 91.6 98.9
34 102.0 111.5 113.7
35 103.4 76.6 87.9
36 107.0 83.4 89.6
37 104.0 113.5 114.4
38 105.4 106.4 108.8
39 107.9 104.1 104.3
40 110.1 108.4 97.0
41 111.0 91.0 100.4
42 98.5 108.3 105.3
43 101.9 121.0 122.3
44 103.4 95.4 105.6
45 102.9 109.9 116.9
46 101.0 101.4 110.5
47 103.4 86.0 88.6
48 107.2 96.5 94.5
49 104.5 124.6 115.7
50 104.7 109.3 107.8
51 107.0 104.5 106.6
52 110.3 101.8 100.7
53 107.9 101.5 100.4
54 97.1 103.4 101.7
55 98.6 125.9 115.2
56 95.3 96.8 100.9
57 101.7 104.4 105.3
58 96.3 121.1 109.8
59 99.0 83.7 92.1
60 104.0 91.5 92.5
non-durable_consumer_goods
1 103.9
2 110.9
3 106.6
4 103.2
5 108.7
6 104.6
7 108.7
8 109.5
9 102.6
10 109.5
11 108.1
12 96.1
13 118.5
14 116.3
15 105.9
16 120.7
17 97.0
18 99.6
19 114.2
20 103.3
21 99.6
22 110.1
23 105.7
24 97.8
25 111.1
26 112.0
27 107.2
28 112.7
29 102.5
30 114.9
31 126.4
32 107.3
33 103.8
34 124.5
35 110.1
36 107.1
37 117.3
38 113.8
39 119.0
40 119.1
41 106.9
42 115.0
43 141.7
44 115.2
45 117.9
46 116.5
47 107.4
48 122.2
49 126.9
50 117.7
51 114.4
52 113.6
53 107.0
54 107.5
55 121.4
56 101.3
57 102.9
58 109.5
59 110.2
60 118.2
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) durable_consumer_goods
101.87583 0.04434
intermediate_and_capital_goods `non-durable_consumer_goods`
-0.15587 0.11776
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8.7110 -3.0096 -0.2484 2.4685 10.3249
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 101.87583 8.39451 12.136 <2e-16 ***
durable_consumer_goods 0.04434 0.07360 0.602 0.5493
intermediate_and_capital_goods -0.15587 0.07174 -2.173 0.0341 *
`non-durable_consumer_goods` 0.11776 0.08496 1.386 0.1712
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.615 on 56 degrees of freedom
Multiple R-squared: 0.1082, Adjusted R-squared: 0.06044
F-statistic: 2.265 on 3 and 56 DF, p-value: 0.09084
> 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.9080130 0.1839740 0.09198702
[2,] 0.9392187 0.1215627 0.06078133
[3,] 0.9039247 0.1921506 0.09607532
[4,] 0.8684701 0.2630599 0.13152995
[5,] 0.9286156 0.1427689 0.07138443
[6,] 0.8847087 0.2305827 0.11529134
[7,] 0.9471778 0.1056444 0.05282219
[8,] 0.9283725 0.1432551 0.07162755
[9,] 0.9242058 0.1515884 0.07579420
[10,] 0.8866522 0.2266956 0.11334779
[11,] 0.8555065 0.2889871 0.14449353
[12,] 0.8890817 0.2218365 0.11091826
[13,] 0.8624097 0.2751806 0.13759031
[14,] 0.9316769 0.1366462 0.06832312
[15,] 0.9041900 0.1916199 0.09580997
[16,] 0.9085061 0.1829877 0.09149387
[17,] 0.9240052 0.1519896 0.07599480
[18,] 0.9071887 0.1856226 0.09281132
[19,] 0.8780330 0.2439340 0.12196699
[20,] 0.8425112 0.3149776 0.15748881
[21,] 0.7938699 0.4122602 0.20613009
[22,] 0.8088739 0.3822523 0.19112615
[23,] 0.9442897 0.1114206 0.05571032
[24,] 0.9213493 0.1573013 0.07865067
[25,] 0.8981783 0.2036434 0.10182168
[26,] 0.8669464 0.2661072 0.13305361
[27,] 0.8286552 0.3426896 0.17134482
[28,] 0.7795636 0.4408727 0.22043637
[29,] 0.7349318 0.5301364 0.26506822
[30,] 0.6819753 0.6360495 0.31802473
[31,] 0.6161092 0.7677815 0.38389075
[32,] 0.5644693 0.8710615 0.43553074
[33,] 0.5309030 0.9381940 0.46909700
[34,] 0.5548931 0.8902138 0.44510691
[35,] 0.7359224 0.5281552 0.26407759
[36,] 0.7303321 0.5393359 0.26966794
[37,] 0.7236185 0.5527630 0.27638150
[38,] 0.6380607 0.7238786 0.36193931
[39,] 0.5443025 0.9113949 0.45569747
[40,] 0.4773046 0.9546093 0.52269536
[41,] 0.3789672 0.7579345 0.62103276
[42,] 0.2886558 0.5773116 0.71134421
[43,] 0.2020337 0.4040674 0.79796632
[44,] 0.1310113 0.2620226 0.86898872
[45,] 0.1161462 0.2322924 0.88385378
[46,] 0.2687793 0.5375586 0.73122072
[47,] 0.6182143 0.7635714 0.38178572
> postscript(file="/var/wessaorg/rcomp/tmp/1v76v1353445960.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/2yc9m1353445960.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/35cun1353445960.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/4kgy21353445960.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/5i9w21353445960.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 = 60
Frequency = 1
1 2 3 4 5 6
-1.420757354 7.999278346 8.904377960 6.948041847 8.421024086 0.593382253
7 8 9 10 11 12
-2.848538551 -5.514591682 -5.382017568 -4.778209381 -1.176825985 -0.944842373
13 14 15 16 17 18
-4.396780722 4.199574958 -2.954776114 -0.003154706 3.839387433 -5.884743032
19 20 21 22 23 24
-3.146523024 -8.710983584 -2.156907197 -5.356206148 -5.676282659 1.739219115
25 26 27 28 29 30
2.183952224 1.760929486 -0.493646369 5.408263141 10.324929117 -1.459810503
31 32 33 34 35 36
2.330583962 1.980736873 1.955187465 -1.757836405 -1.136154014 2.780582780
37 38 39 40 41 42
1.110435652 2.364518468 3.652751656 4.512459034 8.150554132 -5.306581540
43 44 45 46 47 48
-2.963978368 0.188612438 0.489082173 -1.866743083 -1.125898217 1.385383330
49 50 51 52 53 54
0.190438863 0.920814375 3.635195222 6.229480683 4.573210206 -6.167281591
55 56 57 58 59 60
-5.197480951 -7.069248737 -0.508810090 -6.725056306 -5.208085766 -1.433635260
> postscript(file="/var/wessaorg/rcomp/tmp/6iin31353445960.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -1.420757354 NA
1 7.999278346 -1.420757354
2 8.904377960 7.999278346
3 6.948041847 8.904377960
4 8.421024086 6.948041847
5 0.593382253 8.421024086
6 -2.848538551 0.593382253
7 -5.514591682 -2.848538551
8 -5.382017568 -5.514591682
9 -4.778209381 -5.382017568
10 -1.176825985 -4.778209381
11 -0.944842373 -1.176825985
12 -4.396780722 -0.944842373
13 4.199574958 -4.396780722
14 -2.954776114 4.199574958
15 -0.003154706 -2.954776114
16 3.839387433 -0.003154706
17 -5.884743032 3.839387433
18 -3.146523024 -5.884743032
19 -8.710983584 -3.146523024
20 -2.156907197 -8.710983584
21 -5.356206148 -2.156907197
22 -5.676282659 -5.356206148
23 1.739219115 -5.676282659
24 2.183952224 1.739219115
25 1.760929486 2.183952224
26 -0.493646369 1.760929486
27 5.408263141 -0.493646369
28 10.324929117 5.408263141
29 -1.459810503 10.324929117
30 2.330583962 -1.459810503
31 1.980736873 2.330583962
32 1.955187465 1.980736873
33 -1.757836405 1.955187465
34 -1.136154014 -1.757836405
35 2.780582780 -1.136154014
36 1.110435652 2.780582780
37 2.364518468 1.110435652
38 3.652751656 2.364518468
39 4.512459034 3.652751656
40 8.150554132 4.512459034
41 -5.306581540 8.150554132
42 -2.963978368 -5.306581540
43 0.188612438 -2.963978368
44 0.489082173 0.188612438
45 -1.866743083 0.489082173
46 -1.125898217 -1.866743083
47 1.385383330 -1.125898217
48 0.190438863 1.385383330
49 0.920814375 0.190438863
50 3.635195222 0.920814375
51 6.229480683 3.635195222
52 4.573210206 6.229480683
53 -6.167281591 4.573210206
54 -5.197480951 -6.167281591
55 -7.069248737 -5.197480951
56 -0.508810090 -7.069248737
57 -6.725056306 -0.508810090
58 -5.208085766 -6.725056306
59 -1.433635260 -5.208085766
60 NA -1.433635260
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 7.999278346 -1.420757354
[2,] 8.904377960 7.999278346
[3,] 6.948041847 8.904377960
[4,] 8.421024086 6.948041847
[5,] 0.593382253 8.421024086
[6,] -2.848538551 0.593382253
[7,] -5.514591682 -2.848538551
[8,] -5.382017568 -5.514591682
[9,] -4.778209381 -5.382017568
[10,] -1.176825985 -4.778209381
[11,] -0.944842373 -1.176825985
[12,] -4.396780722 -0.944842373
[13,] 4.199574958 -4.396780722
[14,] -2.954776114 4.199574958
[15,] -0.003154706 -2.954776114
[16,] 3.839387433 -0.003154706
[17,] -5.884743032 3.839387433
[18,] -3.146523024 -5.884743032
[19,] -8.710983584 -3.146523024
[20,] -2.156907197 -8.710983584
[21,] -5.356206148 -2.156907197
[22,] -5.676282659 -5.356206148
[23,] 1.739219115 -5.676282659
[24,] 2.183952224 1.739219115
[25,] 1.760929486 2.183952224
[26,] -0.493646369 1.760929486
[27,] 5.408263141 -0.493646369
[28,] 10.324929117 5.408263141
[29,] -1.459810503 10.324929117
[30,] 2.330583962 -1.459810503
[31,] 1.980736873 2.330583962
[32,] 1.955187465 1.980736873
[33,] -1.757836405 1.955187465
[34,] -1.136154014 -1.757836405
[35,] 2.780582780 -1.136154014
[36,] 1.110435652 2.780582780
[37,] 2.364518468 1.110435652
[38,] 3.652751656 2.364518468
[39,] 4.512459034 3.652751656
[40,] 8.150554132 4.512459034
[41,] -5.306581540 8.150554132
[42,] -2.963978368 -5.306581540
[43,] 0.188612438 -2.963978368
[44,] 0.489082173 0.188612438
[45,] -1.866743083 0.489082173
[46,] -1.125898217 -1.866743083
[47,] 1.385383330 -1.125898217
[48,] 0.190438863 1.385383330
[49,] 0.920814375 0.190438863
[50,] 3.635195222 0.920814375
[51,] 6.229480683 3.635195222
[52,] 4.573210206 6.229480683
[53,] -6.167281591 4.573210206
[54,] -5.197480951 -6.167281591
[55,] -7.069248737 -5.197480951
[56,] -0.508810090 -7.069248737
[57,] -6.725056306 -0.508810090
[58,] -5.208085766 -6.725056306
[59,] -1.433635260 -5.208085766
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 7.999278346 -1.420757354
2 8.904377960 7.999278346
3 6.948041847 8.904377960
4 8.421024086 6.948041847
5 0.593382253 8.421024086
6 -2.848538551 0.593382253
7 -5.514591682 -2.848538551
8 -5.382017568 -5.514591682
9 -4.778209381 -5.382017568
10 -1.176825985 -4.778209381
11 -0.944842373 -1.176825985
12 -4.396780722 -0.944842373
13 4.199574958 -4.396780722
14 -2.954776114 4.199574958
15 -0.003154706 -2.954776114
16 3.839387433 -0.003154706
17 -5.884743032 3.839387433
18 -3.146523024 -5.884743032
19 -8.710983584 -3.146523024
20 -2.156907197 -8.710983584
21 -5.356206148 -2.156907197
22 -5.676282659 -5.356206148
23 1.739219115 -5.676282659
24 2.183952224 1.739219115
25 1.760929486 2.183952224
26 -0.493646369 1.760929486
27 5.408263141 -0.493646369
28 10.324929117 5.408263141
29 -1.459810503 10.324929117
30 2.330583962 -1.459810503
31 1.980736873 2.330583962
32 1.955187465 1.980736873
33 -1.757836405 1.955187465
34 -1.136154014 -1.757836405
35 2.780582780 -1.136154014
36 1.110435652 2.780582780
37 2.364518468 1.110435652
38 3.652751656 2.364518468
39 4.512459034 3.652751656
40 8.150554132 4.512459034
41 -5.306581540 8.150554132
42 -2.963978368 -5.306581540
43 0.188612438 -2.963978368
44 0.489082173 0.188612438
45 -1.866743083 0.489082173
46 -1.125898217 -1.866743083
47 1.385383330 -1.125898217
48 0.190438863 1.385383330
49 0.920814375 0.190438863
50 3.635195222 0.920814375
51 6.229480683 3.635195222
52 4.573210206 6.229480683
53 -6.167281591 4.573210206
54 -5.197480951 -6.167281591
55 -7.069248737 -5.197480951
56 -0.508810090 -7.069248737
57 -6.725056306 -0.508810090
58 -5.208085766 -6.725056306
59 -1.433635260 -5.208085766
> 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/7ugg11353445960.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/862vh1353445960.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/9k2g61353445960.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/1016t91353445960.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/11xfrg1353445960.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/12r8aw1353445960.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/13taez1353445960.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/144i0w1353445960.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/15jsdv1353445960.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/16e0j21353445960.tab")
+ }
>
> try(system("convert tmp/1v76v1353445960.ps tmp/1v76v1353445960.png",intern=TRUE))
character(0)
> try(system("convert tmp/2yc9m1353445960.ps tmp/2yc9m1353445960.png",intern=TRUE))
character(0)
> try(system("convert tmp/35cun1353445960.ps tmp/35cun1353445960.png",intern=TRUE))
character(0)
> try(system("convert tmp/4kgy21353445960.ps tmp/4kgy21353445960.png",intern=TRUE))
character(0)
> try(system("convert tmp/5i9w21353445960.ps tmp/5i9w21353445960.png",intern=TRUE))
character(0)
> try(system("convert tmp/6iin31353445960.ps tmp/6iin31353445960.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ugg11353445960.ps tmp/7ugg11353445960.png",intern=TRUE))
character(0)
> try(system("convert tmp/862vh1353445960.ps tmp/862vh1353445960.png",intern=TRUE))
character(0)
> try(system("convert tmp/9k2g61353445960.ps tmp/9k2g61353445960.png",intern=TRUE))
character(0)
> try(system("convert tmp/1016t91353445960.ps tmp/1016t91353445960.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.052 0.825 6.878