R version 2.12.1 (2010-12-16)
Copyright (C) 2010 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(101.82
+ ,107.34
+ ,93.63
+ ,99.85
+ ,101.76
+ ,101.68
+ ,107.34
+ ,93.63
+ ,99.91
+ ,102.37
+ ,101.68
+ ,107.34
+ ,93.63
+ ,99.87
+ ,102.38
+ ,102.45
+ ,107.34
+ ,96.13
+ ,99.86
+ ,102.86
+ ,102.45
+ ,107.34
+ ,96.13
+ ,100.10
+ ,102.87
+ ,102.45
+ ,107.34
+ ,96.13
+ ,100.10
+ ,102.92
+ ,102.45
+ ,107.34
+ ,96.13
+ ,100.12
+ ,102.95
+ ,102.45
+ ,107.34
+ ,96.13
+ ,99.95
+ ,103.02
+ ,102.45
+ ,112.60
+ ,96.13
+ ,99.94
+ ,104.08
+ ,102.52
+ ,112.60
+ ,96.13
+ ,100.18
+ ,104.16
+ ,102.52
+ ,112.60
+ ,96.13
+ ,100.31
+ ,104.24
+ ,102.85
+ ,112.60
+ ,96.13
+ ,100.65
+ ,104.33
+ ,102.85
+ ,112.61
+ ,96.13
+ ,100.65
+ ,104.73
+ ,102.85
+ ,112.61
+ ,96.13
+ ,100.69
+ ,104.86
+ ,103.25
+ ,112.61
+ ,96.13
+ ,101.26
+ ,105.03
+ ,103.25
+ ,112.61
+ ,98.73
+ ,101.26
+ ,105.62
+ ,103.25
+ ,112.61
+ ,98.73
+ ,101.38
+ ,105.63
+ ,103.25
+ ,112.61
+ ,98.73
+ ,101.38
+ ,105.63
+ ,104.45
+ ,112.61
+ ,98.73
+ ,101.38
+ ,105.94
+ ,104.45
+ ,112.61
+ ,98.73
+ ,101.44
+ ,106.61
+ ,104.45
+ ,118.65
+ ,98.73
+ ,101.40
+ ,107.69
+ ,104.80
+ ,118.65
+ ,98.73
+ ,101.40
+ ,107.78
+ ,104.80
+ ,118.65
+ ,98.73
+ ,100.58
+ ,107.93
+ ,105.29
+ ,118.65
+ ,98.73
+ ,100.58
+ ,108.48
+ ,105.29
+ ,114.29
+ ,98.73
+ ,100.58
+ ,108.14
+ ,105.29
+ ,114.29
+ ,98.73
+ ,100.59
+ ,108.48
+ ,105.29
+ ,114.29
+ ,98.73
+ ,100.81
+ ,108.48
+ ,106.04
+ ,114.29
+ ,101.67
+ ,100.75
+ ,108.89
+ ,105.94
+ ,114.29
+ ,101.67
+ ,100.75
+ ,108.93
+ ,105.94
+ ,114.29
+ ,101.67
+ ,100.96
+ ,109.21
+ ,105.94
+ ,114.29
+ ,101.67
+ ,101.31
+ ,109.47
+ ,106.28
+ ,114.29
+ ,101.67
+ ,101.64
+ ,109.80
+ ,106.48
+ ,123.33
+ ,101.67
+ ,101.46
+ ,111.73
+ ,107.19
+ ,123.33
+ ,101.67
+ ,101.73
+ ,111.85
+ ,108.14
+ ,123.33
+ ,101.67
+ ,101.73
+ ,112.12
+ ,108.22
+ ,123.33
+ ,101.67
+ ,101.64
+ ,112.15
+ ,108.22
+ ,123.33
+ ,101.67
+ ,101.77
+ ,112.17
+ ,108.61
+ ,123.33
+ ,101.67
+ ,101.74
+ ,112.67
+ ,108.61
+ ,123.33
+ ,101.67
+ ,101.89
+ ,112.80
+ ,108.61
+ ,123.33
+ ,107.94
+ ,101.89
+ ,113.44
+ ,108.61
+ ,123.33
+ ,107.94
+ ,101.93
+ ,113.53
+ ,109.06
+ ,123.33
+ ,107.94
+ ,101.93
+ ,114.53
+ ,109.06
+ ,123.33
+ ,107.94
+ ,102.32
+ ,114.51
+ ,112.93
+ ,123.33
+ ,107.94
+ ,102.41
+ ,115.05
+ ,115.84
+ ,129.03
+ ,107.94
+ ,103.58
+ ,116.67
+ ,118.57
+ ,128.76
+ ,107.94
+ ,104.12
+ ,117.07
+ ,118.57
+ ,128.76
+ ,107.94
+ ,104.10
+ ,116.92
+ ,118.86
+ ,128.76
+ ,107.94
+ ,104.15
+ ,117.00
+ ,118.98
+ ,128.76
+ ,107.94
+ ,104.15
+ ,117.02
+ ,119.27
+ ,128.76
+ ,107.94
+ ,104.16
+ ,117.35
+ ,119.39
+ ,128.76
+ ,107.94
+ ,102.94
+ ,117.36
+ ,119.49
+ ,128.76
+ ,110.30
+ ,103.07
+ ,117.82
+ ,119.59
+ ,128.76
+ ,110.30
+ ,103.04
+ ,117.88
+ ,120.12
+ ,128.76
+ ,110.30
+ ,103.06
+ ,118.24
+ ,120.14
+ ,128.76
+ ,110.30
+ ,103.05
+ ,118.50
+ ,120.14
+ ,128.76
+ ,110.30
+ ,102.95
+ ,118.80
+ ,120.14
+ ,132.63
+ ,110.30
+ ,102.95
+ ,119.76
+ ,120.14
+ ,132.63
+ ,110.30
+ ,103.05
+ ,120.09)
+ ,dim=c(5
+ ,58)
+ ,dimnames=list(c('Bioscoop'
+ ,'Schouwburgabonn'
+ ,'Daguitstap'
+ ,'HuurDVD'
+ ,'Vrijetijdsbesteding')
+ ,1:58))
> y <- array(NA,dim=c(5,58),dimnames=list(c('Bioscoop','Schouwburgabonn','Daguitstap','HuurDVD','Vrijetijdsbesteding'),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 = '5'
> 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
Vrijetijdsbesteding Bioscoop Schouwburgabonn Daguitstap HuurDVD
1 101.76 101.82 107.34 93.63 99.85
2 102.37 101.68 107.34 93.63 99.91
3 102.38 101.68 107.34 93.63 99.87
4 102.86 102.45 107.34 96.13 99.86
5 102.87 102.45 107.34 96.13 100.10
6 102.92 102.45 107.34 96.13 100.10
7 102.95 102.45 107.34 96.13 100.12
8 103.02 102.45 107.34 96.13 99.95
9 104.08 102.45 112.60 96.13 99.94
10 104.16 102.52 112.60 96.13 100.18
11 104.24 102.52 112.60 96.13 100.31
12 104.33 102.85 112.60 96.13 100.65
13 104.73 102.85 112.61 96.13 100.65
14 104.86 102.85 112.61 96.13 100.69
15 105.03 103.25 112.61 96.13 101.26
16 105.62 103.25 112.61 98.73 101.26
17 105.63 103.25 112.61 98.73 101.38
18 105.63 103.25 112.61 98.73 101.38
19 105.94 104.45 112.61 98.73 101.38
20 106.61 104.45 112.61 98.73 101.44
21 107.69 104.45 118.65 98.73 101.40
22 107.78 104.80 118.65 98.73 101.40
23 107.93 104.80 118.65 98.73 100.58
24 108.48 105.29 118.65 98.73 100.58
25 108.14 105.29 114.29 98.73 100.58
26 108.48 105.29 114.29 98.73 100.59
27 108.48 105.29 114.29 98.73 100.81
28 108.89 106.04 114.29 101.67 100.75
29 108.93 105.94 114.29 101.67 100.75
30 109.21 105.94 114.29 101.67 100.96
31 109.47 105.94 114.29 101.67 101.31
32 109.80 106.28 114.29 101.67 101.64
33 111.73 106.48 123.33 101.67 101.46
34 111.85 107.19 123.33 101.67 101.73
35 112.12 108.14 123.33 101.67 101.73
36 112.15 108.22 123.33 101.67 101.64
37 112.17 108.22 123.33 101.67 101.77
38 112.67 108.61 123.33 101.67 101.74
39 112.80 108.61 123.33 101.67 101.89
40 113.44 108.61 123.33 107.94 101.89
41 113.53 108.61 123.33 107.94 101.93
42 114.53 109.06 123.33 107.94 101.93
43 114.51 109.06 123.33 107.94 102.32
44 115.05 112.93 123.33 107.94 102.41
45 116.67 115.84 129.03 107.94 103.58
46 117.07 118.57 128.76 107.94 104.12
47 116.92 118.57 128.76 107.94 104.10
48 117.00 118.86 128.76 107.94 104.15
49 117.02 118.98 128.76 107.94 104.15
50 117.35 119.27 128.76 107.94 104.16
51 117.36 119.39 128.76 107.94 102.94
52 117.82 119.49 128.76 110.30 103.07
53 117.88 119.59 128.76 110.30 103.04
54 118.24 120.12 128.76 110.30 103.06
55 118.50 120.14 128.76 110.30 103.05
56 118.80 120.14 128.76 110.30 102.95
57 119.76 120.14 132.63 110.30 102.95
58 120.09 120.14 132.63 110.30 103.05
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Bioscoop Schouwburgabonn Daguitstap
29.7574 0.1174 0.3512 0.4417
HuurDVD
-0.1869
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.0357 -0.4095 -0.1323 0.2958 1.5117
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 29.75744 15.83978 1.879 0.06580 .
Bioscoop 0.11738 0.04311 2.723 0.00875 **
Schouwburgabonn 0.35123 0.03276 10.722 6.94e-15 ***
Daguitstap 0.44173 0.04837 9.133 1.82e-12 ***
HuurDVD -0.18690 0.19257 -0.971 0.33617
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6484 on 53 degrees of freedom
Multiple R-squared: 0.9874, Adjusted R-squared: 0.9865
F-statistic: 1040 on 4 and 53 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,] 1.840348e-03 3.680697e-03 0.9981596516
[2,] 1.876940e-04 3.753881e-04 0.9998123060
[3,] 3.725810e-04 7.451620e-04 0.9996274190
[4,] 8.972825e-05 1.794565e-04 0.9999102717
[5,] 6.695647e-04 1.339129e-03 0.9993304353
[6,] 1.408724e-03 2.817448e-03 0.9985912758
[7,] 1.240587e-03 2.481174e-03 0.9987594128
[8,] 4.178812e-04 8.357624e-04 0.9995821188
[9,] 2.159165e-04 4.318331e-04 0.9997840835
[10,] 1.004503e-04 2.009005e-04 0.9998995497
[11,] 5.138960e-05 1.027792e-04 0.9999486104
[12,] 3.297438e-05 6.594876e-05 0.9999670256
[13,] 2.094104e-04 4.188209e-04 0.9997905896
[14,] 2.442620e-04 4.885239e-04 0.9997557380
[15,] 4.199272e-04 8.398543e-04 0.9995800728
[16,] 5.405556e-03 1.081111e-02 0.9945944445
[17,] 3.643364e-02 7.286727e-02 0.9635663649
[18,] 1.701773e-01 3.403545e-01 0.8298227408
[19,] 3.466113e-01 6.932226e-01 0.6533887187
[20,] 3.997315e-01 7.994629e-01 0.6002685406
[21,] 3.943345e-01 7.886690e-01 0.6056655225
[22,] 4.105484e-01 8.210968e-01 0.5894516129
[23,] 3.965752e-01 7.931503e-01 0.6034248276
[24,] 3.761418e-01 7.522836e-01 0.6238581769
[25,] 3.967398e-01 7.934796e-01 0.6032602079
[26,] 3.960056e-01 7.920112e-01 0.6039943873
[27,] 3.805606e-01 7.611213e-01 0.6194393622
[28,] 5.430714e-01 9.138572e-01 0.4569286140
[29,] 6.072943e-01 7.854114e-01 0.3927056982
[30,] 6.331484e-01 7.337032e-01 0.3668516186
[31,] 5.621935e-01 8.756130e-01 0.4378064962
[32,] 4.994737e-01 9.989475e-01 0.5005262534
[33,] 7.528878e-01 4.942244e-01 0.2471121858
[34,] 9.430628e-01 1.138744e-01 0.0569371853
[35,] 9.075564e-01 1.848871e-01 0.0924435526
[36,] 8.585824e-01 2.828351e-01 0.1414175555
[37,] 9.903491e-01 1.930172e-02 0.0096508582
[38,] 9.990859e-01 1.828227e-03 0.0009141133
[39,] 9.989890e-01 2.022062e-03 0.0010110309
[40,] 9.982415e-01 3.517037e-03 0.0017585187
[41,] 9.940266e-01 1.194674e-02 0.0059733717
[42,] 9.786145e-01 4.277097e-02 0.0213854865
[43,] 9.372173e-01 1.255653e-01 0.0627826533
> postscript(file="/var/www/rcomp/tmp/174511321903298.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/rcomp/tmp/2fk441321903298.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/rcomp/tmp/3u42f1321903298.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/rcomp/tmp/4jxoq1321903298.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/rcomp/tmp/52l6i1321903298.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.34677940 0.29086810 0.29339208 -0.42317854 -0.36832242 -0.31832242
7 8 9 10 11 12
-0.28458441 -0.24635749 -1.03567650 -0.91903712 -0.81474006 -0.69992992
13 14 15 16 17 18
-0.30344218 -0.16596616 0.06361436 -0.49487591 -0.46244785 -0.46244785
19 20 21 22 23 24
-0.29330610 0.38790793 -0.66097456 -0.61205822 -0.61531661 -0.12283373
25 26 27 28 29 30
1.06851267 1.41038167 1.45149978 0.46357188 0.51531007 0.83455917
31 32 33 34 35 36
1.15997434 1.51174166 0.20953802 0.29666001 0.45514723 0.45893564
37 38 39 40 41 42
0.50323270 0.95184675 1.10988183 -1.01974663 -0.92227061 0.02490755
43 44 45 46 47 48
0.07779873 0.18035192 -0.32454533 -0.04924051 -0.20297852 -0.14767424
49 50 51 52 53 54
-0.14176007 0.15606819 -0.07603621 -0.64595312 -0.60329832 -0.30177270
55 56 57 58
-0.04598934 0.23532061 -0.16392493 0.18476512
> postscript(file="/var/www/rcomp/tmp/6kdqy1321903298.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.34677940 NA
1 0.29086810 -0.34677940
2 0.29339208 0.29086810
3 -0.42317854 0.29339208
4 -0.36832242 -0.42317854
5 -0.31832242 -0.36832242
6 -0.28458441 -0.31832242
7 -0.24635749 -0.28458441
8 -1.03567650 -0.24635749
9 -0.91903712 -1.03567650
10 -0.81474006 -0.91903712
11 -0.69992992 -0.81474006
12 -0.30344218 -0.69992992
13 -0.16596616 -0.30344218
14 0.06361436 -0.16596616
15 -0.49487591 0.06361436
16 -0.46244785 -0.49487591
17 -0.46244785 -0.46244785
18 -0.29330610 -0.46244785
19 0.38790793 -0.29330610
20 -0.66097456 0.38790793
21 -0.61205822 -0.66097456
22 -0.61531661 -0.61205822
23 -0.12283373 -0.61531661
24 1.06851267 -0.12283373
25 1.41038167 1.06851267
26 1.45149978 1.41038167
27 0.46357188 1.45149978
28 0.51531007 0.46357188
29 0.83455917 0.51531007
30 1.15997434 0.83455917
31 1.51174166 1.15997434
32 0.20953802 1.51174166
33 0.29666001 0.20953802
34 0.45514723 0.29666001
35 0.45893564 0.45514723
36 0.50323270 0.45893564
37 0.95184675 0.50323270
38 1.10988183 0.95184675
39 -1.01974663 1.10988183
40 -0.92227061 -1.01974663
41 0.02490755 -0.92227061
42 0.07779873 0.02490755
43 0.18035192 0.07779873
44 -0.32454533 0.18035192
45 -0.04924051 -0.32454533
46 -0.20297852 -0.04924051
47 -0.14767424 -0.20297852
48 -0.14176007 -0.14767424
49 0.15606819 -0.14176007
50 -0.07603621 0.15606819
51 -0.64595312 -0.07603621
52 -0.60329832 -0.64595312
53 -0.30177270 -0.60329832
54 -0.04598934 -0.30177270
55 0.23532061 -0.04598934
56 -0.16392493 0.23532061
57 0.18476512 -0.16392493
58 NA 0.18476512
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.29086810 -0.34677940
[2,] 0.29339208 0.29086810
[3,] -0.42317854 0.29339208
[4,] -0.36832242 -0.42317854
[5,] -0.31832242 -0.36832242
[6,] -0.28458441 -0.31832242
[7,] -0.24635749 -0.28458441
[8,] -1.03567650 -0.24635749
[9,] -0.91903712 -1.03567650
[10,] -0.81474006 -0.91903712
[11,] -0.69992992 -0.81474006
[12,] -0.30344218 -0.69992992
[13,] -0.16596616 -0.30344218
[14,] 0.06361436 -0.16596616
[15,] -0.49487591 0.06361436
[16,] -0.46244785 -0.49487591
[17,] -0.46244785 -0.46244785
[18,] -0.29330610 -0.46244785
[19,] 0.38790793 -0.29330610
[20,] -0.66097456 0.38790793
[21,] -0.61205822 -0.66097456
[22,] -0.61531661 -0.61205822
[23,] -0.12283373 -0.61531661
[24,] 1.06851267 -0.12283373
[25,] 1.41038167 1.06851267
[26,] 1.45149978 1.41038167
[27,] 0.46357188 1.45149978
[28,] 0.51531007 0.46357188
[29,] 0.83455917 0.51531007
[30,] 1.15997434 0.83455917
[31,] 1.51174166 1.15997434
[32,] 0.20953802 1.51174166
[33,] 0.29666001 0.20953802
[34,] 0.45514723 0.29666001
[35,] 0.45893564 0.45514723
[36,] 0.50323270 0.45893564
[37,] 0.95184675 0.50323270
[38,] 1.10988183 0.95184675
[39,] -1.01974663 1.10988183
[40,] -0.92227061 -1.01974663
[41,] 0.02490755 -0.92227061
[42,] 0.07779873 0.02490755
[43,] 0.18035192 0.07779873
[44,] -0.32454533 0.18035192
[45,] -0.04924051 -0.32454533
[46,] -0.20297852 -0.04924051
[47,] -0.14767424 -0.20297852
[48,] -0.14176007 -0.14767424
[49,] 0.15606819 -0.14176007
[50,] -0.07603621 0.15606819
[51,] -0.64595312 -0.07603621
[52,] -0.60329832 -0.64595312
[53,] -0.30177270 -0.60329832
[54,] -0.04598934 -0.30177270
[55,] 0.23532061 -0.04598934
[56,] -0.16392493 0.23532061
[57,] 0.18476512 -0.16392493
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.29086810 -0.34677940
2 0.29339208 0.29086810
3 -0.42317854 0.29339208
4 -0.36832242 -0.42317854
5 -0.31832242 -0.36832242
6 -0.28458441 -0.31832242
7 -0.24635749 -0.28458441
8 -1.03567650 -0.24635749
9 -0.91903712 -1.03567650
10 -0.81474006 -0.91903712
11 -0.69992992 -0.81474006
12 -0.30344218 -0.69992992
13 -0.16596616 -0.30344218
14 0.06361436 -0.16596616
15 -0.49487591 0.06361436
16 -0.46244785 -0.49487591
17 -0.46244785 -0.46244785
18 -0.29330610 -0.46244785
19 0.38790793 -0.29330610
20 -0.66097456 0.38790793
21 -0.61205822 -0.66097456
22 -0.61531661 -0.61205822
23 -0.12283373 -0.61531661
24 1.06851267 -0.12283373
25 1.41038167 1.06851267
26 1.45149978 1.41038167
27 0.46357188 1.45149978
28 0.51531007 0.46357188
29 0.83455917 0.51531007
30 1.15997434 0.83455917
31 1.51174166 1.15997434
32 0.20953802 1.51174166
33 0.29666001 0.20953802
34 0.45514723 0.29666001
35 0.45893564 0.45514723
36 0.50323270 0.45893564
37 0.95184675 0.50323270
38 1.10988183 0.95184675
39 -1.01974663 1.10988183
40 -0.92227061 -1.01974663
41 0.02490755 -0.92227061
42 0.07779873 0.02490755
43 0.18035192 0.07779873
44 -0.32454533 0.18035192
45 -0.04924051 -0.32454533
46 -0.20297852 -0.04924051
47 -0.14767424 -0.20297852
48 -0.14176007 -0.14767424
49 0.15606819 -0.14176007
50 -0.07603621 0.15606819
51 -0.64595312 -0.07603621
52 -0.60329832 -0.64595312
53 -0.30177270 -0.60329832
54 -0.04598934 -0.30177270
55 0.23532061 -0.04598934
56 -0.16392493 0.23532061
57 0.18476512 -0.16392493
> 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/rcomp/tmp/7w5hw1321903298.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/rcomp/tmp/84shb1321903298.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/rcomp/tmp/9xmdc1321903298.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/rcomp/tmp/106hh11321903298.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11ljg51321903298.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/rcomp/tmp/12h5vk1321903298.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/rcomp/tmp/13r13c1321903298.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/rcomp/tmp/14ftbc1321903298.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/rcomp/tmp/15zt831321903298.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/rcomp/tmp/16lpfi1321903298.tab")
+ }
>
> try(system("convert tmp/174511321903298.ps tmp/174511321903298.png",intern=TRUE))
character(0)
> try(system("convert tmp/2fk441321903298.ps tmp/2fk441321903298.png",intern=TRUE))
character(0)
> try(system("convert tmp/3u42f1321903298.ps tmp/3u42f1321903298.png",intern=TRUE))
character(0)
> try(system("convert tmp/4jxoq1321903298.ps tmp/4jxoq1321903298.png",intern=TRUE))
character(0)
> try(system("convert tmp/52l6i1321903298.ps tmp/52l6i1321903298.png",intern=TRUE))
character(0)
> try(system("convert tmp/6kdqy1321903298.ps tmp/6kdqy1321903298.png",intern=TRUE))
character(0)
> try(system("convert tmp/7w5hw1321903298.ps tmp/7w5hw1321903298.png",intern=TRUE))
character(0)
> try(system("convert tmp/84shb1321903298.ps tmp/84shb1321903298.png",intern=TRUE))
character(0)
> try(system("convert tmp/9xmdc1321903298.ps tmp/9xmdc1321903298.png",intern=TRUE))
character(0)
> try(system("convert tmp/106hh11321903298.ps tmp/106hh11321903298.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.336 0.684 7.450