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
+ ,1
+ ,14
+ ,4
+ ,2
+ ,8
+ ,3
+ ,82
+ ,1
+ ,3
+ ,8
+ ,2
+ ,14
+ ,3
+ ,4
+ ,8
+ ,1
+ ,16
+ ,5
+ ,5
+ ,8
+ ,5
+ ,140
+ ,7
+ ,6
+ ,8
+ ,8
+ ,173
+ ,2
+ ,7
+ ,8
+ ,3
+ ,9
+ ,8
+ ,8
+ ,8
+ ,8
+ ,13
+ ,6
+ ,1
+ ,12
+ ,12
+ ,17
+ ,4
+ ,2
+ ,12
+ ,3
+ ,16
+ ,9
+ ,3
+ ,12
+ ,8
+ ,21
+ ,7
+ ,4
+ ,12
+ ,3
+ ,14
+ ,2
+ ,5
+ ,12
+ ,3
+ ,15
+ ,12
+ ,6
+ ,12
+ ,3
+ ,10
+ ,8
+ ,7
+ ,12
+ ,3
+ ,14
+ ,1
+ ,8
+ ,12
+ ,1
+ ,16
+ ,6
+ ,9
+ ,12
+ ,2
+ ,14
+ ,10
+ ,10
+ ,12
+ ,20
+ ,17
+ ,3
+ ,11
+ ,12
+ ,2
+ ,10
+ ,5
+ ,12
+ ,12
+ ,1
+ ,23
+ ,11
+ ,1
+ ,9
+ ,1
+ ,21
+ ,2
+ ,2
+ ,9
+ ,6
+ ,14
+ ,4
+ ,3
+ ,9
+ ,8
+ ,14
+ ,7
+ ,4
+ ,9
+ ,5
+ ,14
+ ,11
+ ,5
+ ,9
+ ,1
+ ,16
+ ,5
+ ,6
+ ,9
+ ,7
+ ,14
+ ,1
+ ,7
+ ,9
+ ,7
+ ,14
+ ,9
+ ,8
+ ,9
+ ,5
+ ,7
+ ,3
+ ,9
+ ,9
+ ,8
+ ,17
+ ,10
+ ,1
+ ,14
+ ,2
+ ,14
+ ,3
+ ,2
+ ,14
+ ,5
+ ,21
+ ,4
+ ,3
+ ,14
+ ,2
+ ,24
+ ,7
+ ,4
+ ,14
+ ,5
+ ,7
+ ,6
+ ,5
+ ,14
+ ,1
+ ,30
+ ,13
+ ,6
+ ,14
+ ,2
+ ,93
+ ,16
+ ,7
+ ,14
+ ,6
+ ,14
+ ,9
+ ,8
+ ,14
+ ,3
+ ,14
+ ,1
+ ,9
+ ,14
+ ,6
+ ,107
+ ,10
+ ,10
+ ,14
+ ,6
+ ,231
+ ,5
+ ,11
+ ,14
+ ,1
+ ,385
+ ,2
+ ,12
+ ,14
+ ,2
+ ,14
+ ,11
+ ,13
+ ,14
+ ,10
+ ,29
+ ,14
+ ,14
+ ,14
+ ,1
+ ,16
+ ,15
+ ,1
+ ,13
+ ,2
+ ,7
+ ,10
+ ,2
+ ,13
+ ,1
+ ,21
+ ,3
+ ,3
+ ,13
+ ,1
+ ,14
+ ,2
+ ,4
+ ,13
+ ,1
+ ,17
+ ,13
+ ,5
+ ,13
+ ,6
+ ,14
+ ,4
+ ,6
+ ,13
+ ,4
+ ,21
+ ,1
+ ,7
+ ,13
+ ,9
+ ,15
+ ,9
+ ,8
+ ,13
+ ,10
+ ,10
+ ,5
+ ,9
+ ,13
+ ,6
+ ,15
+ ,8
+ ,10
+ ,13
+ ,1
+ ,7
+ ,7
+ ,11
+ ,13
+ ,6
+ ,12
+ ,12
+ ,12
+ ,13
+ ,18
+ ,84
+ ,6
+ ,13
+ ,13
+ ,3
+ ,17
+ ,11
+ ,1
+ ,19
+ ,4
+ ,14
+ ,4
+ ,2
+ ,19
+ ,1
+ ,10
+ ,9
+ ,3
+ ,19
+ ,3
+ ,17
+ ,15
+ ,4
+ ,19
+ ,5
+ ,91
+ ,14
+ ,5
+ ,19
+ ,4
+ ,21
+ ,17
+ ,6
+ ,19
+ ,4
+ ,21
+ ,3
+ ,7
+ ,19
+ ,1
+ ,16
+ ,7
+ ,8
+ ,19
+ ,17
+ ,35
+ ,1
+ ,9
+ ,19
+ ,2
+ ,17
+ ,16
+ ,10
+ ,19
+ ,1
+ ,15
+ ,13
+ ,11
+ ,19
+ ,6
+ ,14
+ ,5
+ ,12
+ ,19
+ ,10
+ ,28
+ ,18
+ ,13
+ ,19
+ ,9
+ ,14
+ ,6
+ ,14
+ ,19
+ ,5
+ ,14
+ ,10
+ ,15
+ ,19
+ ,1
+ ,20
+ ,12
+ ,16
+ ,19
+ ,13
+ ,35
+ ,20
+ ,17
+ ,19
+ ,11
+ ,28
+ ,8
+ ,18
+ ,19
+ ,9
+ ,17
+ ,11
+ ,19
+ ,19
+ ,4
+ ,14
+ ,19
+ ,1
+ ,13
+ ,4
+ ,10
+ ,4
+ ,2
+ ,13
+ ,5
+ ,10
+ ,1
+ ,3
+ ,13
+ ,2
+ ,14
+ ,3
+ ,4
+ ,13
+ ,1
+ ,7
+ ,9
+ ,5
+ ,13
+ ,2
+ ,14
+ ,11
+ ,6
+ ,13
+ ,4
+ ,14
+ ,12
+ ,7
+ ,13
+ ,12
+ ,10
+ ,2
+ ,8
+ ,13
+ ,14
+ ,10
+ ,7
+ ,9
+ ,13
+ ,2
+ ,21
+ ,6
+ ,10
+ ,13
+ ,7
+ ,10
+ ,5
+ ,11
+ ,13
+ ,4
+ ,17
+ ,8
+ ,12
+ ,13
+ ,1
+ ,17
+ ,10
+ ,13
+ ,13
+ ,6
+ ,24
+ ,13
+ ,1
+ ,14
+ ,2
+ ,16
+ ,2
+ ,2
+ ,14
+ ,1
+ ,63
+ ,9
+ ,3
+ ,14
+ ,4
+ ,17
+ ,4
+ ,4
+ ,14
+ ,6
+ ,21
+ ,1
+ ,5
+ ,14
+ ,7
+ ,7
+ ,14
+ ,6
+ ,14
+ ,9
+ ,49
+ ,7
+ ,7
+ ,14
+ ,1
+ ,7
+ ,10
+ ,8
+ ,14
+ ,3
+ ,14
+ ,6
+ ,9
+ ,14
+ ,6
+ ,210
+ ,11
+ ,10
+ ,14
+ ,8
+ ,35
+ ,5
+ ,11
+ ,14
+ ,8
+ ,14
+ ,3
+ ,12
+ ,14
+ ,4
+ ,28
+ ,13
+ ,13
+ ,14
+ ,8
+ ,56
+ ,12
+ ,14
+ ,14
+ ,7
+ ,31
+ ,15)
+ ,dim=c(5
+ ,102)
+ ,dimnames=list(c('position'
+ ,'starters'
+ ,'last'
+ ,'since'
+ ,'number')
+ ,1:102))
> y <- array(NA,dim=c(5,102),dimnames=list(c('position','starters','last','since','number'),1:102))
> 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 = '3'
> 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
last position starters since number
1 1 1 8 14 4
2 3 2 8 82 1
3 2 3 8 14 3
4 1 4 8 16 5
5 5 5 8 140 7
6 8 6 8 173 2
7 3 7 8 9 8
8 8 8 8 13 6
9 12 1 12 17 4
10 3 2 12 16 9
11 8 3 12 21 7
12 3 4 12 14 2
13 3 5 12 15 12
14 3 6 12 10 8
15 3 7 12 14 1
16 1 8 12 16 6
17 2 9 12 14 10
18 20 10 12 17 3
19 2 11 12 10 5
20 1 12 12 23 11
21 1 1 9 21 2
22 6 2 9 14 4
23 8 3 9 14 7
24 5 4 9 14 11
25 1 5 9 16 5
26 7 6 9 14 1
27 7 7 9 14 9
28 5 8 9 7 3
29 8 9 9 17 10
30 2 1 14 14 3
31 5 2 14 21 4
32 2 3 14 24 7
33 5 4 14 7 6
34 1 5 14 30 13
35 2 6 14 93 16
36 6 7 14 14 9
37 3 8 14 14 1
38 6 9 14 107 10
39 6 10 14 231 5
40 1 11 14 385 2
41 2 12 14 14 11
42 10 13 14 29 14
43 1 14 14 16 15
44 2 1 13 7 10
45 1 2 13 21 3
46 1 3 13 14 2
47 1 4 13 17 13
48 6 5 13 14 4
49 4 6 13 21 1
50 9 7 13 15 9
51 10 8 13 10 5
52 6 9 13 15 8
53 1 10 13 7 7
54 6 11 13 12 12
55 18 12 13 84 6
56 3 13 13 17 11
57 4 1 19 14 4
58 1 2 19 10 9
59 3 3 19 17 15
60 5 4 19 91 14
61 4 5 19 21 17
62 4 6 19 21 3
63 1 7 19 16 7
64 17 8 19 35 1
65 2 9 19 17 16
66 1 10 19 15 13
67 6 11 19 14 5
68 10 12 19 28 18
69 9 13 19 14 6
70 5 14 19 14 10
71 1 15 19 20 12
72 13 16 19 35 20
73 11 17 19 28 8
74 9 18 19 17 11
75 4 19 19 14 19
76 4 1 13 10 4
77 5 2 13 10 1
78 2 3 13 14 3
79 1 4 13 7 9
80 2 5 13 14 11
81 4 6 13 14 12
82 12 7 13 10 2
83 14 8 13 10 7
84 2 9 13 21 6
85 7 10 13 10 5
86 4 11 13 17 8
87 1 12 13 17 10
88 6 13 13 24 13
89 2 1 14 16 2
90 1 2 14 63 9
91 4 3 14 17 4
92 6 4 14 21 1
93 7 5 14 7 14
94 9 6 14 49 7
95 1 7 14 7 10
96 3 8 14 14 6
97 6 9 14 210 11
98 8 10 14 35 5
99 8 11 14 14 3
100 4 12 14 28 13
101 8 13 14 56 12
102 7 14 14 31 15
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) position starters since number
3.426639 0.352334 0.044940 -0.001409 -0.204184
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-6.0871 -2.4022 -0.2642 1.4399 13.1472
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.426639 1.639668 2.090 0.039249 *
position 0.352334 0.099933 3.526 0.000647 ***
starters 0.044940 0.130371 0.345 0.731057
since -0.001409 0.007271 -0.194 0.846787
number -0.204184 0.093056 -2.194 0.030609 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.749 on 97 degrees of freedom
Multiple R-squared: 0.1293, Adjusted R-squared: 0.09341
F-statistic: 3.602 on 4 and 97 DF, p-value: 0.008809
> 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.1034299 0.20685990 0.89657005
[2,] 0.0403502 0.08070041 0.95964980
[3,] 0.1659999 0.33199973 0.83400013
[4,] 0.0973738 0.19474759 0.90262620
[5,] 0.4408827 0.88176533 0.55911734
[6,] 0.3946057 0.78921142 0.60539429
[7,] 0.3646714 0.72934270 0.63532865
[8,] 0.3498836 0.69976729 0.65011635
[9,] 0.3484179 0.69683582 0.65158209
[10,] 0.2832323 0.56646460 0.71676770
[11,] 0.9506448 0.09871034 0.04935517
[12,] 0.9635615 0.07287705 0.03643852
[13,] 0.9617796 0.07644088 0.03822044
[14,] 0.9545488 0.09090246 0.04545123
[15,] 0.9427868 0.11442645 0.05721322
[16,] 0.9536478 0.09270444 0.04635222
[17,] 0.9435458 0.11290840 0.05645420
[18,] 0.9351450 0.12970998 0.06485499
[19,] 0.9151608 0.16967837 0.08483919
[20,] 0.9080375 0.18392505 0.09196252
[21,] 0.8777207 0.24455868 0.12227934
[22,] 0.8755905 0.24881904 0.12440952
[23,] 0.8546823 0.29063547 0.14531774
[24,] 0.8161265 0.36774707 0.18387354
[25,] 0.7821535 0.43569306 0.21784653
[26,] 0.7339776 0.53204484 0.26602242
[27,] 0.6957722 0.60845552 0.30422776
[28,] 0.6404028 0.71919433 0.35959716
[29,] 0.5906738 0.81865239 0.40932619
[30,] 0.5769334 0.84613310 0.42306655
[31,] 0.5195943 0.96081132 0.48040566
[32,] 0.4608120 0.92162407 0.53918797
[33,] 0.6079723 0.78405542 0.39202771
[34,] 0.5921587 0.81568258 0.40784129
[35,] 0.6311814 0.73763717 0.36881858
[36,] 0.6464773 0.70704534 0.35352267
[37,] 0.5921493 0.81570135 0.40785068
[38,] 0.5705966 0.85880680 0.42940340
[39,] 0.5595342 0.88093159 0.44046579
[40,] 0.5101610 0.97967802 0.48983901
[41,] 0.4598055 0.91961097 0.54019452
[42,] 0.4160788 0.83215760 0.58392120
[43,] 0.4514339 0.90286788 0.54856606
[44,] 0.4809783 0.96195658 0.51902171
[45,] 0.4249406 0.84988118 0.57505941
[46,] 0.4566389 0.91327784 0.54336108
[47,] 0.4028348 0.80566964 0.59716518
[48,] 0.7745729 0.45085419 0.22542710
[49,] 0.7588587 0.48228254 0.24114127
[50,] 0.7115463 0.57690737 0.28845368
[51,] 0.6712350 0.65752998 0.32876499
[52,] 0.6200280 0.75994398 0.37997199
[53,] 0.5771492 0.84570159 0.42285080
[54,] 0.5309933 0.93801335 0.46900667
[55,] 0.4858662 0.97173249 0.51413375
[56,] 0.5089352 0.98212966 0.49106483
[57,] 0.8136034 0.37279326 0.18639663
[58,] 0.7818156 0.43636874 0.21818437
[59,] 0.7955560 0.40888790 0.20444395
[60,] 0.7532953 0.49340936 0.24670468
[61,] 0.7879791 0.42404183 0.21202091
[62,] 0.7451546 0.50969071 0.25484535
[63,] 0.7071180 0.58576395 0.29288198
[64,] 0.8218865 0.35622709 0.17811355
[65,] 0.9142991 0.17140182 0.08570091
[66,] 0.8949921 0.21001571 0.10500785
[67,] 0.8661112 0.26777754 0.13388877
[68,] 0.8368862 0.32622750 0.16311375
[69,] 0.7897860 0.42042806 0.21021403
[70,] 0.7340795 0.53184105 0.26592053
[71,] 0.6983766 0.60324671 0.30162335
[72,] 0.6611944 0.67761118 0.33880559
[73,] 0.6048689 0.79026216 0.39513108
[74,] 0.5276093 0.94478150 0.47239075
[75,] 0.6299163 0.74016738 0.37008369
[76,] 0.9664637 0.06707263 0.03353632
[77,] 0.9543870 0.09122595 0.04561298
[78,] 0.9489336 0.10213272 0.05106636
[79,] 0.9187510 0.16249806 0.08124903
[80,] 0.9279334 0.14413317 0.07206658
[81,] 0.8818202 0.23635954 0.11817977
[82,] 0.8326309 0.33473822 0.16736911
[83,] 0.8075127 0.38497458 0.19248729
[84,] 0.7195960 0.56080802 0.28040401
[85,] 0.5931105 0.81377905 0.40688952
[86,] 0.6316284 0.73674326 0.36837163
[87,] 0.9611701 0.07765978 0.03882989
> postscript(file="/var/wessaorg/rcomp/tmp/1vv5m1322131164.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/222ru1322131164.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/3w76a1322131164.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/4mo9a1322131164.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/5037o1322131164.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 = 102
Frequency = 1
1 2 3 4 5 6
-2.30203960 -1.17113701 -2.21089184 -3.15204128 1.07866442 2.75189679
7 8 9 10 11 12
-1.60635328 2.63857956 8.52242492 0.18960081 4.43594230 -1.94717127
13 14 15 16 17 18
-0.25625930 -1.43237181 -3.20835776 -4.53695592 -3.07537244 13.14723296
19 20 21 22 23 24
-4.80659431 -4.91551357 -2.74548693 2.30068580 4.56090283 2.02530361
25 26 27 28 29 30
-3.54931588 1.27879755 2.55993337 -1.02736396 3.06367458 -1.77586550
31 32 33 34 35 36
1.08584456 -1.64971247 0.76982250 -2.12082652 -0.77186463 1.33523158
37 38 39 40 41 42
-3.65057272 0.96575117 -0.23282943 -5.98078307 -4.01807213 4.26327463
43 44 45 46 47 48
-4.90318829 -0.31149938 -3.07339884 -3.63977739 -1.74186435 1.06392164
49 50 51 52 53 54
-1.89110334 4.38158059 4.20546807 0.47272834 -5.09505885 0.58056894
55 56 57 58 59 60
11.10455485 -3.32124007 0.20361648 -2.13343358 0.74919526 2.29691726
61 62 63 64 65 66
1.45852888 -1.75237796 -4.29502042 10.15430712 -2.16062644 -4.12832926
67 68 69 70 71 72
-1.11554221 5.20623347 1.38397306 -2.15162616 -6.08714099 7.21512455
73 74 75 76 77 78
2.40272468 0.64744657 -3.07564356 0.46762402 0.50273851 -2.43559363
79 80 81 82 83 84
-2.57268587 -1.50679206 0.34505745 5.94525104 8.61383558 -3.92718728
85 86 87 88 89 90
0.50079958 -2.22912285 -5.17308958 0.09698799 -1.97723196 -1.83407342
91 92 93 94 95 96
-0.27212428 0.76862480 4.05095832 4.32850101 -3.47044520 -2.62965393
97 98 99 100 101 102
1.31502573 1.49107544 0.70079206 -1.58998353 1.89294063 1.11794144
> postscript(file="/var/wessaorg/rcomp/tmp/6xf111322131164.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 = 102
Frequency = 1
lag(myerror, k = 1) myerror
0 -2.30203960 NA
1 -1.17113701 -2.30203960
2 -2.21089184 -1.17113701
3 -3.15204128 -2.21089184
4 1.07866442 -3.15204128
5 2.75189679 1.07866442
6 -1.60635328 2.75189679
7 2.63857956 -1.60635328
8 8.52242492 2.63857956
9 0.18960081 8.52242492
10 4.43594230 0.18960081
11 -1.94717127 4.43594230
12 -0.25625930 -1.94717127
13 -1.43237181 -0.25625930
14 -3.20835776 -1.43237181
15 -4.53695592 -3.20835776
16 -3.07537244 -4.53695592
17 13.14723296 -3.07537244
18 -4.80659431 13.14723296
19 -4.91551357 -4.80659431
20 -2.74548693 -4.91551357
21 2.30068580 -2.74548693
22 4.56090283 2.30068580
23 2.02530361 4.56090283
24 -3.54931588 2.02530361
25 1.27879755 -3.54931588
26 2.55993337 1.27879755
27 -1.02736396 2.55993337
28 3.06367458 -1.02736396
29 -1.77586550 3.06367458
30 1.08584456 -1.77586550
31 -1.64971247 1.08584456
32 0.76982250 -1.64971247
33 -2.12082652 0.76982250
34 -0.77186463 -2.12082652
35 1.33523158 -0.77186463
36 -3.65057272 1.33523158
37 0.96575117 -3.65057272
38 -0.23282943 0.96575117
39 -5.98078307 -0.23282943
40 -4.01807213 -5.98078307
41 4.26327463 -4.01807213
42 -4.90318829 4.26327463
43 -0.31149938 -4.90318829
44 -3.07339884 -0.31149938
45 -3.63977739 -3.07339884
46 -1.74186435 -3.63977739
47 1.06392164 -1.74186435
48 -1.89110334 1.06392164
49 4.38158059 -1.89110334
50 4.20546807 4.38158059
51 0.47272834 4.20546807
52 -5.09505885 0.47272834
53 0.58056894 -5.09505885
54 11.10455485 0.58056894
55 -3.32124007 11.10455485
56 0.20361648 -3.32124007
57 -2.13343358 0.20361648
58 0.74919526 -2.13343358
59 2.29691726 0.74919526
60 1.45852888 2.29691726
61 -1.75237796 1.45852888
62 -4.29502042 -1.75237796
63 10.15430712 -4.29502042
64 -2.16062644 10.15430712
65 -4.12832926 -2.16062644
66 -1.11554221 -4.12832926
67 5.20623347 -1.11554221
68 1.38397306 5.20623347
69 -2.15162616 1.38397306
70 -6.08714099 -2.15162616
71 7.21512455 -6.08714099
72 2.40272468 7.21512455
73 0.64744657 2.40272468
74 -3.07564356 0.64744657
75 0.46762402 -3.07564356
76 0.50273851 0.46762402
77 -2.43559363 0.50273851
78 -2.57268587 -2.43559363
79 -1.50679206 -2.57268587
80 0.34505745 -1.50679206
81 5.94525104 0.34505745
82 8.61383558 5.94525104
83 -3.92718728 8.61383558
84 0.50079958 -3.92718728
85 -2.22912285 0.50079958
86 -5.17308958 -2.22912285
87 0.09698799 -5.17308958
88 -1.97723196 0.09698799
89 -1.83407342 -1.97723196
90 -0.27212428 -1.83407342
91 0.76862480 -0.27212428
92 4.05095832 0.76862480
93 4.32850101 4.05095832
94 -3.47044520 4.32850101
95 -2.62965393 -3.47044520
96 1.31502573 -2.62965393
97 1.49107544 1.31502573
98 0.70079206 1.49107544
99 -1.58998353 0.70079206
100 1.89294063 -1.58998353
101 1.11794144 1.89294063
102 NA 1.11794144
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.17113701 -2.30203960
[2,] -2.21089184 -1.17113701
[3,] -3.15204128 -2.21089184
[4,] 1.07866442 -3.15204128
[5,] 2.75189679 1.07866442
[6,] -1.60635328 2.75189679
[7,] 2.63857956 -1.60635328
[8,] 8.52242492 2.63857956
[9,] 0.18960081 8.52242492
[10,] 4.43594230 0.18960081
[11,] -1.94717127 4.43594230
[12,] -0.25625930 -1.94717127
[13,] -1.43237181 -0.25625930
[14,] -3.20835776 -1.43237181
[15,] -4.53695592 -3.20835776
[16,] -3.07537244 -4.53695592
[17,] 13.14723296 -3.07537244
[18,] -4.80659431 13.14723296
[19,] -4.91551357 -4.80659431
[20,] -2.74548693 -4.91551357
[21,] 2.30068580 -2.74548693
[22,] 4.56090283 2.30068580
[23,] 2.02530361 4.56090283
[24,] -3.54931588 2.02530361
[25,] 1.27879755 -3.54931588
[26,] 2.55993337 1.27879755
[27,] -1.02736396 2.55993337
[28,] 3.06367458 -1.02736396
[29,] -1.77586550 3.06367458
[30,] 1.08584456 -1.77586550
[31,] -1.64971247 1.08584456
[32,] 0.76982250 -1.64971247
[33,] -2.12082652 0.76982250
[34,] -0.77186463 -2.12082652
[35,] 1.33523158 -0.77186463
[36,] -3.65057272 1.33523158
[37,] 0.96575117 -3.65057272
[38,] -0.23282943 0.96575117
[39,] -5.98078307 -0.23282943
[40,] -4.01807213 -5.98078307
[41,] 4.26327463 -4.01807213
[42,] -4.90318829 4.26327463
[43,] -0.31149938 -4.90318829
[44,] -3.07339884 -0.31149938
[45,] -3.63977739 -3.07339884
[46,] -1.74186435 -3.63977739
[47,] 1.06392164 -1.74186435
[48,] -1.89110334 1.06392164
[49,] 4.38158059 -1.89110334
[50,] 4.20546807 4.38158059
[51,] 0.47272834 4.20546807
[52,] -5.09505885 0.47272834
[53,] 0.58056894 -5.09505885
[54,] 11.10455485 0.58056894
[55,] -3.32124007 11.10455485
[56,] 0.20361648 -3.32124007
[57,] -2.13343358 0.20361648
[58,] 0.74919526 -2.13343358
[59,] 2.29691726 0.74919526
[60,] 1.45852888 2.29691726
[61,] -1.75237796 1.45852888
[62,] -4.29502042 -1.75237796
[63,] 10.15430712 -4.29502042
[64,] -2.16062644 10.15430712
[65,] -4.12832926 -2.16062644
[66,] -1.11554221 -4.12832926
[67,] 5.20623347 -1.11554221
[68,] 1.38397306 5.20623347
[69,] -2.15162616 1.38397306
[70,] -6.08714099 -2.15162616
[71,] 7.21512455 -6.08714099
[72,] 2.40272468 7.21512455
[73,] 0.64744657 2.40272468
[74,] -3.07564356 0.64744657
[75,] 0.46762402 -3.07564356
[76,] 0.50273851 0.46762402
[77,] -2.43559363 0.50273851
[78,] -2.57268587 -2.43559363
[79,] -1.50679206 -2.57268587
[80,] 0.34505745 -1.50679206
[81,] 5.94525104 0.34505745
[82,] 8.61383558 5.94525104
[83,] -3.92718728 8.61383558
[84,] 0.50079958 -3.92718728
[85,] -2.22912285 0.50079958
[86,] -5.17308958 -2.22912285
[87,] 0.09698799 -5.17308958
[88,] -1.97723196 0.09698799
[89,] -1.83407342 -1.97723196
[90,] -0.27212428 -1.83407342
[91,] 0.76862480 -0.27212428
[92,] 4.05095832 0.76862480
[93,] 4.32850101 4.05095832
[94,] -3.47044520 4.32850101
[95,] -2.62965393 -3.47044520
[96,] 1.31502573 -2.62965393
[97,] 1.49107544 1.31502573
[98,] 0.70079206 1.49107544
[99,] -1.58998353 0.70079206
[100,] 1.89294063 -1.58998353
[101,] 1.11794144 1.89294063
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.17113701 -2.30203960
2 -2.21089184 -1.17113701
3 -3.15204128 -2.21089184
4 1.07866442 -3.15204128
5 2.75189679 1.07866442
6 -1.60635328 2.75189679
7 2.63857956 -1.60635328
8 8.52242492 2.63857956
9 0.18960081 8.52242492
10 4.43594230 0.18960081
11 -1.94717127 4.43594230
12 -0.25625930 -1.94717127
13 -1.43237181 -0.25625930
14 -3.20835776 -1.43237181
15 -4.53695592 -3.20835776
16 -3.07537244 -4.53695592
17 13.14723296 -3.07537244
18 -4.80659431 13.14723296
19 -4.91551357 -4.80659431
20 -2.74548693 -4.91551357
21 2.30068580 -2.74548693
22 4.56090283 2.30068580
23 2.02530361 4.56090283
24 -3.54931588 2.02530361
25 1.27879755 -3.54931588
26 2.55993337 1.27879755
27 -1.02736396 2.55993337
28 3.06367458 -1.02736396
29 -1.77586550 3.06367458
30 1.08584456 -1.77586550
31 -1.64971247 1.08584456
32 0.76982250 -1.64971247
33 -2.12082652 0.76982250
34 -0.77186463 -2.12082652
35 1.33523158 -0.77186463
36 -3.65057272 1.33523158
37 0.96575117 -3.65057272
38 -0.23282943 0.96575117
39 -5.98078307 -0.23282943
40 -4.01807213 -5.98078307
41 4.26327463 -4.01807213
42 -4.90318829 4.26327463
43 -0.31149938 -4.90318829
44 -3.07339884 -0.31149938
45 -3.63977739 -3.07339884
46 -1.74186435 -3.63977739
47 1.06392164 -1.74186435
48 -1.89110334 1.06392164
49 4.38158059 -1.89110334
50 4.20546807 4.38158059
51 0.47272834 4.20546807
52 -5.09505885 0.47272834
53 0.58056894 -5.09505885
54 11.10455485 0.58056894
55 -3.32124007 11.10455485
56 0.20361648 -3.32124007
57 -2.13343358 0.20361648
58 0.74919526 -2.13343358
59 2.29691726 0.74919526
60 1.45852888 2.29691726
61 -1.75237796 1.45852888
62 -4.29502042 -1.75237796
63 10.15430712 -4.29502042
64 -2.16062644 10.15430712
65 -4.12832926 -2.16062644
66 -1.11554221 -4.12832926
67 5.20623347 -1.11554221
68 1.38397306 5.20623347
69 -2.15162616 1.38397306
70 -6.08714099 -2.15162616
71 7.21512455 -6.08714099
72 2.40272468 7.21512455
73 0.64744657 2.40272468
74 -3.07564356 0.64744657
75 0.46762402 -3.07564356
76 0.50273851 0.46762402
77 -2.43559363 0.50273851
78 -2.57268587 -2.43559363
79 -1.50679206 -2.57268587
80 0.34505745 -1.50679206
81 5.94525104 0.34505745
82 8.61383558 5.94525104
83 -3.92718728 8.61383558
84 0.50079958 -3.92718728
85 -2.22912285 0.50079958
86 -5.17308958 -2.22912285
87 0.09698799 -5.17308958
88 -1.97723196 0.09698799
89 -1.83407342 -1.97723196
90 -0.27212428 -1.83407342
91 0.76862480 -0.27212428
92 4.05095832 0.76862480
93 4.32850101 4.05095832
94 -3.47044520 4.32850101
95 -2.62965393 -3.47044520
96 1.31502573 -2.62965393
97 1.49107544 1.31502573
98 0.70079206 1.49107544
99 -1.58998353 0.70079206
100 1.89294063 -1.58998353
101 1.11794144 1.89294063
> 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/7l5sa1322131164.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/86w821322131164.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/9sdpu1322131164.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/10fhpb1322131164.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/11ovma1322131164.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/12kvkl1322131164.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/13ysaw1322131164.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/148byt1322131164.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/15hr4o1322131164.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/16ncqg1322131164.tab")
+ }
>
> try(system("convert tmp/1vv5m1322131164.ps tmp/1vv5m1322131164.png",intern=TRUE))
character(0)
> try(system("convert tmp/222ru1322131164.ps tmp/222ru1322131164.png",intern=TRUE))
character(0)
> try(system("convert tmp/3w76a1322131164.ps tmp/3w76a1322131164.png",intern=TRUE))
character(0)
> try(system("convert tmp/4mo9a1322131164.ps tmp/4mo9a1322131164.png",intern=TRUE))
character(0)
> try(system("convert tmp/5037o1322131164.ps tmp/5037o1322131164.png",intern=TRUE))
character(0)
> try(system("convert tmp/6xf111322131164.ps tmp/6xf111322131164.png",intern=TRUE))
character(0)
> try(system("convert tmp/7l5sa1322131164.ps tmp/7l5sa1322131164.png",intern=TRUE))
character(0)
> try(system("convert tmp/86w821322131164.ps tmp/86w821322131164.png",intern=TRUE))
character(0)
> try(system("convert tmp/9sdpu1322131164.ps tmp/9sdpu1322131164.png",intern=TRUE))
character(0)
> try(system("convert tmp/10fhpb1322131164.ps tmp/10fhpb1322131164.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.628 0.532 4.230