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(-6
+ ,-4
+ ,38
+ ,2
+ ,-3
+ ,-2
+ ,37
+ ,2.3
+ ,-2
+ ,2
+ ,32
+ ,2.8
+ ,-5
+ ,-5
+ ,32
+ ,2.4
+ ,-11
+ ,-15
+ ,44
+ ,2.3
+ ,-11
+ ,-16
+ ,43
+ ,2.7
+ ,-11
+ ,-18
+ ,42
+ ,2.7
+ ,-10
+ ,-13
+ ,38
+ ,2.9
+ ,-14
+ ,-23
+ ,37
+ ,3
+ ,-8
+ ,-10
+ ,35
+ ,2.2
+ ,-9
+ ,-10
+ ,37
+ ,2.3
+ ,-5
+ ,-6
+ ,33
+ ,2.8
+ ,-1
+ ,-3
+ ,24
+ ,2.8
+ ,-2
+ ,-4
+ ,24
+ ,2.8
+ ,-5
+ ,-7
+ ,31
+ ,2.2
+ ,-4
+ ,-7
+ ,25
+ ,2.6
+ ,-6
+ ,-7
+ ,28
+ ,2.8
+ ,-2
+ ,-3
+ ,24
+ ,2.5
+ ,-2
+ ,0
+ ,25
+ ,2.4
+ ,-2
+ ,-5
+ ,16
+ ,2.3
+ ,-2
+ ,-3
+ ,17
+ ,1.9
+ ,2
+ ,3
+ ,11
+ ,1.7
+ ,1
+ ,2
+ ,12
+ ,2
+ ,-8
+ ,-7
+ ,39
+ ,2.1
+ ,-1
+ ,-1
+ ,19
+ ,1.7
+ ,1
+ ,0
+ ,14
+ ,1.8
+ ,-1
+ ,-3
+ ,15
+ ,1.8
+ ,2
+ ,4
+ ,7
+ ,1.8
+ ,2
+ ,2
+ ,12
+ ,1.3
+ ,1
+ ,3
+ ,12
+ ,1.3
+ ,-1
+ ,0
+ ,14
+ ,1.3
+ ,-2
+ ,-10
+ ,9
+ ,1.2
+ ,-2
+ ,-10
+ ,8
+ ,1.4
+ ,-1
+ ,-9
+ ,4
+ ,2.2
+ ,-8
+ ,-22
+ ,7
+ ,2.9
+ ,-4
+ ,-16
+ ,3
+ ,3.1
+ ,-6
+ ,-18
+ ,5
+ ,3.5
+ ,-3
+ ,-14
+ ,0
+ ,3.6
+ ,-3
+ ,-12
+ ,-2
+ ,4.4
+ ,-7
+ ,-17
+ ,6
+ ,4.1
+ ,-9
+ ,-23
+ ,11
+ ,5.1
+ ,-11
+ ,-28
+ ,9
+ ,5.8
+ ,-13
+ ,-31
+ ,17
+ ,5.9
+ ,-11
+ ,-21
+ ,21
+ ,5.4
+ ,-9
+ ,-19
+ ,21
+ ,5.5
+ ,-17
+ ,-22
+ ,41
+ ,4.8
+ ,-22
+ ,-22
+ ,57
+ ,3.2
+ ,-25
+ ,-25
+ ,65
+ ,2.7
+ ,-20
+ ,-16
+ ,68
+ ,2.1
+ ,-24
+ ,-22
+ ,73
+ ,1.9
+ ,-24
+ ,-21
+ ,71
+ ,0.6
+ ,-22
+ ,-10
+ ,71
+ ,0.7
+ ,-19
+ ,-7
+ ,70
+ ,-0.2
+ ,-18
+ ,-5
+ ,69
+ ,-1
+ ,-17
+ ,-4
+ ,65
+ ,-1.7
+ ,-11
+ ,7
+ ,57
+ ,-0.7
+ ,-11
+ ,6
+ ,57
+ ,-1
+ ,-12
+ ,3
+ ,57
+ ,-0.9
+ ,-10
+ ,10
+ ,55
+ ,0
+ ,-15
+ ,0
+ ,65
+ ,0.3
+ ,-15
+ ,-2
+ ,65
+ ,0.8
+ ,-15
+ ,-1
+ ,64
+ ,0.8
+ ,-13
+ ,2
+ ,60
+ ,1.9
+ ,-8
+ ,8
+ ,43
+ ,2.1
+ ,-13
+ ,-6
+ ,47
+ ,2.5
+ ,-9
+ ,-4
+ ,40
+ ,2.7
+ ,-7
+ ,4
+ ,31
+ ,2.4
+ ,-4
+ ,7
+ ,27
+ ,2.4
+ ,-4
+ ,3
+ ,24
+ ,2.9
+ ,-2
+ ,3
+ ,23
+ ,3.1
+ ,0
+ ,8
+ ,17
+ ,3
+ ,-2
+ ,3
+ ,16
+ ,3.4
+ ,-3
+ ,-3
+ ,15
+ ,3.7
+ ,1
+ ,4
+ ,8
+ ,3.5
+ ,-2
+ ,-5
+ ,5
+ ,3.5
+ ,-1
+ ,-1
+ ,6
+ ,3.3
+ ,1
+ ,5
+ ,5
+ ,3.1
+ ,-3
+ ,0
+ ,12
+ ,3.4
+ ,-4
+ ,-6
+ ,8
+ ,4
+ ,-9
+ ,-13
+ ,17
+ ,3.4
+ ,-9
+ ,-15
+ ,22
+ ,3.4
+ ,-7
+ ,-8
+ ,24
+ ,3.4)
+ ,dim=c(4
+ ,82)
+ ,dimnames=list(c('Consumentenvertrouwen'
+ ,'Economische_situatie'
+ ,'Werkloosheid'
+ ,'HICP')
+ ,1:82))
> y <- array(NA,dim=c(4,82),dimnames=list(c('Consumentenvertrouwen','Economische_situatie','Werkloosheid','HICP'),1:82))
> 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 Economische_situatie Werkloosheid HICP
1 -6 -4 38 2.0
2 -3 -2 37 2.3
3 -2 2 32 2.8
4 -5 -5 32 2.4
5 -11 -15 44 2.3
6 -11 -16 43 2.7
7 -11 -18 42 2.7
8 -10 -13 38 2.9
9 -14 -23 37 3.0
10 -8 -10 35 2.2
11 -9 -10 37 2.3
12 -5 -6 33 2.8
13 -1 -3 24 2.8
14 -2 -4 24 2.8
15 -5 -7 31 2.2
16 -4 -7 25 2.6
17 -6 -7 28 2.8
18 -2 -3 24 2.5
19 -2 0 25 2.4
20 -2 -5 16 2.3
21 -2 -3 17 1.9
22 2 3 11 1.7
23 1 2 12 2.0
24 -8 -7 39 2.1
25 -1 -1 19 1.7
26 1 0 14 1.8
27 -1 -3 15 1.8
28 2 4 7 1.8
29 2 2 12 1.3
30 1 3 12 1.3
31 -1 0 14 1.3
32 -2 -10 9 1.2
33 -2 -10 8 1.4
34 -1 -9 4 2.2
35 -8 -22 7 2.9
36 -4 -16 3 3.1
37 -6 -18 5 3.5
38 -3 -14 0 3.6
39 -3 -12 -2 4.4
40 -7 -17 6 4.1
41 -9 -23 11 5.1
42 -11 -28 9 5.8
43 -13 -31 17 5.9
44 -11 -21 21 5.4
45 -9 -19 21 5.5
46 -17 -22 41 4.8
47 -22 -22 57 3.2
48 -25 -25 65 2.7
49 -20 -16 68 2.1
50 -24 -22 73 1.9
51 -24 -21 71 0.6
52 -22 -10 71 0.7
53 -19 -7 70 -0.2
54 -18 -5 69 -1.0
55 -17 -4 65 -1.7
56 -11 7 57 -0.7
57 -11 6 57 -1.0
58 -12 3 57 -0.9
59 -10 10 55 0.0
60 -15 0 65 0.3
61 -15 -2 65 0.8
62 -15 -1 64 0.8
63 -13 2 60 1.9
64 -8 8 43 2.1
65 -13 -6 47 2.5
66 -9 -4 40 2.7
67 -7 4 31 2.4
68 -4 7 27 2.4
69 -4 3 24 2.9
70 -2 3 23 3.1
71 0 8 17 3.0
72 -2 3 16 3.4
73 -3 -3 15 3.7
74 1 4 8 3.5
75 -2 -5 5 3.5
76 -1 -1 6 3.3
77 1 5 5 3.1
78 -3 0 12 3.4
79 -4 -6 8 4.0
80 -9 -13 17 3.4
81 -9 -15 22 3.4
82 -7 -8 24 3.4
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Economische_situatie Werkloosheid
3.5337 0.3295 -0.2721
HICP
-0.2379
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.8456 -1.3569 -0.4669 1.2923 4.7405
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.53366 0.77262 4.574 1.78e-05 ***
Economische_situatie 0.32954 0.02739 12.030 < 2e-16 ***
Werkloosheid -0.27211 0.01322 -20.590 < 2e-16 ***
HICP -0.23789 0.21604 -1.101 0.274
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.836 on 78 degrees of freedom
Multiple R-squared: 0.9299, Adjusted R-squared: 0.9272
F-statistic: 344.7 on 3 and 78 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.18982645 0.3796529040 0.8101735480
[2,] 0.15451851 0.3090370141 0.8454814930
[3,] 0.07543955 0.1508790917 0.9245604542
[4,] 0.03607973 0.0721594636 0.9639202682
[5,] 0.03077039 0.0615407807 0.9692296097
[6,] 0.02303185 0.0460637012 0.9769681494
[7,] 0.08752628 0.1750525655 0.9124737172
[8,] 0.07773638 0.1554727677 0.9222636161
[9,] 0.06893910 0.1378781951 0.9310609025
[10,] 0.05394949 0.1078989877 0.9460505062
[11,] 0.08181323 0.1636264515 0.9181867743
[12,] 0.08204861 0.1640972198 0.9179513901
[13,] 0.09463554 0.1892710820 0.9053644590
[14,] 0.07178624 0.1435724880 0.9282137560
[15,] 0.05419254 0.1083850857 0.9458074571
[16,] 0.03973350 0.0794669939 0.9602665030
[17,] 0.03082295 0.0616458966 0.9691770517
[18,] 0.04971143 0.0994228609 0.9502885695
[19,] 0.04986386 0.0997277297 0.9501361352
[20,] 0.07711098 0.1542219603 0.9228890198
[21,] 0.07638526 0.1527705297 0.9236147351
[22,] 0.09025794 0.1805158726 0.9097420637
[23,] 0.14787476 0.2957495279 0.8521252360
[24,] 0.14481995 0.2896399009 0.8551800495
[25,] 0.13973848 0.2794769595 0.8602615203
[26,] 0.26417233 0.5283446615 0.7358276693
[27,] 0.33831789 0.6766357790 0.6616821105
[28,] 0.39146962 0.7829392431 0.6085303785
[29,] 0.38594348 0.7718869653 0.6140565173
[30,] 0.37376027 0.7475205476 0.6262397262
[31,] 0.35119843 0.7023968657 0.6488015671
[32,] 0.32647241 0.6529448199 0.6735275901
[33,] 0.33442454 0.6688490878 0.6655754561
[34,] 0.37468022 0.7493604301 0.6253197849
[35,] 0.32497370 0.6499473928 0.6750263036
[36,] 0.27208014 0.5441602796 0.7279198602
[37,] 0.24404635 0.4880926969 0.7559536516
[38,] 0.22670554 0.4534110819 0.7732944591
[39,] 0.28022382 0.5604476420 0.7197761790
[40,] 0.55536168 0.8892766432 0.4446383216
[41,] 0.95050778 0.0989844414 0.0494922207
[42,] 0.99237172 0.0152565632 0.0076282816
[43,] 0.99396289 0.0120742295 0.0060371148
[44,] 0.99294782 0.0141043652 0.0070521826
[45,] 0.99139543 0.0172091308 0.0086045654
[46,] 0.99958364 0.0008327159 0.0004163580
[47,] 0.99952148 0.0009570426 0.0004785213
[48,] 0.99929449 0.0014110250 0.0007055125
[49,] 0.99891509 0.0021698150 0.0010849075
[50,] 0.99894676 0.0021064751 0.0010532375
[51,] 0.99860296 0.0027940811 0.0013970405
[52,] 0.99820757 0.0035848557 0.0017924278
[53,] 0.99828266 0.0034346783 0.0017173392
[54,] 0.99707901 0.0058419896 0.0029209948
[55,] 0.99654505 0.0069098927 0.0034549464
[56,] 0.99706601 0.0058679826 0.0029339913
[57,] 0.99509407 0.0098118530 0.0049059265
[58,] 0.99609135 0.0078172960 0.0039086480
[59,] 0.99412333 0.0117533323 0.0058766661
[60,] 0.99217185 0.0156562925 0.0078281462
[61,] 0.99647701 0.0070459897 0.0035229948
[62,] 0.99681876 0.0063624783 0.0031812392
[63,] 0.99612296 0.0077540718 0.0038770359
[64,] 0.99371304 0.0125739216 0.0062869608
[65,] 0.98479181 0.0304163851 0.0152081926
[66,] 0.97077343 0.0584531408 0.0292265704
[67,] 0.94484289 0.1103142151 0.0551571076
[68,] 0.91823873 0.1635225372 0.0817612686
[69,] 0.86709970 0.2658006091 0.1329003046
> postscript(file="/var/www/rcomp/tmp/12ldw1321983953.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/2r7xv1321983953.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/3w1fu1321983953.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/4cxzg1321983954.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/50rxn1321983954.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 = 82
Frequency = 1
1 2 3 4 5 6
2.600347875 4.740524388 3.180768060 2.392404689 2.929315119 3.081905929
7 8 9 10 11 12
3.468882690 1.780324509 0.827424464 1.808856070 1.376858626 3.089209097
13 14 15 16 17 18
3.651622565 2.981164344 2.731803534 2.194318568 1.058216881 2.580255693
19 20 21 22 23 24
1.839948195 1.014906943 0.532774353 0.875304965 0.548320416 1.884868970
25 26 27 28 29 30
1.370326478 1.704039659 0.964771797 -0.518875055 1.381797714 0.052255934
31 32 33 34 35 36
-0.414905128 0.496189716 0.271660831 0.044003514 -1.689110250 -0.707210211
37 38 39 40 41 42
-1.408757223 -1.063669381 -2.076654879 -2.323458459 -0.747784210 -1.477766208
43 44 45 46 47 48
-0.288497518 -0.614432905 0.750272493 -0.985488885 -2.012403416 -1.965868470
49 50 51 52 53 54
0.741842166 0.032048926 -1.150962899 -2.752133518 -1.226966274 -1.348468292
55 56 57 58 59 60
-1.932959970 -1.496884367 -1.238709460 -1.226295163 -1.863200603 -0.775347941
61 62 63 64 65 66
0.002680406 -0.598968173 -0.414342178 -1.969830529 -1.172662586 0.311084175
67 68 69 70 71 72
-2.845578128 -1.922630664 -1.301839156 0.473631959 -0.830506692 -1.359748763
73 74 75 76 77 78
-0.583238012 -0.842355979 -1.692800359 -1.786438594 -2.083373986 -2.459550621
79 80 81 82
-2.427993395 -2.814973488 -0.795355932 -0.557934792
> postscript(file="/var/www/rcomp/tmp/667ne1321983954.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 = 82
Frequency = 1
lag(myerror, k = 1) myerror
0 2.600347875 NA
1 4.740524388 2.600347875
2 3.180768060 4.740524388
3 2.392404689 3.180768060
4 2.929315119 2.392404689
5 3.081905929 2.929315119
6 3.468882690 3.081905929
7 1.780324509 3.468882690
8 0.827424464 1.780324509
9 1.808856070 0.827424464
10 1.376858626 1.808856070
11 3.089209097 1.376858626
12 3.651622565 3.089209097
13 2.981164344 3.651622565
14 2.731803534 2.981164344
15 2.194318568 2.731803534
16 1.058216881 2.194318568
17 2.580255693 1.058216881
18 1.839948195 2.580255693
19 1.014906943 1.839948195
20 0.532774353 1.014906943
21 0.875304965 0.532774353
22 0.548320416 0.875304965
23 1.884868970 0.548320416
24 1.370326478 1.884868970
25 1.704039659 1.370326478
26 0.964771797 1.704039659
27 -0.518875055 0.964771797
28 1.381797714 -0.518875055
29 0.052255934 1.381797714
30 -0.414905128 0.052255934
31 0.496189716 -0.414905128
32 0.271660831 0.496189716
33 0.044003514 0.271660831
34 -1.689110250 0.044003514
35 -0.707210211 -1.689110250
36 -1.408757223 -0.707210211
37 -1.063669381 -1.408757223
38 -2.076654879 -1.063669381
39 -2.323458459 -2.076654879
40 -0.747784210 -2.323458459
41 -1.477766208 -0.747784210
42 -0.288497518 -1.477766208
43 -0.614432905 -0.288497518
44 0.750272493 -0.614432905
45 -0.985488885 0.750272493
46 -2.012403416 -0.985488885
47 -1.965868470 -2.012403416
48 0.741842166 -1.965868470
49 0.032048926 0.741842166
50 -1.150962899 0.032048926
51 -2.752133518 -1.150962899
52 -1.226966274 -2.752133518
53 -1.348468292 -1.226966274
54 -1.932959970 -1.348468292
55 -1.496884367 -1.932959970
56 -1.238709460 -1.496884367
57 -1.226295163 -1.238709460
58 -1.863200603 -1.226295163
59 -0.775347941 -1.863200603
60 0.002680406 -0.775347941
61 -0.598968173 0.002680406
62 -0.414342178 -0.598968173
63 -1.969830529 -0.414342178
64 -1.172662586 -1.969830529
65 0.311084175 -1.172662586
66 -2.845578128 0.311084175
67 -1.922630664 -2.845578128
68 -1.301839156 -1.922630664
69 0.473631959 -1.301839156
70 -0.830506692 0.473631959
71 -1.359748763 -0.830506692
72 -0.583238012 -1.359748763
73 -0.842355979 -0.583238012
74 -1.692800359 -0.842355979
75 -1.786438594 -1.692800359
76 -2.083373986 -1.786438594
77 -2.459550621 -2.083373986
78 -2.427993395 -2.459550621
79 -2.814973488 -2.427993395
80 -0.795355932 -2.814973488
81 -0.557934792 -0.795355932
82 NA -0.557934792
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4.740524388 2.600347875
[2,] 3.180768060 4.740524388
[3,] 2.392404689 3.180768060
[4,] 2.929315119 2.392404689
[5,] 3.081905929 2.929315119
[6,] 3.468882690 3.081905929
[7,] 1.780324509 3.468882690
[8,] 0.827424464 1.780324509
[9,] 1.808856070 0.827424464
[10,] 1.376858626 1.808856070
[11,] 3.089209097 1.376858626
[12,] 3.651622565 3.089209097
[13,] 2.981164344 3.651622565
[14,] 2.731803534 2.981164344
[15,] 2.194318568 2.731803534
[16,] 1.058216881 2.194318568
[17,] 2.580255693 1.058216881
[18,] 1.839948195 2.580255693
[19,] 1.014906943 1.839948195
[20,] 0.532774353 1.014906943
[21,] 0.875304965 0.532774353
[22,] 0.548320416 0.875304965
[23,] 1.884868970 0.548320416
[24,] 1.370326478 1.884868970
[25,] 1.704039659 1.370326478
[26,] 0.964771797 1.704039659
[27,] -0.518875055 0.964771797
[28,] 1.381797714 -0.518875055
[29,] 0.052255934 1.381797714
[30,] -0.414905128 0.052255934
[31,] 0.496189716 -0.414905128
[32,] 0.271660831 0.496189716
[33,] 0.044003514 0.271660831
[34,] -1.689110250 0.044003514
[35,] -0.707210211 -1.689110250
[36,] -1.408757223 -0.707210211
[37,] -1.063669381 -1.408757223
[38,] -2.076654879 -1.063669381
[39,] -2.323458459 -2.076654879
[40,] -0.747784210 -2.323458459
[41,] -1.477766208 -0.747784210
[42,] -0.288497518 -1.477766208
[43,] -0.614432905 -0.288497518
[44,] 0.750272493 -0.614432905
[45,] -0.985488885 0.750272493
[46,] -2.012403416 -0.985488885
[47,] -1.965868470 -2.012403416
[48,] 0.741842166 -1.965868470
[49,] 0.032048926 0.741842166
[50,] -1.150962899 0.032048926
[51,] -2.752133518 -1.150962899
[52,] -1.226966274 -2.752133518
[53,] -1.348468292 -1.226966274
[54,] -1.932959970 -1.348468292
[55,] -1.496884367 -1.932959970
[56,] -1.238709460 -1.496884367
[57,] -1.226295163 -1.238709460
[58,] -1.863200603 -1.226295163
[59,] -0.775347941 -1.863200603
[60,] 0.002680406 -0.775347941
[61,] -0.598968173 0.002680406
[62,] -0.414342178 -0.598968173
[63,] -1.969830529 -0.414342178
[64,] -1.172662586 -1.969830529
[65,] 0.311084175 -1.172662586
[66,] -2.845578128 0.311084175
[67,] -1.922630664 -2.845578128
[68,] -1.301839156 -1.922630664
[69,] 0.473631959 -1.301839156
[70,] -0.830506692 0.473631959
[71,] -1.359748763 -0.830506692
[72,] -0.583238012 -1.359748763
[73,] -0.842355979 -0.583238012
[74,] -1.692800359 -0.842355979
[75,] -1.786438594 -1.692800359
[76,] -2.083373986 -1.786438594
[77,] -2.459550621 -2.083373986
[78,] -2.427993395 -2.459550621
[79,] -2.814973488 -2.427993395
[80,] -0.795355932 -2.814973488
[81,] -0.557934792 -0.795355932
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4.740524388 2.600347875
2 3.180768060 4.740524388
3 2.392404689 3.180768060
4 2.929315119 2.392404689
5 3.081905929 2.929315119
6 3.468882690 3.081905929
7 1.780324509 3.468882690
8 0.827424464 1.780324509
9 1.808856070 0.827424464
10 1.376858626 1.808856070
11 3.089209097 1.376858626
12 3.651622565 3.089209097
13 2.981164344 3.651622565
14 2.731803534 2.981164344
15 2.194318568 2.731803534
16 1.058216881 2.194318568
17 2.580255693 1.058216881
18 1.839948195 2.580255693
19 1.014906943 1.839948195
20 0.532774353 1.014906943
21 0.875304965 0.532774353
22 0.548320416 0.875304965
23 1.884868970 0.548320416
24 1.370326478 1.884868970
25 1.704039659 1.370326478
26 0.964771797 1.704039659
27 -0.518875055 0.964771797
28 1.381797714 -0.518875055
29 0.052255934 1.381797714
30 -0.414905128 0.052255934
31 0.496189716 -0.414905128
32 0.271660831 0.496189716
33 0.044003514 0.271660831
34 -1.689110250 0.044003514
35 -0.707210211 -1.689110250
36 -1.408757223 -0.707210211
37 -1.063669381 -1.408757223
38 -2.076654879 -1.063669381
39 -2.323458459 -2.076654879
40 -0.747784210 -2.323458459
41 -1.477766208 -0.747784210
42 -0.288497518 -1.477766208
43 -0.614432905 -0.288497518
44 0.750272493 -0.614432905
45 -0.985488885 0.750272493
46 -2.012403416 -0.985488885
47 -1.965868470 -2.012403416
48 0.741842166 -1.965868470
49 0.032048926 0.741842166
50 -1.150962899 0.032048926
51 -2.752133518 -1.150962899
52 -1.226966274 -2.752133518
53 -1.348468292 -1.226966274
54 -1.932959970 -1.348468292
55 -1.496884367 -1.932959970
56 -1.238709460 -1.496884367
57 -1.226295163 -1.238709460
58 -1.863200603 -1.226295163
59 -0.775347941 -1.863200603
60 0.002680406 -0.775347941
61 -0.598968173 0.002680406
62 -0.414342178 -0.598968173
63 -1.969830529 -0.414342178
64 -1.172662586 -1.969830529
65 0.311084175 -1.172662586
66 -2.845578128 0.311084175
67 -1.922630664 -2.845578128
68 -1.301839156 -1.922630664
69 0.473631959 -1.301839156
70 -0.830506692 0.473631959
71 -1.359748763 -0.830506692
72 -0.583238012 -1.359748763
73 -0.842355979 -0.583238012
74 -1.692800359 -0.842355979
75 -1.786438594 -1.692800359
76 -2.083373986 -1.786438594
77 -2.459550621 -2.083373986
78 -2.427993395 -2.459550621
79 -2.814973488 -2.427993395
80 -0.795355932 -2.814973488
81 -0.557934792 -0.795355932
> 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/7mhes1321983954.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/8xy201321983954.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/9jgwb1321983954.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/10nu6k1321983954.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/11vle71321983954.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/12jdoj1321983954.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/135l7b1321983954.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/14tvqg1321983954.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/15d58a1321983954.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/163ban1321983954.tab")
+ }
>
> try(system("convert tmp/12ldw1321983953.ps tmp/12ldw1321983953.png",intern=TRUE))
character(0)
> try(system("convert tmp/2r7xv1321983953.ps tmp/2r7xv1321983953.png",intern=TRUE))
character(0)
> try(system("convert tmp/3w1fu1321983953.ps tmp/3w1fu1321983953.png",intern=TRUE))
character(0)
> try(system("convert tmp/4cxzg1321983954.ps tmp/4cxzg1321983954.png",intern=TRUE))
character(0)
> try(system("convert tmp/50rxn1321983954.ps tmp/50rxn1321983954.png",intern=TRUE))
character(0)
> try(system("convert tmp/667ne1321983954.ps tmp/667ne1321983954.png",intern=TRUE))
character(0)
> try(system("convert tmp/7mhes1321983954.ps tmp/7mhes1321983954.png",intern=TRUE))
character(0)
> try(system("convert tmp/8xy201321983954.ps tmp/8xy201321983954.png",intern=TRUE))
character(0)
> try(system("convert tmp/9jgwb1321983954.ps tmp/9jgwb1321983954.png",intern=TRUE))
character(0)
> try(system("convert tmp/10nu6k1321983954.ps tmp/10nu6k1321983954.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.592 0.688 5.253