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(-1
+ ,8.5
+ ,6
+ ,1.01
+ ,-2
+ ,8.6
+ ,6
+ ,1.02
+ ,-5
+ ,8.6
+ ,5
+ ,1.04
+ ,-4
+ ,8.6
+ ,5
+ ,1.06
+ ,-6
+ ,8.6
+ ,3
+ ,1.06
+ ,-2
+ ,8.4
+ ,5
+ ,1.06
+ ,-2
+ ,8
+ ,5
+ ,1.06
+ ,-2
+ ,7.9
+ ,5
+ ,1.06
+ ,-2
+ ,8
+ ,3
+ ,1.02
+ ,2
+ ,8
+ ,6
+ ,0.98
+ ,1
+ ,8
+ ,6
+ ,0.99
+ ,-8
+ ,8
+ ,4
+ ,0.99
+ ,-1
+ ,7.9
+ ,6
+ ,0.94
+ ,1
+ ,7.9
+ ,5
+ ,0.96
+ ,-1
+ ,7.9
+ ,4
+ ,0.98
+ ,2
+ ,8
+ ,5
+ ,1.01
+ ,2
+ ,7.9
+ ,5
+ ,1.01
+ ,1
+ ,7.5
+ ,4
+ ,1.02
+ ,-1
+ ,7.2
+ ,3
+ ,1.04
+ ,-2
+ ,7
+ ,2
+ ,1.03
+ ,-2
+ ,6.9
+ ,3
+ ,1.05
+ ,-1
+ ,7.1
+ ,2
+ ,1.08
+ ,-8
+ ,7.1
+ ,-1
+ ,1.17
+ ,-4
+ ,7.2
+ ,0
+ ,1.11
+ ,-6
+ ,7.1
+ ,-2
+ ,1.11
+ ,-3
+ ,6.9
+ ,1
+ ,1.11
+ ,-3
+ ,6.8
+ ,-2
+ ,1.2
+ ,-7
+ ,6.7
+ ,-2
+ ,1.21
+ ,-9
+ ,6.7
+ ,-2
+ ,1.31
+ ,-11
+ ,6.9
+ ,-6
+ ,1.37
+ ,-13
+ ,7.3
+ ,-4
+ ,1.37
+ ,-11
+ ,7.4
+ ,-2
+ ,1.26
+ ,-9
+ ,7.3
+ ,0
+ ,1.23
+ ,-17
+ ,7.1
+ ,-5
+ ,1.17
+ ,-22
+ ,7
+ ,-4
+ ,1.06
+ ,-25
+ ,7.1
+ ,-5
+ ,0.95
+ ,-20
+ ,7.5
+ ,-1
+ ,0.92
+ ,-24
+ ,7.7
+ ,-2
+ ,0.92
+ ,-24
+ ,7.8
+ ,-4
+ ,0.9
+ ,-22
+ ,7.7
+ ,-1
+ ,0.93
+ ,-19
+ ,7.7
+ ,1
+ ,0.93
+ ,-18
+ ,7.8
+ ,1
+ ,0.97
+ ,-17
+ ,8
+ ,-2
+ ,0.96
+ ,-11
+ ,8.1
+ ,1
+ ,0.99
+ ,-11
+ ,8.1
+ ,1
+ ,0.98
+ ,-12
+ ,8
+ ,3
+ ,0.96
+ ,-10
+ ,8.1
+ ,3
+ ,1
+ ,-15
+ ,8.2
+ ,1
+ ,0.99
+ ,-15
+ ,8.3
+ ,1
+ ,1.03
+ ,-15
+ ,8.4
+ ,0
+ ,1.02
+ ,-13
+ ,8.5
+ ,2
+ ,1.07
+ ,-8
+ ,8.5
+ ,2
+ ,1.13
+ ,-13
+ ,8.5
+ ,-1
+ ,1.15
+ ,-9
+ ,8.5
+ ,1
+ ,1.16
+ ,-7
+ ,8.5
+ ,0
+ ,1.14
+ ,-4
+ ,8.3
+ ,1
+ ,1.15
+ ,-4
+ ,8.2
+ ,1
+ ,1.15
+ ,-2
+ ,8.1
+ ,3
+ ,1.16
+ ,0
+ ,7.9
+ ,2
+ ,1.17
+ ,-2
+ ,7.6
+ ,0
+ ,1.22
+ ,-3
+ ,7.3
+ ,0
+ ,1.26
+ ,1
+ ,7.1
+ ,3
+ ,1.29
+ ,-2
+ ,7
+ ,-2
+ ,1.36
+ ,-1
+ ,7
+ ,0
+ ,1.38
+ ,1
+ ,7
+ ,1
+ ,1.37
+ ,-3
+ ,7
+ ,-1
+ ,1.37
+ ,-4
+ ,6.9
+ ,-2
+ ,1.37)
+ ,dim=c(4
+ ,67)
+ ,dimnames=list(c('Consumentenvertrouwen'
+ ,'Werkloosheidsgraad'
+ ,'Financiële_situatie_gezinnen'
+ ,'Diesel')
+ ,1:67))
> y <- array(NA,dim=c(4,67),dimnames=list(c('Consumentenvertrouwen','Werkloosheidsgraad','Financiële_situatie_gezinnen','Diesel'),1:67))
> 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 = '1'
> 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
Consumentenvertrouwen Werkloosheidsgraad Financi\353le_situatie_gezinnen
1 -1 8.5 6
2 -2 8.6 6
3 -5 8.6 5
4 -4 8.6 5
5 -6 8.6 3
6 -2 8.4 5
7 -2 8.0 5
8 -2 7.9 5
9 -2 8.0 3
10 2 8.0 6
11 1 8.0 6
12 -8 8.0 4
13 -1 7.9 6
14 1 7.9 5
15 -1 7.9 4
16 2 8.0 5
17 2 7.9 5
18 1 7.5 4
19 -1 7.2 3
20 -2 7.0 2
21 -2 6.9 3
22 -1 7.1 2
23 -8 7.1 -1
24 -4 7.2 0
25 -6 7.1 -2
26 -3 6.9 1
27 -3 6.8 -2
28 -7 6.7 -2
29 -9 6.7 -2
30 -11 6.9 -6
31 -13 7.3 -4
32 -11 7.4 -2
33 -9 7.3 0
34 -17 7.1 -5
35 -22 7.0 -4
36 -25 7.1 -5
37 -20 7.5 -1
38 -24 7.7 -2
39 -24 7.8 -4
40 -22 7.7 -1
41 -19 7.7 1
42 -18 7.8 1
43 -17 8.0 -2
44 -11 8.1 1
45 -11 8.1 1
46 -12 8.0 3
47 -10 8.1 3
48 -15 8.2 1
49 -15 8.3 1
50 -15 8.4 0
51 -13 8.5 2
52 -8 8.5 2
53 -13 8.5 -1
54 -9 8.5 1
55 -7 8.5 0
56 -4 8.3 1
57 -4 8.2 1
58 -2 8.1 3
59 0 7.9 2
60 -2 7.6 0
61 -3 7.3 0
62 1 7.1 3
63 -2 7.0 -2
64 -1 7.0 0
65 1 7.0 1
66 -3 7.0 -1
67 -4 6.9 -2
Diesel
1 1.01
2 1.02
3 1.04
4 1.06
5 1.06
6 1.06
7 1.06
8 1.06
9 1.02
10 0.98
11 0.99
12 0.99
13 0.94
14 0.96
15 0.98
16 1.01
17 1.01
18 1.02
19 1.04
20 1.03
21 1.05
22 1.08
23 1.17
24 1.11
25 1.11
26 1.11
27 1.20
28 1.21
29 1.31
30 1.37
31 1.37
32 1.26
33 1.23
34 1.17
35 1.06
36 0.95
37 0.92
38 0.92
39 0.90
40 0.93
41 0.93
42 0.97
43 0.96
44 0.99
45 0.98
46 0.96
47 1.00
48 0.99
49 1.03
50 1.02
51 1.07
52 1.13
53 1.15
54 1.16
55 1.14
56 1.15
57 1.15
58 1.16
59 1.17
60 1.22
61 1.26
62 1.29
63 1.36
64 1.38
65 1.37
66 1.37
67 1.37
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Werkloosheidsgraad
-15.041 -3.786
`Financi\353le_situatie_gezinnen` Diesel
2.479 31.124
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-6.2284 -2.4201 -0.0112 2.4424 6.3357
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -15.0409 8.3629 -1.799 0.0769 .
Werkloosheidsgraad -3.7864 0.8293 -4.566 2.36e-05 ***
`Financi\353le_situatie_gezinnen` 2.4794 0.1504 16.482 < 2e-16 ***
Diesel 31.1240 3.3049 9.418 1.23e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.088 on 63 degrees of freedom
Multiple R-squared: 0.8326, Adjusted R-squared: 0.8246
F-statistic: 104.4 on 3 and 63 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.0536462161 0.1072924322 0.946353784
[2,] 0.0177384202 0.0354768404 0.982261580
[3,] 0.0064652069 0.0129304139 0.993534793
[4,] 0.0017425589 0.0034851178 0.998257441
[5,] 0.0004671154 0.0009342307 0.999532885
[6,] 0.1319424153 0.2638848305 0.868057585
[7,] 0.0844321857 0.1688643714 0.915567814
[8,] 0.0790242793 0.1580485586 0.920975721
[9,] 0.0585546512 0.1171093024 0.941445349
[10,] 0.0731796427 0.1463592854 0.926820357
[11,] 0.0705409938 0.1410819876 0.929459006
[12,] 0.0504652491 0.1009304981 0.949534751
[13,] 0.0357923930 0.0715847859 0.964207607
[14,] 0.0270994467 0.0541988935 0.972900553
[15,] 0.0232815565 0.0465631129 0.976718444
[16,] 0.0202604154 0.0405208308 0.979739585
[17,] 0.0122749668 0.0245499337 0.987725033
[18,] 0.0133815691 0.0267631382 0.986618431
[19,] 0.0224443461 0.0448886922 0.977555654
[20,] 0.0209286409 0.0418572818 0.979071359
[21,] 0.0756225570 0.1512451139 0.924377443
[22,] 0.0761549390 0.1523098780 0.923845061
[23,] 0.0701900009 0.1403800019 0.929809999
[24,] 0.0559388709 0.1118777418 0.944061129
[25,] 0.0870427813 0.1740855626 0.912957219
[26,] 0.0921466219 0.1842932439 0.907853378
[27,] 0.1073149504 0.2146299007 0.892685050
[28,] 0.2125415692 0.4250831385 0.787458431
[29,] 0.6274095825 0.7451808351 0.372590418
[30,] 0.6303279490 0.7393441019 0.369672051
[31,] 0.6408841320 0.7182317360 0.359115868
[32,] 0.6485378241 0.7029243518 0.351462176
[33,] 0.6231620836 0.7536758329 0.376837916
[34,] 0.6505509224 0.6988981552 0.349449078
[35,] 0.7490611483 0.5018777034 0.250938852
[36,] 0.8513823084 0.2972353832 0.148617692
[37,] 0.8480827793 0.3038344415 0.151917221
[38,] 0.8261432771 0.3477134459 0.173856723
[39,] 0.8313240372 0.3373519255 0.168675963
[40,] 0.7987235461 0.4025529077 0.201276454
[41,] 0.7449561753 0.5100876494 0.255043825
[42,] 0.6796691357 0.6406617287 0.320330864
[43,] 0.6967803224 0.6064393552 0.303219678
[44,] 0.7218450820 0.5563098361 0.278154918
[45,] 0.9752481483 0.0495037034 0.024751852
[46,] 0.9707008786 0.0585982427 0.029299121
[47,] 0.9922933873 0.0154132254 0.007706613
[48,] 0.9979323143 0.0041353714 0.002067686
[49,] 0.9972727294 0.0054545413 0.002727271
[50,] 0.9933541580 0.0132916840 0.006645842
[51,] 0.9863213785 0.0273572431 0.013678622
[52,] 0.9941141050 0.0117717899 0.005885895
[53,] 0.9794496844 0.0411006312 0.020550316
[54,] 0.9373838057 0.1252323887 0.062616194
> postscript(file="/var/wessaorg/rcomp/tmp/15wwn1321976149.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/2x88u1321976149.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/3hyfx1321976149.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/41s4t1321976149.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/5tfyn1321976149.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 = 67
Frequency = 1
1 2 3 4 5 6
-0.08619139 -1.01878888 -2.16186091 -1.78434125 1.17447536 -0.54162662
7 8 9 10 11 12
-2.05619736 -2.43484005 4.14757994 1.95431570 0.64307553 -3.39810786
13 14 15 16 17 18
-0.17936630 3.67756166 3.53448963 3.50000349 3.12136081 2.77495821
19 20 21 22 23 24
1.49595812 2.52932123 -0.95121011 2.35176306 -0.01117356 3.75550185
25 26 27 28 29 30
6.33567578 1.14016548 5.39858619 0.70870333 -4.40369838 2.40377920
31 32 33 34 35 36
-3.04046668 -2.19699873 -4.60073752 0.90645968 -3.52794943 -0.24625656
37 38 39 40 41 42
-2.71559854 -3.47890486 2.48103478 -4.26955334 -6.22836996 -6.09468796
43 44 45 46 47 48
3.41206251 1.41875975 1.72999992 -3.98497903 -2.85129703 -2.20259756
49 50 51 52 53 54
-3.06891556 0.10037560 -4.03599918 -0.90344021 0.91230437 -0.35775241
55 56 57 58 59 60
4.74413624 4.19620239 3.81755970 0.16886023 3.57974300 3.84643070
61 62 63 64 65 66
0.46554197 -4.66368884 2.17602882 -2.40526814 -2.57343628 -1.61461966
67
-0.51385404
> postscript(file="/var/wessaorg/rcomp/tmp/6jg7p1321976149.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 = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.08619139 NA
1 -1.01878888 -0.08619139
2 -2.16186091 -1.01878888
3 -1.78434125 -2.16186091
4 1.17447536 -1.78434125
5 -0.54162662 1.17447536
6 -2.05619736 -0.54162662
7 -2.43484005 -2.05619736
8 4.14757994 -2.43484005
9 1.95431570 4.14757994
10 0.64307553 1.95431570
11 -3.39810786 0.64307553
12 -0.17936630 -3.39810786
13 3.67756166 -0.17936630
14 3.53448963 3.67756166
15 3.50000349 3.53448963
16 3.12136081 3.50000349
17 2.77495821 3.12136081
18 1.49595812 2.77495821
19 2.52932123 1.49595812
20 -0.95121011 2.52932123
21 2.35176306 -0.95121011
22 -0.01117356 2.35176306
23 3.75550185 -0.01117356
24 6.33567578 3.75550185
25 1.14016548 6.33567578
26 5.39858619 1.14016548
27 0.70870333 5.39858619
28 -4.40369838 0.70870333
29 2.40377920 -4.40369838
30 -3.04046668 2.40377920
31 -2.19699873 -3.04046668
32 -4.60073752 -2.19699873
33 0.90645968 -4.60073752
34 -3.52794943 0.90645968
35 -0.24625656 -3.52794943
36 -2.71559854 -0.24625656
37 -3.47890486 -2.71559854
38 2.48103478 -3.47890486
39 -4.26955334 2.48103478
40 -6.22836996 -4.26955334
41 -6.09468796 -6.22836996
42 3.41206251 -6.09468796
43 1.41875975 3.41206251
44 1.72999992 1.41875975
45 -3.98497903 1.72999992
46 -2.85129703 -3.98497903
47 -2.20259756 -2.85129703
48 -3.06891556 -2.20259756
49 0.10037560 -3.06891556
50 -4.03599918 0.10037560
51 -0.90344021 -4.03599918
52 0.91230437 -0.90344021
53 -0.35775241 0.91230437
54 4.74413624 -0.35775241
55 4.19620239 4.74413624
56 3.81755970 4.19620239
57 0.16886023 3.81755970
58 3.57974300 0.16886023
59 3.84643070 3.57974300
60 0.46554197 3.84643070
61 -4.66368884 0.46554197
62 2.17602882 -4.66368884
63 -2.40526814 2.17602882
64 -2.57343628 -2.40526814
65 -1.61461966 -2.57343628
66 -0.51385404 -1.61461966
67 NA -0.51385404
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.01878888 -0.08619139
[2,] -2.16186091 -1.01878888
[3,] -1.78434125 -2.16186091
[4,] 1.17447536 -1.78434125
[5,] -0.54162662 1.17447536
[6,] -2.05619736 -0.54162662
[7,] -2.43484005 -2.05619736
[8,] 4.14757994 -2.43484005
[9,] 1.95431570 4.14757994
[10,] 0.64307553 1.95431570
[11,] -3.39810786 0.64307553
[12,] -0.17936630 -3.39810786
[13,] 3.67756166 -0.17936630
[14,] 3.53448963 3.67756166
[15,] 3.50000349 3.53448963
[16,] 3.12136081 3.50000349
[17,] 2.77495821 3.12136081
[18,] 1.49595812 2.77495821
[19,] 2.52932123 1.49595812
[20,] -0.95121011 2.52932123
[21,] 2.35176306 -0.95121011
[22,] -0.01117356 2.35176306
[23,] 3.75550185 -0.01117356
[24,] 6.33567578 3.75550185
[25,] 1.14016548 6.33567578
[26,] 5.39858619 1.14016548
[27,] 0.70870333 5.39858619
[28,] -4.40369838 0.70870333
[29,] 2.40377920 -4.40369838
[30,] -3.04046668 2.40377920
[31,] -2.19699873 -3.04046668
[32,] -4.60073752 -2.19699873
[33,] 0.90645968 -4.60073752
[34,] -3.52794943 0.90645968
[35,] -0.24625656 -3.52794943
[36,] -2.71559854 -0.24625656
[37,] -3.47890486 -2.71559854
[38,] 2.48103478 -3.47890486
[39,] -4.26955334 2.48103478
[40,] -6.22836996 -4.26955334
[41,] -6.09468796 -6.22836996
[42,] 3.41206251 -6.09468796
[43,] 1.41875975 3.41206251
[44,] 1.72999992 1.41875975
[45,] -3.98497903 1.72999992
[46,] -2.85129703 -3.98497903
[47,] -2.20259756 -2.85129703
[48,] -3.06891556 -2.20259756
[49,] 0.10037560 -3.06891556
[50,] -4.03599918 0.10037560
[51,] -0.90344021 -4.03599918
[52,] 0.91230437 -0.90344021
[53,] -0.35775241 0.91230437
[54,] 4.74413624 -0.35775241
[55,] 4.19620239 4.74413624
[56,] 3.81755970 4.19620239
[57,] 0.16886023 3.81755970
[58,] 3.57974300 0.16886023
[59,] 3.84643070 3.57974300
[60,] 0.46554197 3.84643070
[61,] -4.66368884 0.46554197
[62,] 2.17602882 -4.66368884
[63,] -2.40526814 2.17602882
[64,] -2.57343628 -2.40526814
[65,] -1.61461966 -2.57343628
[66,] -0.51385404 -1.61461966
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.01878888 -0.08619139
2 -2.16186091 -1.01878888
3 -1.78434125 -2.16186091
4 1.17447536 -1.78434125
5 -0.54162662 1.17447536
6 -2.05619736 -0.54162662
7 -2.43484005 -2.05619736
8 4.14757994 -2.43484005
9 1.95431570 4.14757994
10 0.64307553 1.95431570
11 -3.39810786 0.64307553
12 -0.17936630 -3.39810786
13 3.67756166 -0.17936630
14 3.53448963 3.67756166
15 3.50000349 3.53448963
16 3.12136081 3.50000349
17 2.77495821 3.12136081
18 1.49595812 2.77495821
19 2.52932123 1.49595812
20 -0.95121011 2.52932123
21 2.35176306 -0.95121011
22 -0.01117356 2.35176306
23 3.75550185 -0.01117356
24 6.33567578 3.75550185
25 1.14016548 6.33567578
26 5.39858619 1.14016548
27 0.70870333 5.39858619
28 -4.40369838 0.70870333
29 2.40377920 -4.40369838
30 -3.04046668 2.40377920
31 -2.19699873 -3.04046668
32 -4.60073752 -2.19699873
33 0.90645968 -4.60073752
34 -3.52794943 0.90645968
35 -0.24625656 -3.52794943
36 -2.71559854 -0.24625656
37 -3.47890486 -2.71559854
38 2.48103478 -3.47890486
39 -4.26955334 2.48103478
40 -6.22836996 -4.26955334
41 -6.09468796 -6.22836996
42 3.41206251 -6.09468796
43 1.41875975 3.41206251
44 1.72999992 1.41875975
45 -3.98497903 1.72999992
46 -2.85129703 -3.98497903
47 -2.20259756 -2.85129703
48 -3.06891556 -2.20259756
49 0.10037560 -3.06891556
50 -4.03599918 0.10037560
51 -0.90344021 -4.03599918
52 0.91230437 -0.90344021
53 -0.35775241 0.91230437
54 4.74413624 -0.35775241
55 4.19620239 4.74413624
56 3.81755970 4.19620239
57 0.16886023 3.81755970
58 3.57974300 0.16886023
59 3.84643070 3.57974300
60 0.46554197 3.84643070
61 -4.66368884 0.46554197
62 2.17602882 -4.66368884
63 -2.40526814 2.17602882
64 -2.57343628 -2.40526814
65 -1.61461966 -2.57343628
66 -0.51385404 -1.61461966
> 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/7n2281321976149.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/8mtds1321976149.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/9qagn1321976149.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/10joi01321976149.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/11myxo1321976149.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/12q2ik1321976149.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/13ig711321976149.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/14xfrl1321976149.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/15iuu41321976149.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/16o7xe1321976149.tab")
+ }
>
> try(system("convert tmp/15wwn1321976149.ps tmp/15wwn1321976149.png",intern=TRUE))
character(0)
> try(system("convert tmp/2x88u1321976149.ps tmp/2x88u1321976149.png",intern=TRUE))
character(0)
> try(system("convert tmp/3hyfx1321976149.ps tmp/3hyfx1321976149.png",intern=TRUE))
character(0)
> try(system("convert tmp/41s4t1321976149.ps tmp/41s4t1321976149.png",intern=TRUE))
character(0)
> try(system("convert tmp/5tfyn1321976149.ps tmp/5tfyn1321976149.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jg7p1321976149.ps tmp/6jg7p1321976149.png",intern=TRUE))
character(0)
> try(system("convert tmp/7n2281321976149.ps tmp/7n2281321976149.png",intern=TRUE))
character(0)
> try(system("convert tmp/8mtds1321976149.ps tmp/8mtds1321976149.png",intern=TRUE))
character(0)
> try(system("convert tmp/9qagn1321976149.ps tmp/9qagn1321976149.png",intern=TRUE))
character(0)
> try(system("convert tmp/10joi01321976149.ps tmp/10joi01321976149.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.282 0.488 3.802