R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(6
+ ,101.82
+ ,107.34
+ ,93.63
+ ,99.85
+ ,101.76
+ ,6
+ ,101.68
+ ,107.34
+ ,93.63
+ ,99.91
+ ,102.37
+ ,6
+ ,101.68
+ ,107.34
+ ,93.63
+ ,99.87
+ ,102.38
+ ,6
+ ,102.45
+ ,107.34
+ ,96.13
+ ,99.86
+ ,102.86
+ ,6
+ ,102.45
+ ,107.34
+ ,96.13
+ ,100.10
+ ,102.87
+ ,6
+ ,102.45
+ ,107.34
+ ,96.13
+ ,100.10
+ ,102.92
+ ,6
+ ,102.45
+ ,107.34
+ ,96.13
+ ,100.12
+ ,102.95
+ ,6
+ ,102.45
+ ,107.34
+ ,96.13
+ ,99.95
+ ,103.02
+ ,6
+ ,102.45
+ ,112.60
+ ,96.13
+ ,99.94
+ ,104.08
+ ,6
+ ,102.52
+ ,112.60
+ ,96.13
+ ,100.18
+ ,104.16
+ ,6
+ ,102.52
+ ,112.60
+ ,96.13
+ ,100.31
+ ,104.24
+ ,6
+ ,102.85
+ ,112.60
+ ,96.13
+ ,100.65
+ ,104.33
+ ,7
+ ,102.85
+ ,112.61
+ ,96.13
+ ,100.65
+ ,104.73
+ ,7
+ ,102.85
+ ,112.61
+ ,96.13
+ ,100.69
+ ,104.86
+ ,7
+ ,103.25
+ ,112.61
+ ,96.13
+ ,101.26
+ ,105.03
+ ,7
+ ,103.25
+ ,112.61
+ ,98.73
+ ,101.26
+ ,105.62
+ ,7
+ ,103.25
+ ,112.61
+ ,98.73
+ ,101.38
+ ,105.63
+ ,7
+ ,103.25
+ ,112.61
+ ,98.73
+ ,101.38
+ ,105.63
+ ,7
+ ,104.45
+ ,112.61
+ ,98.73
+ ,101.38
+ ,105.94
+ ,7
+ ,104.45
+ ,112.61
+ ,98.73
+ ,101.44
+ ,106.61
+ ,7
+ ,104.45
+ ,118.65
+ ,98.73
+ ,101.40
+ ,107.69
+ ,7
+ ,104.80
+ ,118.65
+ ,98.73
+ ,101.40
+ ,107.78
+ ,7
+ ,104.80
+ ,118.65
+ ,98.73
+ ,100.58
+ ,107.93
+ ,7
+ ,105.29
+ ,118.65
+ ,98.73
+ ,100.58
+ ,108.48
+ ,8
+ ,105.29
+ ,114.29
+ ,98.73
+ ,100.58
+ ,108.14
+ ,8
+ ,105.29
+ ,114.29
+ ,98.73
+ ,100.59
+ ,108.48
+ ,8
+ ,105.29
+ ,114.29
+ ,98.73
+ ,100.81
+ ,108.48
+ ,8
+ ,106.04
+ ,114.29
+ ,101.67
+ ,100.75
+ ,108.89
+ ,8
+ ,105.94
+ ,114.29
+ ,101.67
+ ,100.75
+ ,108.93
+ ,8
+ ,105.94
+ ,114.29
+ ,101.67
+ ,100.96
+ ,109.21
+ ,8
+ ,105.94
+ ,114.29
+ ,101.67
+ ,101.31
+ ,109.47
+ ,8
+ ,106.28
+ ,114.29
+ ,101.67
+ ,101.64
+ ,109.80
+ ,8
+ ,106.48
+ ,123.33
+ ,101.67
+ ,101.46
+ ,111.73
+ ,8
+ ,107.19
+ ,123.33
+ ,101.67
+ ,101.73
+ ,111.85
+ ,8
+ ,108.14
+ ,123.33
+ ,101.67
+ ,101.73
+ ,112.12
+ ,8
+ ,108.22
+ ,123.33
+ ,101.67
+ ,101.64
+ ,112.15
+ ,9
+ ,108.22
+ ,123.33
+ ,101.67
+ ,101.77
+ ,112.17
+ ,9
+ ,108.61
+ ,123.33
+ ,101.67
+ ,101.74
+ ,112.67
+ ,9
+ ,108.61
+ ,123.33
+ ,101.67
+ ,101.89
+ ,112.80
+ ,9
+ ,108.61
+ ,123.33
+ ,107.94
+ ,101.89
+ ,113.44
+ ,9
+ ,108.61
+ ,123.33
+ ,107.94
+ ,101.93
+ ,113.53
+ ,9
+ ,109.06
+ ,123.33
+ ,107.94
+ ,101.93
+ ,114.53
+ ,9
+ ,109.06
+ ,123.33
+ ,107.94
+ ,102.32
+ ,114.51
+ ,9
+ ,112.93
+ ,123.33
+ ,107.94
+ ,102.41
+ ,115.05
+ ,9
+ ,115.84
+ ,129.03
+ ,107.94
+ ,103.58
+ ,116.67
+ ,9
+ ,118.57
+ ,128.76
+ ,107.94
+ ,104.12
+ ,117.07
+ ,9
+ ,118.57
+ ,128.76
+ ,107.94
+ ,104.10
+ ,116.92
+ ,9
+ ,118.86
+ ,128.76
+ ,107.94
+ ,104.15
+ ,117.00
+ ,10
+ ,118.98
+ ,128.76
+ ,107.94
+ ,104.15
+ ,117.02
+ ,10
+ ,119.27
+ ,128.76
+ ,107.94
+ ,104.16
+ ,117.35
+ ,10
+ ,119.39
+ ,128.76
+ ,107.94
+ ,102.94
+ ,117.36
+ ,10
+ ,119.49
+ ,128.76
+ ,110.30
+ ,103.07
+ ,117.82
+ ,10
+ ,119.59
+ ,128.76
+ ,110.30
+ ,103.04
+ ,117.88
+ ,10
+ ,120.12
+ ,128.76
+ ,110.30
+ ,103.06
+ ,118.24
+ ,10
+ ,120.14
+ ,128.76
+ ,110.30
+ ,103.05
+ ,118.50
+ ,10
+ ,120.14
+ ,128.76
+ ,110.30
+ ,102.95
+ ,118.80
+ ,10
+ ,120.14
+ ,132.63
+ ,110.30
+ ,102.95
+ ,119.76
+ ,10
+ ,120.14
+ ,132.63
+ ,110.30
+ ,103.05
+ ,120.09)
+ ,dim=c(6
+ ,58)
+ ,dimnames=list(c('jaar'
+ ,'bioscoop'
+ ,'schouwburgabonnement'
+ ,'eendagsattractie'
+ ,'huurDVD'
+ ,'Cultuurenvrijetijdsbesteding')
+ ,1:58))
> y <- array(NA,dim=c(6,58),dimnames=list(c('jaar','bioscoop','schouwburgabonnement','eendagsattractie','huurDVD','Cultuurenvrijetijdsbesteding'),1:58))
> 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 = '6'
> library(lattice)
> library(lmtest)
Loading required package: zoo
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Cultuurenvrijetijdsbesteding jaar bioscoop schouwburgabonnement
1 101.76 6 101.82 107.34
2 102.37 6 101.68 107.34
3 102.38 6 101.68 107.34
4 102.86 6 102.45 107.34
5 102.87 6 102.45 107.34
6 102.92 6 102.45 107.34
7 102.95 6 102.45 107.34
8 103.02 6 102.45 107.34
9 104.08 6 102.45 112.60
10 104.16 6 102.52 112.60
11 104.24 6 102.52 112.60
12 104.33 6 102.85 112.60
13 104.73 7 102.85 112.61
14 104.86 7 102.85 112.61
15 105.03 7 103.25 112.61
16 105.62 7 103.25 112.61
17 105.63 7 103.25 112.61
18 105.63 7 103.25 112.61
19 105.94 7 104.45 112.61
20 106.61 7 104.45 112.61
21 107.69 7 104.45 118.65
22 107.78 7 104.80 118.65
23 107.93 7 104.80 118.65
24 108.48 7 105.29 118.65
25 108.14 8 105.29 114.29
26 108.48 8 105.29 114.29
27 108.48 8 105.29 114.29
28 108.89 8 106.04 114.29
29 108.93 8 105.94 114.29
30 109.21 8 105.94 114.29
31 109.47 8 105.94 114.29
32 109.80 8 106.28 114.29
33 111.73 8 106.48 123.33
34 111.85 8 107.19 123.33
35 112.12 8 108.14 123.33
36 112.15 8 108.22 123.33
37 112.17 9 108.22 123.33
38 112.67 9 108.61 123.33
39 112.80 9 108.61 123.33
40 113.44 9 108.61 123.33
41 113.53 9 108.61 123.33
42 114.53 9 109.06 123.33
43 114.51 9 109.06 123.33
44 115.05 9 112.93 123.33
45 116.67 9 115.84 129.03
46 117.07 9 118.57 128.76
47 116.92 9 118.57 128.76
48 117.00 9 118.86 128.76
49 117.02 10 118.98 128.76
50 117.35 10 119.27 128.76
51 117.36 10 119.39 128.76
52 117.82 10 119.49 128.76
53 117.88 10 119.59 128.76
54 118.24 10 120.12 128.76
55 118.50 10 120.14 128.76
56 118.80 10 120.14 128.76
57 119.76 10 120.14 132.63
58 120.09 10 120.14 132.63
eendagsattractie huurDVD
1 93.63 99.85
2 93.63 99.91
3 93.63 99.87
4 96.13 99.86
5 96.13 100.10
6 96.13 100.10
7 96.13 100.12
8 96.13 99.95
9 96.13 99.94
10 96.13 100.18
11 96.13 100.31
12 96.13 100.65
13 96.13 100.65
14 96.13 100.69
15 96.13 101.26
16 98.73 101.26
17 98.73 101.38
18 98.73 101.38
19 98.73 101.38
20 98.73 101.44
21 98.73 101.40
22 98.73 101.40
23 98.73 100.58
24 98.73 100.58
25 98.73 100.58
26 98.73 100.59
27 98.73 100.81
28 101.67 100.75
29 101.67 100.75
30 101.67 100.96
31 101.67 101.31
32 101.67 101.64
33 101.67 101.46
34 101.67 101.73
35 101.67 101.73
36 101.67 101.64
37 101.67 101.77
38 101.67 101.74
39 101.67 101.89
40 107.94 101.89
41 107.94 101.93
42 107.94 101.93
43 107.94 102.32
44 107.94 102.41
45 107.94 103.58
46 107.94 104.12
47 107.94 104.10
48 107.94 104.15
49 107.94 104.15
50 107.94 104.16
51 107.94 102.94
52 110.30 103.07
53 110.30 103.04
54 110.30 103.06
55 110.30 103.05
56 110.30 102.95
57 110.30 102.95
58 110.30 103.05
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) jaar bioscoop
43.11965 1.06911 0.09551
schouwburgabonnement eendagsattractie huurDVD
0.29383 0.27101 -0.14029
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.8718 -0.3487 0.0531 0.3250 1.1003
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 43.11965 11.20716 3.848 0.000328 ***
jaar 1.06911 0.14211 7.523 7.26e-10 ***
bioscoop 0.09551 0.03026 3.157 0.002656 **
schouwburgabonnement 0.29383 0.02412 12.181 < 2e-16 ***
eendagsattractie 0.27101 0.04070 6.659 1.73e-08 ***
huurDVD -0.14029 0.13467 -1.042 0.302359
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.453 on 52 degrees of freedom
Multiple R-squared: 0.994, Adjusted R-squared: 0.9934
F-statistic: 1716 on 5 and 52 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.0050207249 0.0100414498 0.99497928
[2,] 0.0066968809 0.0133937618 0.99330312
[3,] 0.0017401582 0.0034803164 0.99825984
[4,] 0.0074195965 0.0148391930 0.99258040
[5,] 0.0027960297 0.0055920594 0.99720397
[6,] 0.0011019431 0.0022038861 0.99889806
[7,] 0.0004766657 0.0009533314 0.99952333
[8,] 0.0006499257 0.0012998514 0.99935007
[9,] 0.0003208096 0.0006416191 0.99967919
[10,] 0.0001728282 0.0003456564 0.99982717
[11,] 0.0001194581 0.0002389161 0.99988054
[12,] 0.0018615106 0.0037230213 0.99813849
[13,] 0.0017155412 0.0034310823 0.99828446
[14,] 0.0017987286 0.0035974573 0.99820127
[15,] 0.0038416732 0.0076833464 0.99615833
[16,] 0.0096527818 0.0193055636 0.99034722
[17,] 0.0079365883 0.0158731766 0.99206341
[18,] 0.0101081416 0.0202162832 0.98989186
[19,] 0.0092300943 0.0184601886 0.99076991
[20,] 0.0122062673 0.0244125346 0.98779373
[21,] 0.0124905485 0.0249810970 0.98750945
[22,] 0.0114283558 0.0228567115 0.98857164
[23,] 0.0190947554 0.0381895109 0.98090524
[24,] 0.0535432595 0.1070865191 0.94645674
[25,] 0.0592902143 0.1185804286 0.94070979
[26,] 0.0416617971 0.0833235942 0.95833820
[27,] 0.0409780146 0.0819560291 0.95902199
[28,] 0.0371994480 0.0743988960 0.96280055
[29,] 0.1469171704 0.2938343408 0.85308283
[30,] 0.1330948611 0.2661897222 0.86690514
[31,] 0.1075052513 0.2150105025 0.89249475
[32,] 0.3926227862 0.7852455725 0.60737721
[33,] 0.7698056739 0.4603886522 0.23019433
[34,] 0.6929113646 0.6141772709 0.30708864
[35,] 0.6016050555 0.7967898890 0.39839494
[36,] 0.8458962986 0.3082074028 0.15410370
[37,] 0.9699615593 0.0600768813 0.03003844
[38,] 0.9594854534 0.0810290933 0.04051455
[39,] 0.9334781827 0.1330436346 0.06652182
[40,] 0.8575499838 0.2849000325 0.14245002
[41,] 0.7774457775 0.4451084450 0.22255422
> postscript(file="/var/wessaorg/rcomp/tmp/1mdpa1321906094.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/20y311321906094.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/3ekaj1321906094.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/4bbtb1321906094.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/5xmdf1321906094.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 = 58
Frequency = 1
1 2 3 4 5 6
-0.405616120 0.226172428 0.230560930 -0.041912378 0.001756608 0.051756608
7 8 9 10 11 12
0.084562356 0.130713492 -0.356222280 -0.249238945 -0.151001578 -0.044821916
13 14 15 16 17 18
-0.716871975 -0.581260477 -0.369500355 -0.484129762 -0.457295269 -0.457295269
19 20 21 22 23 24
-0.261906425 0.416510822 -0.283819059 -0.227247313 -0.192283014 0.310917431
25 26 27 28 29 30
0.182893752 0.524296626 0.555159863 0.088337392 0.137888322 0.447348684
31 32 33 34 35 36
0.756449288 1.100270983 0.329716360 0.419782368 0.599048537 0.608781923
37 38 39 40 41 42
-0.422092493 0.036450258 0.187493374 -0.871747542 -0.776136044 0.180884772
43 44 45 46 47 48
0.215596874 0.398601766 0.229989001 0.524337277 0.371531529 0.430848205
49 50 51 52 53 54
-0.629724694 -0.326019516 -0.498631309 -0.669531564 -0.623291116 -0.311105295
55 56 57 58
-0.054418355 0.231552901 0.054440294 0.398469038
> postscript(file="/var/wessaorg/rcomp/tmp/643l51321906094.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 = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.405616120 NA
1 0.226172428 -0.405616120
2 0.230560930 0.226172428
3 -0.041912378 0.230560930
4 0.001756608 -0.041912378
5 0.051756608 0.001756608
6 0.084562356 0.051756608
7 0.130713492 0.084562356
8 -0.356222280 0.130713492
9 -0.249238945 -0.356222280
10 -0.151001578 -0.249238945
11 -0.044821916 -0.151001578
12 -0.716871975 -0.044821916
13 -0.581260477 -0.716871975
14 -0.369500355 -0.581260477
15 -0.484129762 -0.369500355
16 -0.457295269 -0.484129762
17 -0.457295269 -0.457295269
18 -0.261906425 -0.457295269
19 0.416510822 -0.261906425
20 -0.283819059 0.416510822
21 -0.227247313 -0.283819059
22 -0.192283014 -0.227247313
23 0.310917431 -0.192283014
24 0.182893752 0.310917431
25 0.524296626 0.182893752
26 0.555159863 0.524296626
27 0.088337392 0.555159863
28 0.137888322 0.088337392
29 0.447348684 0.137888322
30 0.756449288 0.447348684
31 1.100270983 0.756449288
32 0.329716360 1.100270983
33 0.419782368 0.329716360
34 0.599048537 0.419782368
35 0.608781923 0.599048537
36 -0.422092493 0.608781923
37 0.036450258 -0.422092493
38 0.187493374 0.036450258
39 -0.871747542 0.187493374
40 -0.776136044 -0.871747542
41 0.180884772 -0.776136044
42 0.215596874 0.180884772
43 0.398601766 0.215596874
44 0.229989001 0.398601766
45 0.524337277 0.229989001
46 0.371531529 0.524337277
47 0.430848205 0.371531529
48 -0.629724694 0.430848205
49 -0.326019516 -0.629724694
50 -0.498631309 -0.326019516
51 -0.669531564 -0.498631309
52 -0.623291116 -0.669531564
53 -0.311105295 -0.623291116
54 -0.054418355 -0.311105295
55 0.231552901 -0.054418355
56 0.054440294 0.231552901
57 0.398469038 0.054440294
58 NA 0.398469038
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.226172428 -0.405616120
[2,] 0.230560930 0.226172428
[3,] -0.041912378 0.230560930
[4,] 0.001756608 -0.041912378
[5,] 0.051756608 0.001756608
[6,] 0.084562356 0.051756608
[7,] 0.130713492 0.084562356
[8,] -0.356222280 0.130713492
[9,] -0.249238945 -0.356222280
[10,] -0.151001578 -0.249238945
[11,] -0.044821916 -0.151001578
[12,] -0.716871975 -0.044821916
[13,] -0.581260477 -0.716871975
[14,] -0.369500355 -0.581260477
[15,] -0.484129762 -0.369500355
[16,] -0.457295269 -0.484129762
[17,] -0.457295269 -0.457295269
[18,] -0.261906425 -0.457295269
[19,] 0.416510822 -0.261906425
[20,] -0.283819059 0.416510822
[21,] -0.227247313 -0.283819059
[22,] -0.192283014 -0.227247313
[23,] 0.310917431 -0.192283014
[24,] 0.182893752 0.310917431
[25,] 0.524296626 0.182893752
[26,] 0.555159863 0.524296626
[27,] 0.088337392 0.555159863
[28,] 0.137888322 0.088337392
[29,] 0.447348684 0.137888322
[30,] 0.756449288 0.447348684
[31,] 1.100270983 0.756449288
[32,] 0.329716360 1.100270983
[33,] 0.419782368 0.329716360
[34,] 0.599048537 0.419782368
[35,] 0.608781923 0.599048537
[36,] -0.422092493 0.608781923
[37,] 0.036450258 -0.422092493
[38,] 0.187493374 0.036450258
[39,] -0.871747542 0.187493374
[40,] -0.776136044 -0.871747542
[41,] 0.180884772 -0.776136044
[42,] 0.215596874 0.180884772
[43,] 0.398601766 0.215596874
[44,] 0.229989001 0.398601766
[45,] 0.524337277 0.229989001
[46,] 0.371531529 0.524337277
[47,] 0.430848205 0.371531529
[48,] -0.629724694 0.430848205
[49,] -0.326019516 -0.629724694
[50,] -0.498631309 -0.326019516
[51,] -0.669531564 -0.498631309
[52,] -0.623291116 -0.669531564
[53,] -0.311105295 -0.623291116
[54,] -0.054418355 -0.311105295
[55,] 0.231552901 -0.054418355
[56,] 0.054440294 0.231552901
[57,] 0.398469038 0.054440294
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.226172428 -0.405616120
2 0.230560930 0.226172428
3 -0.041912378 0.230560930
4 0.001756608 -0.041912378
5 0.051756608 0.001756608
6 0.084562356 0.051756608
7 0.130713492 0.084562356
8 -0.356222280 0.130713492
9 -0.249238945 -0.356222280
10 -0.151001578 -0.249238945
11 -0.044821916 -0.151001578
12 -0.716871975 -0.044821916
13 -0.581260477 -0.716871975
14 -0.369500355 -0.581260477
15 -0.484129762 -0.369500355
16 -0.457295269 -0.484129762
17 -0.457295269 -0.457295269
18 -0.261906425 -0.457295269
19 0.416510822 -0.261906425
20 -0.283819059 0.416510822
21 -0.227247313 -0.283819059
22 -0.192283014 -0.227247313
23 0.310917431 -0.192283014
24 0.182893752 0.310917431
25 0.524296626 0.182893752
26 0.555159863 0.524296626
27 0.088337392 0.555159863
28 0.137888322 0.088337392
29 0.447348684 0.137888322
30 0.756449288 0.447348684
31 1.100270983 0.756449288
32 0.329716360 1.100270983
33 0.419782368 0.329716360
34 0.599048537 0.419782368
35 0.608781923 0.599048537
36 -0.422092493 0.608781923
37 0.036450258 -0.422092493
38 0.187493374 0.036450258
39 -0.871747542 0.187493374
40 -0.776136044 -0.871747542
41 0.180884772 -0.776136044
42 0.215596874 0.180884772
43 0.398601766 0.215596874
44 0.229989001 0.398601766
45 0.524337277 0.229989001
46 0.371531529 0.524337277
47 0.430848205 0.371531529
48 -0.629724694 0.430848205
49 -0.326019516 -0.629724694
50 -0.498631309 -0.326019516
51 -0.669531564 -0.498631309
52 -0.623291116 -0.669531564
53 -0.311105295 -0.623291116
54 -0.054418355 -0.311105295
55 0.231552901 -0.054418355
56 0.054440294 0.231552901
57 0.398469038 0.054440294
> 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/79kk71321906094.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/8ymel1321906094.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/9z3w01321906094.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/101z1g1321906094.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/11kc0v1321906094.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/12nu6k1321906094.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/13gk0k1321906094.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/14nlpx1321906094.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/155a441321906094.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/16xjqx1321906094.tab")
+ }
>
> try(system("convert tmp/1mdpa1321906094.ps tmp/1mdpa1321906094.png",intern=TRUE))
character(0)
> try(system("convert tmp/20y311321906094.ps tmp/20y311321906094.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ekaj1321906094.ps tmp/3ekaj1321906094.png",intern=TRUE))
character(0)
> try(system("convert tmp/4bbtb1321906094.ps tmp/4bbtb1321906094.png",intern=TRUE))
character(0)
> try(system("convert tmp/5xmdf1321906094.ps tmp/5xmdf1321906094.png",intern=TRUE))
character(0)
> try(system("convert tmp/643l51321906094.ps tmp/643l51321906094.png",intern=TRUE))
character(0)
> try(system("convert tmp/79kk71321906094.ps tmp/79kk71321906094.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ymel1321906094.ps tmp/8ymel1321906094.png",intern=TRUE))
character(0)
> try(system("convert tmp/9z3w01321906094.ps tmp/9z3w01321906094.png",intern=TRUE))
character(0)
> try(system("convert tmp/101z1g1321906094.ps tmp/101z1g1321906094.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.578 0.559 4.177