R version 2.12.0 (2010-10-15)
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(2011-11
+ ,-14
+ ,-20
+ ,36
+ ,-2
+ ,3
+ ,2011-10
+ ,-7
+ ,-8
+ ,24
+ ,1
+ ,5
+ ,2011-09
+ ,-9
+ ,-15
+ ,22
+ ,-1
+ ,4
+ ,2011-08
+ ,-9
+ ,-13
+ ,17
+ ,-1
+ ,-4
+ ,2011-07
+ ,-4
+ ,-6
+ ,8
+ ,-2
+ ,-1
+ ,2011-06
+ ,-3
+ ,0
+ ,12
+ ,-1
+ ,3
+ ,2011-05
+ ,1
+ ,5
+ ,5
+ ,1
+ ,2
+ ,2011-04
+ ,-1
+ ,-1
+ ,6
+ ,0
+ ,2
+ ,2011-03
+ ,-2
+ ,-5
+ ,5
+ ,-2
+ ,2
+ ,2011-02
+ ,1
+ ,4
+ ,8
+ ,3
+ ,6
+ ,2011-01
+ ,-3
+ ,-3
+ ,15
+ ,0
+ ,6
+ ,2010-12
+ ,-2
+ ,3
+ ,16
+ ,0
+ ,6
+ ,2010-11
+ ,0
+ ,8
+ ,17
+ ,2
+ ,6
+ ,2010-10
+ ,-2
+ ,3
+ ,23
+ ,3
+ ,7
+ ,2010-09
+ ,-4
+ ,3
+ ,24
+ ,1
+ ,4
+ ,2010-08
+ ,-4
+ ,7
+ ,27
+ ,1
+ ,3
+ ,2010-07
+ ,-7
+ ,4
+ ,31
+ ,0
+ ,0
+ ,2010-06
+ ,-9
+ ,-4
+ ,40
+ ,1
+ ,6
+ ,2010-05
+ ,-13
+ ,-6
+ ,47
+ ,-1
+ ,3
+ ,2010-04
+ ,-8
+ ,8
+ ,43
+ ,2
+ ,1
+ ,2010-03
+ ,-13
+ ,2
+ ,60
+ ,2
+ ,6
+ ,2010-02
+ ,-15
+ ,-1
+ ,64
+ ,0
+ ,5
+ ,2010-01
+ ,-15
+ ,-2
+ ,65
+ ,1
+ ,7
+ ,2009-12
+ ,-15
+ ,0
+ ,65
+ ,1
+ ,4
+ ,2009-11
+ ,-10
+ ,10
+ ,55
+ ,3
+ ,3
+ ,2009-10
+ ,-12
+ ,3
+ ,57
+ ,3
+ ,6
+ ,2009-09
+ ,-11
+ ,6
+ ,57
+ ,1
+ ,6
+ ,2009-08
+ ,-11
+ ,7
+ ,57
+ ,1
+ ,5
+ ,2009-07
+ ,-17
+ ,-4
+ ,65
+ ,-2
+ ,2
+ ,2009-06
+ ,-18
+ ,-5
+ ,69
+ ,1
+ ,3
+ ,2009-05
+ ,-19
+ ,-7
+ ,70
+ ,1
+ ,-2
+ ,2009-04
+ ,-22
+ ,-10
+ ,71
+ ,-1
+ ,-4
+ ,2009-03
+ ,-24
+ ,-21
+ ,71
+ ,-4
+ ,0
+ ,2009-02
+ ,-24
+ ,-22
+ ,73
+ ,-2
+ ,1
+ ,2009-01
+ ,-20
+ ,-16
+ ,68
+ ,-1
+ ,4
+ ,2008-12
+ ,-25
+ ,-25
+ ,65
+ ,-5
+ ,-3
+ ,2008-11
+ ,-22
+ ,-22
+ ,57
+ ,-4
+ ,-3
+ ,2008-10
+ ,-17
+ ,-22
+ ,41
+ ,-5
+ ,0
+ ,2008-09
+ ,-9
+ ,-19
+ ,21
+ ,0
+ ,6
+ ,2008-08
+ ,-11
+ ,-21
+ ,21
+ ,-2
+ ,-1
+ ,2008-07
+ ,-13
+ ,-31
+ ,17
+ ,-4
+ ,0
+ ,2008-06
+ ,-11
+ ,-28
+ ,9
+ ,-6
+ ,-1
+ ,2008-05
+ ,-9
+ ,-23
+ ,11
+ ,-2
+ ,1
+ ,2008-04
+ ,-7
+ ,-17
+ ,6
+ ,-2
+ ,-4
+ ,2008-03
+ ,-3
+ ,-12
+ ,-2
+ ,-2
+ ,-1
+ ,2008-02
+ ,-3
+ ,-14
+ ,0
+ ,1
+ ,-1
+ ,2008-01
+ ,-6
+ ,-18
+ ,5
+ ,-2
+ ,0
+ ,2007-12
+ ,-4
+ ,-16
+ ,3
+ ,0
+ ,3
+ ,2007-11
+ ,-8
+ ,-22
+ ,7
+ ,-1
+ ,0
+ ,2007-10
+ ,-1
+ ,-9
+ ,4
+ ,2
+ ,8
+ ,2007-09
+ ,-2
+ ,-10
+ ,8
+ ,3
+ ,8
+ ,2007-08
+ ,-2
+ ,-10
+ ,9
+ ,2
+ ,8
+ ,2007-07
+ ,-1
+ ,0
+ ,14
+ ,3
+ ,8
+ ,2007-06
+ ,1
+ ,3
+ ,12
+ ,4
+ ,11
+ ,2007-05
+ ,2
+ ,2
+ ,12
+ ,5
+ ,13
+ ,2007-04
+ ,2
+ ,4
+ ,7
+ ,5
+ ,5
+ ,2007-03
+ ,-1
+ ,-3
+ ,15
+ ,4
+ ,12
+ ,2007-02
+ ,1
+ ,0
+ ,14
+ ,5
+ ,13
+ ,2007-01
+ ,-1
+ ,-1
+ ,19
+ ,6
+ ,9
+ ,2006-12
+ ,-8
+ ,-7
+ ,39
+ ,4
+ ,11
+ ,2006-11
+ ,1
+ ,2
+ ,12
+ ,6
+ ,7
+ ,2006-10
+ ,2
+ ,3
+ ,11
+ ,6
+ ,12
+ ,2006-09
+ ,-2
+ ,-3
+ ,17
+ ,3
+ ,11
+ ,2006-08
+ ,-2
+ ,-5
+ ,16
+ ,5
+ ,10
+ ,2006-07
+ ,-2
+ ,0
+ ,25
+ ,5
+ ,13
+ ,2006-06
+ ,-2
+ ,-3
+ ,24
+ ,5
+ ,14
+ ,2006-05
+ ,-6
+ ,-7
+ ,28
+ ,3
+ ,10
+ ,2006-04
+ ,-4
+ ,-7
+ ,25
+ ,5
+ ,13
+ ,2006-03
+ ,-5
+ ,-7
+ ,31
+ ,5
+ ,12
+ ,2006-02
+ ,-2
+ ,-4
+ ,24
+ ,6
+ ,13
+ ,2006-01
+ ,-1
+ ,-3
+ ,24
+ ,6
+ ,17
+ ,2005-12
+ ,-5
+ ,-6
+ ,33
+ ,5
+ ,15
+ ,2005-11
+ ,-9
+ ,-10
+ ,37
+ ,4
+ ,6)
+ ,dim=c(6
+ ,73)
+ ,dimnames=list(c('JAARTAL'
+ ,'CONSUMENTENVERTROUWEN'
+ ,'ALGEMENEECONOMISCHSITUATIE'
+ ,'WERKLOOSHEIDINBELGIË'
+ ,'FINANCIËLESITUATIEVANDEGEZINNEN'
+ ,'SPAARVERMOGENVANDEGEZINNEN')
+ ,1:73))
> y <- array(NA,dim=c(6,73),dimnames=list(c('JAARTAL','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 = '3'
> #'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 JAARTAL CONSUMENTENVERTROUWEN
1 -20 2000 -14
2 -8 2001 -7
3 -15 2002 -9
4 -13 2003 -9
5 -6 2004 -4
6 0 2005 -3
7 5 2006 1
8 -1 2007 -1
9 -5 2008 -2
10 4 2009 1
11 -3 2010 -3
12 3 1998 -2
13 8 1999 0
14 3 2000 -2
15 3 2001 -4
16 7 2002 -4
17 4 2003 -7
18 -4 2004 -9
19 -6 2005 -13
20 8 2006 -8
21 2 2007 -13
22 -1 2008 -15
23 -2 2009 -15
24 0 1997 -15
25 10 1998 -10
26 3 1999 -12
27 6 2000 -11
28 7 2001 -11
29 -4 2002 -17
30 -5 2003 -18
31 -7 2004 -19
32 -10 2005 -22
33 -21 2006 -24
34 -22 2007 -24
35 -16 2008 -20
36 -25 1996 -25
37 -22 1997 -22
38 -22 1998 -17
39 -19 1999 -9
40 -21 2000 -11
41 -31 2001 -13
42 -28 2002 -11
43 -23 2003 -9
44 -17 2004 -7
45 -12 2005 -3
46 -14 2006 -3
47 -18 2007 -6
48 -16 1995 -4
49 -22 1996 -8
50 -9 1997 -1
51 -10 1998 -2
52 -10 1999 -2
53 0 2000 -1
54 3 2001 1
55 2 2002 2
56 4 2003 2
57 -3 2004 -1
58 0 2005 1
59 -1 2006 -1
60 -7 1994 -8
61 2 1995 1
62 3 1996 2
63 -3 1997 -2
64 -5 1998 -2
65 0 1999 -2
66 -3 2000 -2
67 -7 2001 -6
68 -7 2002 -4
69 -7 2003 -5
70 -4 2004 -2
71 -3 2005 -1
72 -6 1993 -5
73 -10 1994 -9
WERKLOOSHEIDINBELGI\303\213 FINANCI\303\213LESITUATIEVANDEGEZINNEN
1 36 -2
2 24 1
3 22 -1
4 17 -1
5 8 -2
6 12 -1
7 5 1
8 6 0
9 5 -2
10 8 3
11 15 0
12 16 0
13 17 2
14 23 3
15 24 1
16 27 1
17 31 0
18 40 1
19 47 -1
20 43 2
21 60 2
22 64 0
23 65 1
24 65 1
25 55 3
26 57 3
27 57 1
28 57 1
29 65 -2
30 69 1
31 70 1
32 71 -1
33 71 -4
34 73 -2
35 68 -1
36 65 -5
37 57 -4
38 41 -5
39 21 0
40 21 -2
41 17 -4
42 9 -6
43 11 -2
44 6 -2
45 -2 -2
46 0 1
47 5 -2
48 3 0
49 7 -1
50 4 2
51 8 3
52 9 2
53 14 3
54 12 4
55 12 5
56 7 5
57 15 4
58 14 5
59 19 6
60 39 4
61 12 6
62 11 6
63 17 3
64 16 5
65 25 5
66 24 5
67 28 3
68 25 5
69 31 5
70 24 6
71 24 6
72 33 5
73 37 4
SPAARVERMOGENVANDEGEZINNEN t
1 3 1
2 5 2
3 4 3
4 -4 4
5 -1 5
6 3 6
7 2 7
8 2 8
9 2 9
10 6 10
11 6 11
12 6 12
13 6 13
14 7 14
15 4 15
16 3 16
17 0 17
18 6 18
19 3 19
20 1 20
21 6 21
22 5 22
23 7 23
24 4 24
25 3 25
26 6 26
27 6 27
28 5 28
29 2 29
30 3 30
31 -2 31
32 -4 32
33 0 33
34 1 34
35 4 35
36 -3 36
37 -3 37
38 0 38
39 6 39
40 -1 40
41 0 41
42 -1 42
43 1 43
44 -4 44
45 -1 45
46 -1 46
47 0 47
48 3 48
49 0 49
50 8 50
51 8 51
52 8 52
53 8 53
54 11 54
55 13 55
56 5 56
57 12 57
58 13 58
59 9 59
60 11 60
61 7 61
62 12 62
63 11 63
64 10 64
65 13 65
66 14 66
67 10 67
68 13 68
69 12 69
70 13 70
71 17 71
72 15 72
73 6 73
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept)
177.98700
JAARTAL
-0.08891
CONSUMENTENVERTROUWEN
3.70130
`WERKLOOSHEIDINBELGI\303\213`
0.93737
`FINANCI\303\213LESITUATIEVANDEGEZINNEN`
-0.77659
SPAARVERMOGENVANDEGEZINNEN
-0.82740
t
-0.03131
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.76265 -0.72777 -0.01018 0.72235 2.07013
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 177.98700 71.27599 2.497 0.0150
JAARTAL -0.08891 0.03556 -2.501 0.0149
CONSUMENTENVERTROUWEN 3.70130 0.10161 36.427 < 2e-16
`WERKLOOSHEIDINBELGI\303\213` 0.93737 0.02550 36.767 < 2e-16
`FINANCI\303\213LESITUATIEVANDEGEZINNEN` -0.77659 0.14078 -5.516 6.22e-07
SPAARVERMOGENVANDEGEZINNEN -0.82740 0.05336 -15.507 < 2e-16
t -0.03131 0.01086 -2.884 0.0053
(Intercept) *
JAARTAL *
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.134 on 66 degrees of freedom
Multiple R-squared: 0.9877, Adjusted R-squared: 0.9866
F-statistic: 886.8 on 6 and 66 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.8648400 0.27031991 0.13515996
[2,] 0.7659428 0.46811442 0.23405721
[3,] 0.6450697 0.70986052 0.35493026
[4,] 0.6542766 0.69144674 0.34572337
[5,] 0.7979313 0.40413731 0.20206866
[6,] 0.8296204 0.34075920 0.17037960
[7,] 0.8163856 0.36722884 0.18361442
[8,] 0.8182325 0.36353509 0.18176755
[9,] 0.8705838 0.25883235 0.12941618
[10,] 0.8289187 0.34216269 0.17108135
[11,] 0.7695687 0.46086267 0.23043133
[12,] 0.7338641 0.53227186 0.26613593
[13,] 0.7120592 0.57588162 0.28794081
[14,] 0.6380356 0.72392884 0.36196442
[15,] 0.6094076 0.78118488 0.39059244
[16,] 0.5614305 0.87713895 0.43856947
[17,] 0.6649749 0.67005026 0.33502513
[18,] 0.5959749 0.80805010 0.40402505
[19,] 0.5197813 0.96043736 0.48021868
[20,] 0.5125945 0.97481110 0.48740555
[21,] 0.4791467 0.95829344 0.52085328
[22,] 0.7227081 0.55458386 0.27729193
[23,] 0.8613100 0.27737995 0.13868998
[24,] 0.8184359 0.36312827 0.18156413
[25,] 0.7793339 0.44133212 0.22066606
[26,] 0.7869257 0.42614861 0.21307430
[27,] 0.8076076 0.38478479 0.19239240
[28,] 0.8762266 0.24754682 0.12377341
[29,] 0.8343706 0.33125887 0.16562943
[30,] 0.8056811 0.38863780 0.19431890
[31,] 0.7887081 0.42258380 0.21129190
[32,] 0.7461306 0.50773873 0.25386936
[33,] 0.7130947 0.57381062 0.28690531
[34,] 0.6688237 0.66235250 0.33117625
[35,] 0.6026914 0.79461724 0.39730862
[36,] 0.5301610 0.93967795 0.46983897
[37,] 0.5074225 0.98515495 0.49257748
[38,] 0.4543366 0.90867321 0.54566340
[39,] 0.4792933 0.95858669 0.52070666
[40,] 0.4495230 0.89904595 0.55047703
[41,] 0.3965851 0.79317011 0.60341494
[42,] 0.3274896 0.65497918 0.67251041
[43,] 0.9739207 0.05215850 0.02607925
[44,] 0.9634965 0.07300694 0.03650347
[45,] 0.9673346 0.06533082 0.03266541
[46,] 0.9670970 0.06580595 0.03290297
[47,] 0.9453189 0.10936230 0.05468115
[48,] 0.9183569 0.16328613 0.08164307
[49,] 0.9215236 0.15695289 0.07847644
[50,] 0.8701411 0.25971773 0.12985886
[51,] 0.8103090 0.37938196 0.18969098
[52,] 0.7642059 0.47158826 0.23579413
[53,] 0.7225275 0.55494490 0.27747245
[54,] 0.6713002 0.65739955 0.32869977
> postscript(file="/var/www/rcomp/tmp/1auz11322163453.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/2gc7f1322163453.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/3jjtl1322163453.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/4qdxv1322163453.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/5q1v61322163453.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
-1.132367294 0.311712576 0.328712105 0.516591775 -0.727770981 2.027848947
7 8 9 10 11 12
-0.369787497 -0.560912343 -1.355187617 1.041534731 0.075630826 0.401340435
13 14 15 16 17 18
-1.265244748 -2.762647573 -0.212553662 0.268161168 1.484036074 -1.688481748
19 20 21 22 23 24
0.640003629 0.678141459 1.506596346 -0.100620721 0.513614216 -1.004198692
25 26 27 28 29 30
0.708969132 1.839256095 -0.294996673 -0.002173547 -0.985034803 1.244167732
31 32 33 34 35 36
-2.008669183 2.070127419 -0.427203880 -0.801151773 -1.540517704 0.844362150
37 38 39 40 41 42
1.136212568 -0.546564160 0.557932304 -1.264198992 -0.717663278 0.118336724
43 44 45 46 47 48
0.722349953 -0.010182459 0.285974646 -1.138786645 -0.103857056 -0.631976676
49 50 51 52 53 54
1.285202262 0.257349198 0.105982547 -1.487750361 1.020905495 1.872039328
55 56 57 58 59 60
-0.277660334 -0.089780664 1.650603112 -0.090429808 -0.787456500 -0.559707069
61 62 63 64 65 66
-1.198670994 1.294609494 1.438676362 1.222039446 0.388132320 -0.726877801
67 68 69 70 71 72
1.586317214 1.151405742 -1.478683696 -1.296804051 -0.568291655 -0.666399935
73
-1.713612958
> postscript(file="/var/www/rcomp/tmp/6ci1a1322163453.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 -1.132367294 NA
1 0.311712576 -1.132367294
2 0.328712105 0.311712576
3 0.516591775 0.328712105
4 -0.727770981 0.516591775
5 2.027848947 -0.727770981
6 -0.369787497 2.027848947
7 -0.560912343 -0.369787497
8 -1.355187617 -0.560912343
9 1.041534731 -1.355187617
10 0.075630826 1.041534731
11 0.401340435 0.075630826
12 -1.265244748 0.401340435
13 -2.762647573 -1.265244748
14 -0.212553662 -2.762647573
15 0.268161168 -0.212553662
16 1.484036074 0.268161168
17 -1.688481748 1.484036074
18 0.640003629 -1.688481748
19 0.678141459 0.640003629
20 1.506596346 0.678141459
21 -0.100620721 1.506596346
22 0.513614216 -0.100620721
23 -1.004198692 0.513614216
24 0.708969132 -1.004198692
25 1.839256095 0.708969132
26 -0.294996673 1.839256095
27 -0.002173547 -0.294996673
28 -0.985034803 -0.002173547
29 1.244167732 -0.985034803
30 -2.008669183 1.244167732
31 2.070127419 -2.008669183
32 -0.427203880 2.070127419
33 -0.801151773 -0.427203880
34 -1.540517704 -0.801151773
35 0.844362150 -1.540517704
36 1.136212568 0.844362150
37 -0.546564160 1.136212568
38 0.557932304 -0.546564160
39 -1.264198992 0.557932304
40 -0.717663278 -1.264198992
41 0.118336724 -0.717663278
42 0.722349953 0.118336724
43 -0.010182459 0.722349953
44 0.285974646 -0.010182459
45 -1.138786645 0.285974646
46 -0.103857056 -1.138786645
47 -0.631976676 -0.103857056
48 1.285202262 -0.631976676
49 0.257349198 1.285202262
50 0.105982547 0.257349198
51 -1.487750361 0.105982547
52 1.020905495 -1.487750361
53 1.872039328 1.020905495
54 -0.277660334 1.872039328
55 -0.089780664 -0.277660334
56 1.650603112 -0.089780664
57 -0.090429808 1.650603112
58 -0.787456500 -0.090429808
59 -0.559707069 -0.787456500
60 -1.198670994 -0.559707069
61 1.294609494 -1.198670994
62 1.438676362 1.294609494
63 1.222039446 1.438676362
64 0.388132320 1.222039446
65 -0.726877801 0.388132320
66 1.586317214 -0.726877801
67 1.151405742 1.586317214
68 -1.478683696 1.151405742
69 -1.296804051 -1.478683696
70 -0.568291655 -1.296804051
71 -0.666399935 -0.568291655
72 -1.713612958 -0.666399935
73 NA -1.713612958
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.311712576 -1.132367294
[2,] 0.328712105 0.311712576
[3,] 0.516591775 0.328712105
[4,] -0.727770981 0.516591775
[5,] 2.027848947 -0.727770981
[6,] -0.369787497 2.027848947
[7,] -0.560912343 -0.369787497
[8,] -1.355187617 -0.560912343
[9,] 1.041534731 -1.355187617
[10,] 0.075630826 1.041534731
[11,] 0.401340435 0.075630826
[12,] -1.265244748 0.401340435
[13,] -2.762647573 -1.265244748
[14,] -0.212553662 -2.762647573
[15,] 0.268161168 -0.212553662
[16,] 1.484036074 0.268161168
[17,] -1.688481748 1.484036074
[18,] 0.640003629 -1.688481748
[19,] 0.678141459 0.640003629
[20,] 1.506596346 0.678141459
[21,] -0.100620721 1.506596346
[22,] 0.513614216 -0.100620721
[23,] -1.004198692 0.513614216
[24,] 0.708969132 -1.004198692
[25,] 1.839256095 0.708969132
[26,] -0.294996673 1.839256095
[27,] -0.002173547 -0.294996673
[28,] -0.985034803 -0.002173547
[29,] 1.244167732 -0.985034803
[30,] -2.008669183 1.244167732
[31,] 2.070127419 -2.008669183
[32,] -0.427203880 2.070127419
[33,] -0.801151773 -0.427203880
[34,] -1.540517704 -0.801151773
[35,] 0.844362150 -1.540517704
[36,] 1.136212568 0.844362150
[37,] -0.546564160 1.136212568
[38,] 0.557932304 -0.546564160
[39,] -1.264198992 0.557932304
[40,] -0.717663278 -1.264198992
[41,] 0.118336724 -0.717663278
[42,] 0.722349953 0.118336724
[43,] -0.010182459 0.722349953
[44,] 0.285974646 -0.010182459
[45,] -1.138786645 0.285974646
[46,] -0.103857056 -1.138786645
[47,] -0.631976676 -0.103857056
[48,] 1.285202262 -0.631976676
[49,] 0.257349198 1.285202262
[50,] 0.105982547 0.257349198
[51,] -1.487750361 0.105982547
[52,] 1.020905495 -1.487750361
[53,] 1.872039328 1.020905495
[54,] -0.277660334 1.872039328
[55,] -0.089780664 -0.277660334
[56,] 1.650603112 -0.089780664
[57,] -0.090429808 1.650603112
[58,] -0.787456500 -0.090429808
[59,] -0.559707069 -0.787456500
[60,] -1.198670994 -0.559707069
[61,] 1.294609494 -1.198670994
[62,] 1.438676362 1.294609494
[63,] 1.222039446 1.438676362
[64,] 0.388132320 1.222039446
[65,] -0.726877801 0.388132320
[66,] 1.586317214 -0.726877801
[67,] 1.151405742 1.586317214
[68,] -1.478683696 1.151405742
[69,] -1.296804051 -1.478683696
[70,] -0.568291655 -1.296804051
[71,] -0.666399935 -0.568291655
[72,] -1.713612958 -0.666399935
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.311712576 -1.132367294
2 0.328712105 0.311712576
3 0.516591775 0.328712105
4 -0.727770981 0.516591775
5 2.027848947 -0.727770981
6 -0.369787497 2.027848947
7 -0.560912343 -0.369787497
8 -1.355187617 -0.560912343
9 1.041534731 -1.355187617
10 0.075630826 1.041534731
11 0.401340435 0.075630826
12 -1.265244748 0.401340435
13 -2.762647573 -1.265244748
14 -0.212553662 -2.762647573
15 0.268161168 -0.212553662
16 1.484036074 0.268161168
17 -1.688481748 1.484036074
18 0.640003629 -1.688481748
19 0.678141459 0.640003629
20 1.506596346 0.678141459
21 -0.100620721 1.506596346
22 0.513614216 -0.100620721
23 -1.004198692 0.513614216
24 0.708969132 -1.004198692
25 1.839256095 0.708969132
26 -0.294996673 1.839256095
27 -0.002173547 -0.294996673
28 -0.985034803 -0.002173547
29 1.244167732 -0.985034803
30 -2.008669183 1.244167732
31 2.070127419 -2.008669183
32 -0.427203880 2.070127419
33 -0.801151773 -0.427203880
34 -1.540517704 -0.801151773
35 0.844362150 -1.540517704
36 1.136212568 0.844362150
37 -0.546564160 1.136212568
38 0.557932304 -0.546564160
39 -1.264198992 0.557932304
40 -0.717663278 -1.264198992
41 0.118336724 -0.717663278
42 0.722349953 0.118336724
43 -0.010182459 0.722349953
44 0.285974646 -0.010182459
45 -1.138786645 0.285974646
46 -0.103857056 -1.138786645
47 -0.631976676 -0.103857056
48 1.285202262 -0.631976676
49 0.257349198 1.285202262
50 0.105982547 0.257349198
51 -1.487750361 0.105982547
52 1.020905495 -1.487750361
53 1.872039328 1.020905495
54 -0.277660334 1.872039328
55 -0.089780664 -0.277660334
56 1.650603112 -0.089780664
57 -0.090429808 1.650603112
58 -0.787456500 -0.090429808
59 -0.559707069 -0.787456500
60 -1.198670994 -0.559707069
61 1.294609494 -1.198670994
62 1.438676362 1.294609494
63 1.222039446 1.438676362
64 0.388132320 1.222039446
65 -0.726877801 0.388132320
66 1.586317214 -0.726877801
67 1.151405742 1.586317214
68 -1.478683696 1.151405742
69 -1.296804051 -1.478683696
70 -0.568291655 -1.296804051
71 -0.666399935 -0.568291655
72 -1.713612958 -0.666399935
> 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/70nmu1322163453.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/8yxn91322163453.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/9sna81322163453.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/100iok1322163453.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/111h5s1322163453.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/12v2nl1322163453.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/134isg1322163453.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/141tr11322163453.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/15hdwf1322163453.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/16bjp31322163454.tab")
+ }
>
> try(system("convert tmp/1auz11322163453.ps tmp/1auz11322163453.png",intern=TRUE))
character(0)
> try(system("convert tmp/2gc7f1322163453.ps tmp/2gc7f1322163453.png",intern=TRUE))
character(0)
> try(system("convert tmp/3jjtl1322163453.ps tmp/3jjtl1322163453.png",intern=TRUE))
character(0)
> try(system("convert tmp/4qdxv1322163453.ps tmp/4qdxv1322163453.png",intern=TRUE))
character(0)
> try(system("convert tmp/5q1v61322163453.ps tmp/5q1v61322163453.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ci1a1322163453.ps tmp/6ci1a1322163453.png",intern=TRUE))
character(0)
> try(system("convert tmp/70nmu1322163453.ps tmp/70nmu1322163453.png",intern=TRUE))
character(0)
> try(system("convert tmp/8yxn91322163453.ps tmp/8yxn91322163453.png",intern=TRUE))
character(0)
> try(system("convert tmp/9sna81322163453.ps tmp/9sna81322163453.png",intern=TRUE))
character(0)
> try(system("convert tmp/100iok1322163453.ps tmp/100iok1322163453.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.990 0.210 3.283