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(-14
+ ,-20
+ ,36
+ ,-2
+ ,3
+ ,-7
+ ,-8
+ ,24
+ ,1
+ ,5
+ ,-9
+ ,-15
+ ,22
+ ,-1
+ ,4
+ ,-9
+ ,-13
+ ,17
+ ,-1
+ ,-4
+ ,-4
+ ,-6
+ ,8
+ ,-2
+ ,-1
+ ,-3
+ ,0
+ ,12
+ ,-1
+ ,3
+ ,1
+ ,5
+ ,5
+ ,1
+ ,2
+ ,-1
+ ,-1
+ ,6
+ ,0
+ ,2
+ ,-2
+ ,-5
+ ,5
+ ,-2
+ ,2
+ ,1
+ ,4
+ ,8
+ ,3
+ ,6
+ ,-3
+ ,-3
+ ,15
+ ,0
+ ,6
+ ,-2
+ ,3
+ ,16
+ ,0
+ ,6
+ ,0
+ ,8
+ ,17
+ ,2
+ ,6
+ ,-2
+ ,3
+ ,23
+ ,3
+ ,7
+ ,-4
+ ,3
+ ,24
+ ,1
+ ,4
+ ,-4
+ ,7
+ ,27
+ ,1
+ ,3
+ ,-7
+ ,4
+ ,31
+ ,0
+ ,0
+ ,-9
+ ,-4
+ ,40
+ ,1
+ ,6
+ ,-13
+ ,-6
+ ,47
+ ,-1
+ ,3
+ ,-8
+ ,8
+ ,43
+ ,2
+ ,1
+ ,-13
+ ,2
+ ,60
+ ,2
+ ,6
+ ,-15
+ ,-1
+ ,64
+ ,0
+ ,5
+ ,-15
+ ,-2
+ ,65
+ ,1
+ ,7
+ ,-15
+ ,0
+ ,65
+ ,1
+ ,4
+ ,-10
+ ,10
+ ,55
+ ,3
+ ,3
+ ,-12
+ ,3
+ ,57
+ ,3
+ ,6
+ ,-11
+ ,6
+ ,57
+ ,1
+ ,6
+ ,-11
+ ,7
+ ,57
+ ,1
+ ,5
+ ,-17
+ ,-4
+ ,65
+ ,-2
+ ,2
+ ,-18
+ ,-5
+ ,69
+ ,1
+ ,3
+ ,-19
+ ,-7
+ ,70
+ ,1
+ ,-2
+ ,-22
+ ,-10
+ ,71
+ ,-1
+ ,-4
+ ,-24
+ ,-21
+ ,71
+ ,-4
+ ,0
+ ,-24
+ ,-22
+ ,73
+ ,-2
+ ,1
+ ,-20
+ ,-16
+ ,68
+ ,-1
+ ,4
+ ,-25
+ ,-25
+ ,65
+ ,-5
+ ,-3
+ ,-22
+ ,-22
+ ,57
+ ,-4
+ ,-3
+ ,-17
+ ,-22
+ ,41
+ ,-5
+ ,0
+ ,-9
+ ,-19
+ ,21
+ ,0
+ ,6
+ ,-11
+ ,-21
+ ,21
+ ,-2
+ ,-1
+ ,-13
+ ,-31
+ ,17
+ ,-4
+ ,0
+ ,-11
+ ,-28
+ ,9
+ ,-6
+ ,-1
+ ,-9
+ ,-23
+ ,11
+ ,-2
+ ,1
+ ,-7
+ ,-17
+ ,6
+ ,-2
+ ,-4
+ ,-3
+ ,-12
+ ,-2
+ ,-2
+ ,-1
+ ,-3
+ ,-14
+ ,0
+ ,1
+ ,-1
+ ,-6
+ ,-18
+ ,5
+ ,-2
+ ,0
+ ,-4
+ ,-16
+ ,3
+ ,0
+ ,3
+ ,-8
+ ,-22
+ ,7
+ ,-1
+ ,0
+ ,-1
+ ,-9
+ ,4
+ ,2
+ ,8
+ ,-2
+ ,-10
+ ,8
+ ,3
+ ,8
+ ,-2
+ ,-10
+ ,9
+ ,2
+ ,8
+ ,-1
+ ,0
+ ,14
+ ,3
+ ,8
+ ,1
+ ,3
+ ,12
+ ,4
+ ,11
+ ,2
+ ,2
+ ,12
+ ,5
+ ,13
+ ,2
+ ,4
+ ,7
+ ,5
+ ,5
+ ,-1
+ ,-3
+ ,15
+ ,4
+ ,12
+ ,1
+ ,0
+ ,14
+ ,5
+ ,13
+ ,-1
+ ,-1
+ ,19
+ ,6
+ ,9
+ ,-8
+ ,-7
+ ,39
+ ,4
+ ,11
+ ,1
+ ,2
+ ,12
+ ,6
+ ,7
+ ,2
+ ,3
+ ,11
+ ,6
+ ,12
+ ,-2
+ ,-3
+ ,17
+ ,3
+ ,11
+ ,-2
+ ,-5
+ ,16
+ ,5
+ ,10
+ ,-2
+ ,0
+ ,25
+ ,5
+ ,13
+ ,-2
+ ,-3
+ ,24
+ ,5
+ ,14
+ ,-6
+ ,-7
+ ,28
+ ,3
+ ,10
+ ,-4
+ ,-7
+ ,25
+ ,5
+ ,13
+ ,-5
+ ,-7
+ ,31
+ ,5
+ ,12
+ ,-2
+ ,-4
+ ,24
+ ,6
+ ,13
+ ,-1
+ ,-3
+ ,24
+ ,6
+ ,17
+ ,-5
+ ,-6
+ ,33
+ ,5
+ ,15
+ ,-9
+ ,-10
+ ,37
+ ,4
+ ,6)
+ ,dim=c(5
+ ,73)
+ ,dimnames=list(c('CONSUMENTENVERTROUWEN'
+ ,'ALGEMENEECONOMISCHSITUATIE'
+ ,'WERKLOOSHEIDINBELGIË'
+ ,'FINANCIËLESITUATIEVANDEGEZINNEN'
+ ,'SPAARVERMOGENVANDEGEZINNEN')
+ ,1:73))
> y <- array(NA,dim=c(5,73),dimnames=list(c('CONSUMENTENVERTROUWEN','ALGEMENEECONOMISCHSITUATIE','WERKLOOSHEIDINBELGIË','FINANCIËLESITUATIEVANDEGEZINNEN','SPAARVERMOGENVANDEGEZINNEN'),1:73))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '2'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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
ALGEMENEECONOMISCHSITUATIE CONSUMENTENVERTROUWEN WERKLOOSHEIDINBELGI\303\213
1 -20 -14 36
2 -8 -7 24
3 -15 -9 22
4 -13 -9 17
5 -6 -4 8
6 0 -3 12
7 5 1 5
8 -1 -1 6
9 -5 -2 5
10 4 1 8
11 -3 -3 15
12 3 -2 16
13 8 0 17
14 3 -2 23
15 3 -4 24
16 7 -4 27
17 4 -7 31
18 -4 -9 40
19 -6 -13 47
20 8 -8 43
21 2 -13 60
22 -1 -15 64
23 -2 -15 65
24 0 -15 65
25 10 -10 55
26 3 -12 57
27 6 -11 57
28 7 -11 57
29 -4 -17 65
30 -5 -18 69
31 -7 -19 70
32 -10 -22 71
33 -21 -24 71
34 -22 -24 73
35 -16 -20 68
36 -25 -25 65
37 -22 -22 57
38 -22 -17 41
39 -19 -9 21
40 -21 -11 21
41 -31 -13 17
42 -28 -11 9
43 -23 -9 11
44 -17 -7 6
45 -12 -3 -2
46 -14 -3 0
47 -18 -6 5
48 -16 -4 3
49 -22 -8 7
50 -9 -1 4
51 -10 -2 8
52 -10 -2 9
53 0 -1 14
54 3 1 12
55 2 2 12
56 4 2 7
57 -3 -1 15
58 0 1 14
59 -1 -1 19
60 -7 -8 39
61 2 1 12
62 3 2 11
63 -3 -2 17
64 -5 -2 16
65 0 -2 25
66 -3 -2 24
67 -7 -6 28
68 -7 -4 25
69 -7 -5 31
70 -4 -2 24
71 -3 -1 24
72 -6 -5 33
73 -10 -9 37
FINANCI\303\213LESITUATIEVANDEGEZINNEN SPAARVERMOGENVANDEGEZINNEN t
1 -2 3 1
2 1 5 2
3 -1 4 3
4 -1 -4 4
5 -2 -1 5
6 -1 3 6
7 1 2 7
8 0 2 8
9 -2 2 9
10 3 6 10
11 0 6 11
12 0 6 12
13 2 6 13
14 3 7 14
15 1 4 15
16 1 3 16
17 0 0 17
18 1 6 18
19 -1 3 19
20 2 1 20
21 2 6 21
22 0 5 22
23 1 7 23
24 1 4 24
25 3 3 25
26 3 6 26
27 1 6 27
28 1 5 28
29 -2 2 29
30 1 3 30
31 1 -2 31
32 -1 -4 32
33 -4 0 33
34 -2 1 34
35 -1 4 35
36 -5 -3 36
37 -4 -3 37
38 -5 0 38
39 0 6 39
40 -2 -1 40
41 -4 0 41
42 -6 -1 42
43 -2 1 43
44 -2 -4 44
45 -2 -1 45
46 1 -1 46
47 -2 0 47
48 0 3 48
49 -1 0 49
50 2 8 50
51 3 8 51
52 2 8 52
53 3 8 53
54 4 11 54
55 5 13 55
56 5 5 56
57 4 12 57
58 5 13 58
59 6 9 59
60 4 11 60
61 6 7 61
62 6 12 62
63 3 11 63
64 5 10 64
65 5 13 65
66 5 14 66
67 3 10 67
68 5 13 68
69 5 12 69
70 6 13 70
71 6 17 71
72 5 15 72
73 4 6 73
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept)
-0.23850
CONSUMENTENVERTROUWEN
3.67659
`WERKLOOSHEIDINBELGI\303\213`
0.93127
`FINANCI\303\213LESITUATIEVANDEGEZINNEN`
-0.75361
SPAARVERMOGENVANDEGEZINNEN
-0.82139
t
-0.02641
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.44719 -0.78648 -0.03398 0.90511 2.15104
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.23850 0.38065 -0.627 0.5331
CONSUMENTENVERTROUWEN 3.67659 0.10502 35.009 < 2e-16
`WERKLOOSHEIDINBELGI\303\213` 0.93127 0.02635 35.337 < 2e-16
`FINANCI\303\213LESITUATIEVANDEGEZINNEN` -0.75361 0.14588 -5.166 2.33e-06
SPAARVERMOGENVANDEGEZINNEN -0.82139 0.05535 -14.840 < 2e-16
t -0.02641 0.01109 -2.382 0.0201
(Intercept)
CONSUMENTENVERTROUWEN ***
`WERKLOOSHEIDINBELGI\303\213` ***
`FINANCI\303\213LESITUATIEVANDEGEZINNEN` ***
SPAARVERMOGENVANDEGEZINNEN ***
t *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.177 on 67 degrees of freedom
Multiple R-squared: 0.9866, Adjusted R-squared: 0.9856
F-statistic: 985.6 on 5 and 67 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.8285979 0.3428042 0.17140210
[2,] 0.7326425 0.5347150 0.26735750
[3,] 0.6107226 0.7785548 0.38927742
[4,] 0.5710735 0.8578529 0.42892646
[5,] 0.5528078 0.8943843 0.44719217
[6,] 0.6945694 0.6108613 0.30543063
[7,] 0.6938616 0.6122768 0.30613839
[8,] 0.6593129 0.6813742 0.34068710
[9,] 0.6737110 0.6525779 0.32628897
[10,] 0.7270852 0.5458296 0.27291478
[11,] 0.6679569 0.6640861 0.33204306
[12,] 0.5853462 0.8293075 0.41465377
[13,] 0.5383537 0.9232926 0.46164632
[14,] 0.5104876 0.9790248 0.48951242
[15,] 0.4318969 0.8637938 0.56810312
[16,] 0.3833544 0.7667088 0.61664558
[17,] 0.3421605 0.6843211 0.65783945
[18,] 0.4709348 0.9418697 0.52906516
[19,] 0.3942483 0.7884966 0.60575168
[20,] 0.3236437 0.6472873 0.67635634
[21,] 0.2939306 0.5878613 0.70606937
[22,] 0.2727157 0.5454314 0.72728430
[23,] 0.4950286 0.9900572 0.50497142
[24,] 0.6395807 0.7208386 0.36041928
[25,] 0.5764520 0.8470960 0.42354798
[26,] 0.5373167 0.9253665 0.46268327
[27,] 0.6522755 0.6954490 0.34772451
[28,] 0.7213972 0.5572055 0.27860276
[29,] 0.8124649 0.3750702 0.18753508
[30,] 0.7625536 0.4748928 0.23744639
[31,] 0.7285583 0.5428834 0.27144170
[32,] 0.7075081 0.5849837 0.29249185
[33,] 0.6654429 0.6691143 0.33455714
[34,] 0.6209323 0.7581355 0.37906775
[35,] 0.5556639 0.8886723 0.44433615
[36,] 0.4860388 0.9720777 0.51396117
[37,] 0.4191093 0.8382186 0.58089070
[38,] 0.4641804 0.9283609 0.53581957
[39,] 0.4994931 0.9989862 0.50050692
[40,] 0.5058325 0.9883349 0.49416747
[41,] 0.4902886 0.9805771 0.50971145
[42,] 0.4319238 0.8638476 0.56807622
[43,] 0.3563083 0.7126166 0.64369168
[44,] 0.9174290 0.1651420 0.08257098
[45,] 0.8938882 0.2122236 0.10611181
[46,] 0.8952576 0.2094847 0.10474237
[47,] 0.8752701 0.2494599 0.12472993
[48,] 0.8217117 0.3565766 0.17828828
[49,] 0.7806249 0.4387502 0.21937511
[50,] 0.8495621 0.3008758 0.15043788
[51,] 0.8455170 0.3089659 0.15448296
[52,] 0.7709745 0.4580511 0.22902553
[53,] 0.7066767 0.5866467 0.29332333
[54,] 0.7721896 0.4556208 0.22781038
[55,] 0.6837328 0.6325344 0.31626720
[56,] 0.5382842 0.9234316 0.46171581
> postscript(file="/var/wessaorg/rcomp/tmp/1f6qn1322164125.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/2stmj1322164125.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/3qmuz1322164125.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/4guih1322164125.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/5hxw11322164125.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 = 73
Frequency = 1
1 2 3 4 5 6
-0.831600887 0.537555203 0.451055512 0.562660929 -0.701878529 1.962054664
7 8 9 10 11 12
-0.513167755 -0.818467221 -1.691430045 0.565053942 -0.481914007 0.936637934
13 14 15 16 17 18
-0.814169335 -2.447192814 0.029712500 0.440919931 1.554219388 -1.765643355
19 20 21 22 23 24
0.476821861 0.463424720 1.148163177 -0.525951801 -0.034410721 -0.472182998
25 26 27 28 29 30
1.169814632 2.151041172 -0.006366058 0.198649017 -0.890593984 1.169561757
31 32 33 34 35 36
-2.165679492 1.809207093 -0.786476840 -1.293984458 -2.099784794 1.339157075
37 38 39 40 41 42
1.539570950 -0.206084006 0.729441294 -1.147958881 -0.729132791 0.065632426
43 44 45 46 47 48
0.533571679 -0.243817067 -0.009423952 -1.584711053 -0.624334492 -0.117153090
49 50 51 52 53 54
1.672733963 0.588829176 0.320362856 -1.338111801 1.108977474 1.862544684
55 56 57 58 59 60
-0.391232616 -0.279627200 1.322532158 -0.497957198 -1.306679842 -0.033983251
61 62 63 64 65 66
-0.730940380 1.657118370 1.720025960 1.363538634 0.472705481 -0.748222878
67 68 69 70 71 72
1.466656306 0.905106923 -1.800905695 -1.710367359 -1.074971500 -0.120036845
73
-1.258511848
> postscript(file="/var/wessaorg/rcomp/tmp/6dly81322164125.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.831600887 NA
1 0.537555203 -0.831600887
2 0.451055512 0.537555203
3 0.562660929 0.451055512
4 -0.701878529 0.562660929
5 1.962054664 -0.701878529
6 -0.513167755 1.962054664
7 -0.818467221 -0.513167755
8 -1.691430045 -0.818467221
9 0.565053942 -1.691430045
10 -0.481914007 0.565053942
11 0.936637934 -0.481914007
12 -0.814169335 0.936637934
13 -2.447192814 -0.814169335
14 0.029712500 -2.447192814
15 0.440919931 0.029712500
16 1.554219388 0.440919931
17 -1.765643355 1.554219388
18 0.476821861 -1.765643355
19 0.463424720 0.476821861
20 1.148163177 0.463424720
21 -0.525951801 1.148163177
22 -0.034410721 -0.525951801
23 -0.472182998 -0.034410721
24 1.169814632 -0.472182998
25 2.151041172 1.169814632
26 -0.006366058 2.151041172
27 0.198649017 -0.006366058
28 -0.890593984 0.198649017
29 1.169561757 -0.890593984
30 -2.165679492 1.169561757
31 1.809207093 -2.165679492
32 -0.786476840 1.809207093
33 -1.293984458 -0.786476840
34 -2.099784794 -1.293984458
35 1.339157075 -2.099784794
36 1.539570950 1.339157075
37 -0.206084006 1.539570950
38 0.729441294 -0.206084006
39 -1.147958881 0.729441294
40 -0.729132791 -1.147958881
41 0.065632426 -0.729132791
42 0.533571679 0.065632426
43 -0.243817067 0.533571679
44 -0.009423952 -0.243817067
45 -1.584711053 -0.009423952
46 -0.624334492 -1.584711053
47 -0.117153090 -0.624334492
48 1.672733963 -0.117153090
49 0.588829176 1.672733963
50 0.320362856 0.588829176
51 -1.338111801 0.320362856
52 1.108977474 -1.338111801
53 1.862544684 1.108977474
54 -0.391232616 1.862544684
55 -0.279627200 -0.391232616
56 1.322532158 -0.279627200
57 -0.497957198 1.322532158
58 -1.306679842 -0.497957198
59 -0.033983251 -1.306679842
60 -0.730940380 -0.033983251
61 1.657118370 -0.730940380
62 1.720025960 1.657118370
63 1.363538634 1.720025960
64 0.472705481 1.363538634
65 -0.748222878 0.472705481
66 1.466656306 -0.748222878
67 0.905106923 1.466656306
68 -1.800905695 0.905106923
69 -1.710367359 -1.800905695
70 -1.074971500 -1.710367359
71 -0.120036845 -1.074971500
72 -1.258511848 -0.120036845
73 NA -1.258511848
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.537555203 -0.831600887
[2,] 0.451055512 0.537555203
[3,] 0.562660929 0.451055512
[4,] -0.701878529 0.562660929
[5,] 1.962054664 -0.701878529
[6,] -0.513167755 1.962054664
[7,] -0.818467221 -0.513167755
[8,] -1.691430045 -0.818467221
[9,] 0.565053942 -1.691430045
[10,] -0.481914007 0.565053942
[11,] 0.936637934 -0.481914007
[12,] -0.814169335 0.936637934
[13,] -2.447192814 -0.814169335
[14,] 0.029712500 -2.447192814
[15,] 0.440919931 0.029712500
[16,] 1.554219388 0.440919931
[17,] -1.765643355 1.554219388
[18,] 0.476821861 -1.765643355
[19,] 0.463424720 0.476821861
[20,] 1.148163177 0.463424720
[21,] -0.525951801 1.148163177
[22,] -0.034410721 -0.525951801
[23,] -0.472182998 -0.034410721
[24,] 1.169814632 -0.472182998
[25,] 2.151041172 1.169814632
[26,] -0.006366058 2.151041172
[27,] 0.198649017 -0.006366058
[28,] -0.890593984 0.198649017
[29,] 1.169561757 -0.890593984
[30,] -2.165679492 1.169561757
[31,] 1.809207093 -2.165679492
[32,] -0.786476840 1.809207093
[33,] -1.293984458 -0.786476840
[34,] -2.099784794 -1.293984458
[35,] 1.339157075 -2.099784794
[36,] 1.539570950 1.339157075
[37,] -0.206084006 1.539570950
[38,] 0.729441294 -0.206084006
[39,] -1.147958881 0.729441294
[40,] -0.729132791 -1.147958881
[41,] 0.065632426 -0.729132791
[42,] 0.533571679 0.065632426
[43,] -0.243817067 0.533571679
[44,] -0.009423952 -0.243817067
[45,] -1.584711053 -0.009423952
[46,] -0.624334492 -1.584711053
[47,] -0.117153090 -0.624334492
[48,] 1.672733963 -0.117153090
[49,] 0.588829176 1.672733963
[50,] 0.320362856 0.588829176
[51,] -1.338111801 0.320362856
[52,] 1.108977474 -1.338111801
[53,] 1.862544684 1.108977474
[54,] -0.391232616 1.862544684
[55,] -0.279627200 -0.391232616
[56,] 1.322532158 -0.279627200
[57,] -0.497957198 1.322532158
[58,] -1.306679842 -0.497957198
[59,] -0.033983251 -1.306679842
[60,] -0.730940380 -0.033983251
[61,] 1.657118370 -0.730940380
[62,] 1.720025960 1.657118370
[63,] 1.363538634 1.720025960
[64,] 0.472705481 1.363538634
[65,] -0.748222878 0.472705481
[66,] 1.466656306 -0.748222878
[67,] 0.905106923 1.466656306
[68,] -1.800905695 0.905106923
[69,] -1.710367359 -1.800905695
[70,] -1.074971500 -1.710367359
[71,] -0.120036845 -1.074971500
[72,] -1.258511848 -0.120036845
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.537555203 -0.831600887
2 0.451055512 0.537555203
3 0.562660929 0.451055512
4 -0.701878529 0.562660929
5 1.962054664 -0.701878529
6 -0.513167755 1.962054664
7 -0.818467221 -0.513167755
8 -1.691430045 -0.818467221
9 0.565053942 -1.691430045
10 -0.481914007 0.565053942
11 0.936637934 -0.481914007
12 -0.814169335 0.936637934
13 -2.447192814 -0.814169335
14 0.029712500 -2.447192814
15 0.440919931 0.029712500
16 1.554219388 0.440919931
17 -1.765643355 1.554219388
18 0.476821861 -1.765643355
19 0.463424720 0.476821861
20 1.148163177 0.463424720
21 -0.525951801 1.148163177
22 -0.034410721 -0.525951801
23 -0.472182998 -0.034410721
24 1.169814632 -0.472182998
25 2.151041172 1.169814632
26 -0.006366058 2.151041172
27 0.198649017 -0.006366058
28 -0.890593984 0.198649017
29 1.169561757 -0.890593984
30 -2.165679492 1.169561757
31 1.809207093 -2.165679492
32 -0.786476840 1.809207093
33 -1.293984458 -0.786476840
34 -2.099784794 -1.293984458
35 1.339157075 -2.099784794
36 1.539570950 1.339157075
37 -0.206084006 1.539570950
38 0.729441294 -0.206084006
39 -1.147958881 0.729441294
40 -0.729132791 -1.147958881
41 0.065632426 -0.729132791
42 0.533571679 0.065632426
43 -0.243817067 0.533571679
44 -0.009423952 -0.243817067
45 -1.584711053 -0.009423952
46 -0.624334492 -1.584711053
47 -0.117153090 -0.624334492
48 1.672733963 -0.117153090
49 0.588829176 1.672733963
50 0.320362856 0.588829176
51 -1.338111801 0.320362856
52 1.108977474 -1.338111801
53 1.862544684 1.108977474
54 -0.391232616 1.862544684
55 -0.279627200 -0.391232616
56 1.322532158 -0.279627200
57 -0.497957198 1.322532158
58 -1.306679842 -0.497957198
59 -0.033983251 -1.306679842
60 -0.730940380 -0.033983251
61 1.657118370 -0.730940380
62 1.720025960 1.657118370
63 1.363538634 1.720025960
64 0.472705481 1.363538634
65 -0.748222878 0.472705481
66 1.466656306 -0.748222878
67 0.905106923 1.466656306
68 -1.800905695 0.905106923
69 -1.710367359 -1.800905695
70 -1.074971500 -1.710367359
71 -0.120036845 -1.074971500
72 -1.258511848 -0.120036845
> 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/7zq8q1322164125.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/8o2x41322164125.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/99ns11322164125.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/10hek81322164125.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/111f7h1322164125.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/129pzr1322164125.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/13nvt01322164125.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/14gs9g1322164125.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/15ih771322164125.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/16k4ct1322164125.tab")
+ }
>
> try(system("convert tmp/1f6qn1322164125.ps tmp/1f6qn1322164125.png",intern=TRUE))
character(0)
> try(system("convert tmp/2stmj1322164125.ps tmp/2stmj1322164125.png",intern=TRUE))
character(0)
> try(system("convert tmp/3qmuz1322164125.ps tmp/3qmuz1322164125.png",intern=TRUE))
character(0)
> try(system("convert tmp/4guih1322164125.ps tmp/4guih1322164125.png",intern=TRUE))
character(0)
> try(system("convert tmp/5hxw11322164125.ps tmp/5hxw11322164125.png",intern=TRUE))
character(0)
> try(system("convert tmp/6dly81322164125.ps tmp/6dly81322164125.png",intern=TRUE))
character(0)
> try(system("convert tmp/7zq8q1322164125.ps tmp/7zq8q1322164125.png",intern=TRUE))
character(0)
> try(system("convert tmp/8o2x41322164125.ps tmp/8o2x41322164125.png",intern=TRUE))
character(0)
> try(system("convert tmp/99ns11322164125.ps tmp/99ns11322164125.png",intern=TRUE))
character(0)
> try(system("convert tmp/10hek81322164125.ps tmp/10hek81322164125.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.386 0.495 3.950