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,41,2,39,3,50,4,40,5,43,6,38,7,44,8,35,9,39,10,35,11,29,12,49,1,50,2,59,3,63,4,32,5,39,6,47,7,53,8,60,9,57,10,52,11,70,12,90,1,74,2,62,3,55,4,84,5,94,6,70,7,108,8,139,9,120,10,97,11,126,12,149,1,158,2,124,3,140,4,109,5,114,6,77,7,120,8,133,9,110,10,92,11,97,12,78,1,99,2,107,3,112,4,90,5,98,6,125,7,155,8,190,9,236,10,189,11,174,12,178,1,136,2,161,3,171,4,149,5,184,6,155,7,276,8,224,9,213,10,279,11,268,12,287,1,238,2,213,3,257,4,293,5,212,6,246,7,353,8,339,9,308,10,247,11,257,12,322,1,298,2,273,3,312,4,249,5,286,6,279,7,309,8,401,9,309,10,328,11,353,12,354,1,327,2,324,3,285,4,243,5,241,6,287,7,355,8,460,9,364,10,487,11,452,12,391,1,500,2,451,3,375,4,372,5,302,6,316,7,398,8,394,9,431,10,431),dim=c(2,118),dimnames=list(c('Colombia','USA'),1:118))
> y <- array(NA,dim=c(2,118),dimnames=list(c('Colombia','USA'),1:118))
> 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'
> 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
USA Colombia t
1 41 1 1
2 39 2 2
3 50 3 3
4 40 4 4
5 43 5 5
6 38 6 6
7 44 7 7
8 35 8 8
9 39 9 9
10 35 10 10
11 29 11 11
12 49 12 12
13 50 1 13
14 59 2 14
15 63 3 15
16 32 4 16
17 39 5 17
18 47 6 18
19 53 7 19
20 60 8 20
21 57 9 21
22 52 10 22
23 70 11 23
24 90 12 24
25 74 1 25
26 62 2 26
27 55 3 27
28 84 4 28
29 94 5 29
30 70 6 30
31 108 7 31
32 139 8 32
33 120 9 33
34 97 10 34
35 126 11 35
36 149 12 36
37 158 1 37
38 124 2 38
39 140 3 39
40 109 4 40
41 114 5 41
42 77 6 42
43 120 7 43
44 133 8 44
45 110 9 45
46 92 10 46
47 97 11 47
48 78 12 48
49 99 1 49
50 107 2 50
51 112 3 51
52 90 4 52
53 98 5 53
54 125 6 54
55 155 7 55
56 190 8 56
57 236 9 57
58 189 10 58
59 174 11 59
60 178 12 60
61 136 1 61
62 161 2 62
63 171 3 63
64 149 4 64
65 184 5 65
66 155 6 66
67 276 7 67
68 224 8 68
69 213 9 69
70 279 10 70
71 268 11 71
72 287 12 72
73 238 1 73
74 213 2 74
75 257 3 75
76 293 4 76
77 212 5 77
78 246 6 78
79 353 7 79
80 339 8 80
81 308 9 81
82 247 10 82
83 257 11 83
84 322 12 84
85 298 1 85
86 273 2 86
87 312 3 87
88 249 4 88
89 286 5 89
90 279 6 90
91 309 7 91
92 401 8 92
93 309 9 93
94 328 10 94
95 353 11 95
96 354 12 96
97 327 1 97
98 324 2 98
99 285 3 99
100 243 4 100
101 241 5 101
102 287 6 102
103 355 7 103
104 460 8 104
105 364 9 105
106 487 10 106
107 452 11 107
108 391 12 108
109 500 1 109
110 451 2 110
111 375 3 111
112 372 4 112
113 302 5 113
114 316 6 114
115 398 7 115
116 394 8 116
117 431 9 117
118 431 10 118
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Colombia t
-26.680 2.303 3.499
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-97.238 -26.131 1.727 24.547 142.983
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -26.6803 10.8373 -2.462 0.0153 *
Colombia 2.3033 1.1837 1.946 0.0541 .
t 3.4990 0.1188 29.457 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 43.87 on 115 degrees of freedom
Multiple R-squared: 0.8846, Adjusted R-squared: 0.8826
F-statistic: 440.8 on 2 and 115 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,] 2.924770e-03 5.849541e-03 0.9970752
[2,] 3.531810e-04 7.063620e-04 0.9996468
[3,] 6.654536e-05 1.330907e-04 0.9999335
[4,] 6.859359e-06 1.371872e-05 0.9999931
[5,] 7.744253e-07 1.548851e-06 0.9999992
[6,] 1.819752e-07 3.639505e-07 0.9999998
[7,] 4.367093e-07 8.734186e-07 0.9999996
[8,] 6.050300e-08 1.210060e-07 0.9999999
[9,] 1.553539e-08 3.107078e-08 1.0000000
[10,] 4.626215e-09 9.252429e-09 1.0000000
[11,] 4.666063e-08 9.332126e-08 1.0000000
[12,] 1.539979e-08 3.079957e-08 1.0000000
[13,] 2.794030e-09 5.588060e-09 1.0000000
[14,] 6.755006e-10 1.351001e-09 1.0000000
[15,] 3.477173e-10 6.954345e-10 1.0000000
[16,] 9.139455e-11 1.827891e-10 1.0000000
[17,] 1.633058e-11 3.266116e-11 1.0000000
[18,] 2.468954e-11 4.937908e-11 1.0000000
[19,] 7.693622e-10 1.538724e-09 1.0000000
[20,] 2.221252e-10 4.442505e-10 1.0000000
[21,] 6.372120e-11 1.274424e-10 1.0000000
[22,] 2.763950e-11 5.527901e-11 1.0000000
[23,] 1.686252e-11 3.372504e-11 1.0000000
[24,] 2.300929e-11 4.601858e-11 1.0000000
[25,] 5.998004e-12 1.199601e-11 1.0000000
[26,] 3.162103e-11 6.324206e-11 1.0000000
[27,] 7.742455e-09 1.548491e-08 1.0000000
[28,] 1.050578e-08 2.101155e-08 1.0000000
[29,] 3.612916e-09 7.225832e-09 1.0000000
[30,] 3.827674e-09 7.655348e-09 1.0000000
[31,] 1.786754e-08 3.573508e-08 1.0000000
[32,] 1.500155e-07 3.000309e-07 0.9999998
[33,] 7.939029e-08 1.587806e-07 0.9999999
[34,] 6.232236e-08 1.246447e-07 0.9999999
[35,] 3.778251e-08 7.556501e-08 1.0000000
[36,] 1.966034e-08 3.932069e-08 1.0000000
[37,] 9.421098e-08 1.884220e-07 0.9999999
[38,] 4.224500e-08 8.449000e-08 1.0000000
[39,] 1.886112e-08 3.772225e-08 1.0000000
[40,] 1.129918e-08 2.259835e-08 1.0000000
[41,] 1.894643e-08 3.789286e-08 1.0000000
[42,] 2.138338e-08 4.276677e-08 1.0000000
[43,] 9.756248e-08 1.951250e-07 0.9999999
[44,] 1.036605e-07 2.073210e-07 0.9999999
[45,] 6.911363e-08 1.382273e-07 0.9999999
[46,] 3.897925e-08 7.795851e-08 1.0000000
[47,] 6.349989e-08 1.269998e-07 0.9999999
[48,] 6.932503e-08 1.386501e-07 0.9999999
[49,] 3.659812e-08 7.319625e-08 1.0000000
[50,] 2.610789e-08 5.221577e-08 1.0000000
[51,] 8.686387e-08 1.737277e-07 0.9999999
[52,] 4.643121e-06 9.286241e-06 0.9999954
[53,] 4.519190e-06 9.038379e-06 0.9999955
[54,] 3.002124e-06 6.004248e-06 0.9999970
[55,] 2.110411e-06 4.220821e-06 0.9999979
[56,] 1.389472e-06 2.778944e-06 0.9999986
[57,] 7.395612e-07 1.479122e-06 0.9999993
[58,] 4.063160e-07 8.126320e-07 0.9999996
[59,] 3.114849e-07 6.229699e-07 0.9999997
[60,] 1.967694e-07 3.935388e-07 0.9999998
[61,] 2.068460e-07 4.136921e-07 0.9999998
[62,] 6.585218e-06 1.317044e-05 0.9999934
[63,] 6.309159e-06 1.261832e-05 0.9999937
[64,] 4.937418e-06 9.874836e-06 0.9999951
[65,] 2.098868e-05 4.197735e-05 0.9999790
[66,] 3.221842e-05 6.443685e-05 0.9999678
[67,] 6.275475e-05 1.255095e-04 0.9999372
[68,] 5.056412e-05 1.011282e-04 0.9999494
[69,] 3.109508e-05 6.219016e-05 0.9999689
[70,] 2.748912e-05 5.497824e-05 0.9999725
[71,] 5.776520e-05 1.155304e-04 0.9999422
[72,] 4.454678e-05 8.909356e-05 0.9999555
[73,] 2.809271e-05 5.618543e-05 0.9999719
[74,] 3.456202e-04 6.912404e-04 0.9996544
[75,] 1.064065e-03 2.128130e-03 0.9989359
[76,] 1.029434e-03 2.058868e-03 0.9989706
[77,] 7.809154e-04 1.561831e-03 0.9992191
[78,] 6.048129e-04 1.209626e-03 0.9993952
[79,] 5.033144e-04 1.006629e-03 0.9994967
[80,] 4.322185e-04 8.644371e-04 0.9995678
[81,] 2.582928e-04 5.165855e-04 0.9997417
[82,] 2.282520e-04 4.565039e-04 0.9997717
[83,] 1.659766e-04 3.319531e-04 0.9998340
[84,] 9.443668e-05 1.888734e-04 0.9999056
[85,] 5.747119e-05 1.149424e-04 0.9999425
[86,] 3.258841e-05 6.517682e-05 0.9999674
[87,] 1.835662e-04 3.671324e-04 0.9998164
[88,] 1.040308e-04 2.080616e-04 0.9998960
[89,] 5.798478e-05 1.159696e-04 0.9999420
[90,] 3.602603e-05 7.205207e-05 0.9999640
[91,] 2.052445e-05 4.104890e-05 0.9999795
[92,] 1.213746e-05 2.427492e-05 0.9999879
[93,] 6.398751e-06 1.279750e-05 0.9999936
[94,] 3.967627e-06 7.935255e-06 0.9999960
[95,] 1.952724e-05 3.905448e-05 0.9999805
[96,] 4.433698e-04 8.867397e-04 0.9995566
[97,] 3.884527e-03 7.769054e-03 0.9961155
[98,] 6.898615e-03 1.379723e-02 0.9931014
[99,] 9.612355e-03 1.922471e-02 0.9903876
[100,] 1.598123e-02 3.196246e-02 0.9840188
[101,] 2.723854e-02 5.447708e-02 0.9727615
[102,] 2.667983e-02 5.335966e-02 0.9733202
[103,] 1.491099e-02 2.982198e-02 0.9850890
[104,] 1.173664e-01 2.347327e-01 0.8826336
[105,] 4.408230e-01 8.816460e-01 0.5591770
[106,] 4.995004e-01 9.990008e-01 0.5004996
[107,] 7.901001e-01 4.197997e-01 0.2098999
> postscript(file="/var/wessaorg/rcomp/tmp/1or3i1322055808.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/2k0br1322055808.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/3ihd91322055808.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/4nzqr1322055808.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/5217j1322055808.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 = 118
Frequency = 1
1 2 3 4 5 6 7
61.878007 54.075688 59.273368 43.471049 40.668729 29.866410 30.064090
8 9 10 11 12 13 14
15.261771 13.459451 3.657132 -8.145187 6.052493 28.889649 32.087330
15 16 17 18 19 20 21
30.285011 -6.517309 -5.319628 -3.121948 -2.924267 -1.726587 -10.528906
22 23 24 25 26 27 28
-21.331226 -9.133545 5.064135 10.901292 -6.901028 -19.703347 3.494333
29 30 31 32 33 34 35
7.692014 -22.110306 10.087375 35.285056 10.482736 -18.319583 4.878097
36 37 38 39 40 41 42
22.075778 52.912934 13.110615 23.308295 -13.494024 -14.296344 -57.098663
43 44 45 46 47 48 49
-19.900983 -12.703302 -41.505622 -65.307941 -66.110260 -90.912580 -48.075424
50 51 52 53 54 55 56
-45.877743 -46.680063 -74.482382 -72.284701 -51.087021 -26.889340 2.308340
57 58 59 60 61 62 63
42.506021 -10.296299 -31.098618 -32.900938 -53.063781 -33.866101 -29.668420
64 65 66 67 68 69 70
-57.470740 -28.273059 -63.075379 52.122302 -5.680017 -22.482337 37.715344
71 72 73 74 75 76 77
20.913024 34.110705 6.947861 -23.854458 14.343222 44.540903 -42.261417
78 79 80 81 82 83 84
-14.063736 87.133944 67.331625 30.529305 -36.273014 -32.075333 27.122347
85 86 87 88 89 90 91
24.959503 -5.842816 27.354864 -41.447455 -10.249774 -23.052094 1.145587
92 93 94 95 96 97 98
87.343267 -10.459052 2.738628 21.936309 17.133989 11.971146 3.168826
99 100 101 102 103 104 105
-41.633493 -89.435813 -97.238132 -57.040452 5.157229 104.354910 2.552590
106 107 108 109 110 111 112
119.750271 78.947951 12.145632 142.982788 88.180469 6.378149 -2.424170
113 114 115 116 117 118
-78.226490 -70.028809 6.168871 -3.633448 27.564232 21.761913
> postscript(file="/var/wessaorg/rcomp/tmp/6lh831322055808.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 = 118
Frequency = 1
lag(myerror, k = 1) myerror
0 61.878007 NA
1 54.075688 61.878007
2 59.273368 54.075688
3 43.471049 59.273368
4 40.668729 43.471049
5 29.866410 40.668729
6 30.064090 29.866410
7 15.261771 30.064090
8 13.459451 15.261771
9 3.657132 13.459451
10 -8.145187 3.657132
11 6.052493 -8.145187
12 28.889649 6.052493
13 32.087330 28.889649
14 30.285011 32.087330
15 -6.517309 30.285011
16 -5.319628 -6.517309
17 -3.121948 -5.319628
18 -2.924267 -3.121948
19 -1.726587 -2.924267
20 -10.528906 -1.726587
21 -21.331226 -10.528906
22 -9.133545 -21.331226
23 5.064135 -9.133545
24 10.901292 5.064135
25 -6.901028 10.901292
26 -19.703347 -6.901028
27 3.494333 -19.703347
28 7.692014 3.494333
29 -22.110306 7.692014
30 10.087375 -22.110306
31 35.285056 10.087375
32 10.482736 35.285056
33 -18.319583 10.482736
34 4.878097 -18.319583
35 22.075778 4.878097
36 52.912934 22.075778
37 13.110615 52.912934
38 23.308295 13.110615
39 -13.494024 23.308295
40 -14.296344 -13.494024
41 -57.098663 -14.296344
42 -19.900983 -57.098663
43 -12.703302 -19.900983
44 -41.505622 -12.703302
45 -65.307941 -41.505622
46 -66.110260 -65.307941
47 -90.912580 -66.110260
48 -48.075424 -90.912580
49 -45.877743 -48.075424
50 -46.680063 -45.877743
51 -74.482382 -46.680063
52 -72.284701 -74.482382
53 -51.087021 -72.284701
54 -26.889340 -51.087021
55 2.308340 -26.889340
56 42.506021 2.308340
57 -10.296299 42.506021
58 -31.098618 -10.296299
59 -32.900938 -31.098618
60 -53.063781 -32.900938
61 -33.866101 -53.063781
62 -29.668420 -33.866101
63 -57.470740 -29.668420
64 -28.273059 -57.470740
65 -63.075379 -28.273059
66 52.122302 -63.075379
67 -5.680017 52.122302
68 -22.482337 -5.680017
69 37.715344 -22.482337
70 20.913024 37.715344
71 34.110705 20.913024
72 6.947861 34.110705
73 -23.854458 6.947861
74 14.343222 -23.854458
75 44.540903 14.343222
76 -42.261417 44.540903
77 -14.063736 -42.261417
78 87.133944 -14.063736
79 67.331625 87.133944
80 30.529305 67.331625
81 -36.273014 30.529305
82 -32.075333 -36.273014
83 27.122347 -32.075333
84 24.959503 27.122347
85 -5.842816 24.959503
86 27.354864 -5.842816
87 -41.447455 27.354864
88 -10.249774 -41.447455
89 -23.052094 -10.249774
90 1.145587 -23.052094
91 87.343267 1.145587
92 -10.459052 87.343267
93 2.738628 -10.459052
94 21.936309 2.738628
95 17.133989 21.936309
96 11.971146 17.133989
97 3.168826 11.971146
98 -41.633493 3.168826
99 -89.435813 -41.633493
100 -97.238132 -89.435813
101 -57.040452 -97.238132
102 5.157229 -57.040452
103 104.354910 5.157229
104 2.552590 104.354910
105 119.750271 2.552590
106 78.947951 119.750271
107 12.145632 78.947951
108 142.982788 12.145632
109 88.180469 142.982788
110 6.378149 88.180469
111 -2.424170 6.378149
112 -78.226490 -2.424170
113 -70.028809 -78.226490
114 6.168871 -70.028809
115 -3.633448 6.168871
116 27.564232 -3.633448
117 21.761913 27.564232
118 NA 21.761913
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 54.075688 61.878007
[2,] 59.273368 54.075688
[3,] 43.471049 59.273368
[4,] 40.668729 43.471049
[5,] 29.866410 40.668729
[6,] 30.064090 29.866410
[7,] 15.261771 30.064090
[8,] 13.459451 15.261771
[9,] 3.657132 13.459451
[10,] -8.145187 3.657132
[11,] 6.052493 -8.145187
[12,] 28.889649 6.052493
[13,] 32.087330 28.889649
[14,] 30.285011 32.087330
[15,] -6.517309 30.285011
[16,] -5.319628 -6.517309
[17,] -3.121948 -5.319628
[18,] -2.924267 -3.121948
[19,] -1.726587 -2.924267
[20,] -10.528906 -1.726587
[21,] -21.331226 -10.528906
[22,] -9.133545 -21.331226
[23,] 5.064135 -9.133545
[24,] 10.901292 5.064135
[25,] -6.901028 10.901292
[26,] -19.703347 -6.901028
[27,] 3.494333 -19.703347
[28,] 7.692014 3.494333
[29,] -22.110306 7.692014
[30,] 10.087375 -22.110306
[31,] 35.285056 10.087375
[32,] 10.482736 35.285056
[33,] -18.319583 10.482736
[34,] 4.878097 -18.319583
[35,] 22.075778 4.878097
[36,] 52.912934 22.075778
[37,] 13.110615 52.912934
[38,] 23.308295 13.110615
[39,] -13.494024 23.308295
[40,] -14.296344 -13.494024
[41,] -57.098663 -14.296344
[42,] -19.900983 -57.098663
[43,] -12.703302 -19.900983
[44,] -41.505622 -12.703302
[45,] -65.307941 -41.505622
[46,] -66.110260 -65.307941
[47,] -90.912580 -66.110260
[48,] -48.075424 -90.912580
[49,] -45.877743 -48.075424
[50,] -46.680063 -45.877743
[51,] -74.482382 -46.680063
[52,] -72.284701 -74.482382
[53,] -51.087021 -72.284701
[54,] -26.889340 -51.087021
[55,] 2.308340 -26.889340
[56,] 42.506021 2.308340
[57,] -10.296299 42.506021
[58,] -31.098618 -10.296299
[59,] -32.900938 -31.098618
[60,] -53.063781 -32.900938
[61,] -33.866101 -53.063781
[62,] -29.668420 -33.866101
[63,] -57.470740 -29.668420
[64,] -28.273059 -57.470740
[65,] -63.075379 -28.273059
[66,] 52.122302 -63.075379
[67,] -5.680017 52.122302
[68,] -22.482337 -5.680017
[69,] 37.715344 -22.482337
[70,] 20.913024 37.715344
[71,] 34.110705 20.913024
[72,] 6.947861 34.110705
[73,] -23.854458 6.947861
[74,] 14.343222 -23.854458
[75,] 44.540903 14.343222
[76,] -42.261417 44.540903
[77,] -14.063736 -42.261417
[78,] 87.133944 -14.063736
[79,] 67.331625 87.133944
[80,] 30.529305 67.331625
[81,] -36.273014 30.529305
[82,] -32.075333 -36.273014
[83,] 27.122347 -32.075333
[84,] 24.959503 27.122347
[85,] -5.842816 24.959503
[86,] 27.354864 -5.842816
[87,] -41.447455 27.354864
[88,] -10.249774 -41.447455
[89,] -23.052094 -10.249774
[90,] 1.145587 -23.052094
[91,] 87.343267 1.145587
[92,] -10.459052 87.343267
[93,] 2.738628 -10.459052
[94,] 21.936309 2.738628
[95,] 17.133989 21.936309
[96,] 11.971146 17.133989
[97,] 3.168826 11.971146
[98,] -41.633493 3.168826
[99,] -89.435813 -41.633493
[100,] -97.238132 -89.435813
[101,] -57.040452 -97.238132
[102,] 5.157229 -57.040452
[103,] 104.354910 5.157229
[104,] 2.552590 104.354910
[105,] 119.750271 2.552590
[106,] 78.947951 119.750271
[107,] 12.145632 78.947951
[108,] 142.982788 12.145632
[109,] 88.180469 142.982788
[110,] 6.378149 88.180469
[111,] -2.424170 6.378149
[112,] -78.226490 -2.424170
[113,] -70.028809 -78.226490
[114,] 6.168871 -70.028809
[115,] -3.633448 6.168871
[116,] 27.564232 -3.633448
[117,] 21.761913 27.564232
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 54.075688 61.878007
2 59.273368 54.075688
3 43.471049 59.273368
4 40.668729 43.471049
5 29.866410 40.668729
6 30.064090 29.866410
7 15.261771 30.064090
8 13.459451 15.261771
9 3.657132 13.459451
10 -8.145187 3.657132
11 6.052493 -8.145187
12 28.889649 6.052493
13 32.087330 28.889649
14 30.285011 32.087330
15 -6.517309 30.285011
16 -5.319628 -6.517309
17 -3.121948 -5.319628
18 -2.924267 -3.121948
19 -1.726587 -2.924267
20 -10.528906 -1.726587
21 -21.331226 -10.528906
22 -9.133545 -21.331226
23 5.064135 -9.133545
24 10.901292 5.064135
25 -6.901028 10.901292
26 -19.703347 -6.901028
27 3.494333 -19.703347
28 7.692014 3.494333
29 -22.110306 7.692014
30 10.087375 -22.110306
31 35.285056 10.087375
32 10.482736 35.285056
33 -18.319583 10.482736
34 4.878097 -18.319583
35 22.075778 4.878097
36 52.912934 22.075778
37 13.110615 52.912934
38 23.308295 13.110615
39 -13.494024 23.308295
40 -14.296344 -13.494024
41 -57.098663 -14.296344
42 -19.900983 -57.098663
43 -12.703302 -19.900983
44 -41.505622 -12.703302
45 -65.307941 -41.505622
46 -66.110260 -65.307941
47 -90.912580 -66.110260
48 -48.075424 -90.912580
49 -45.877743 -48.075424
50 -46.680063 -45.877743
51 -74.482382 -46.680063
52 -72.284701 -74.482382
53 -51.087021 -72.284701
54 -26.889340 -51.087021
55 2.308340 -26.889340
56 42.506021 2.308340
57 -10.296299 42.506021
58 -31.098618 -10.296299
59 -32.900938 -31.098618
60 -53.063781 -32.900938
61 -33.866101 -53.063781
62 -29.668420 -33.866101
63 -57.470740 -29.668420
64 -28.273059 -57.470740
65 -63.075379 -28.273059
66 52.122302 -63.075379
67 -5.680017 52.122302
68 -22.482337 -5.680017
69 37.715344 -22.482337
70 20.913024 37.715344
71 34.110705 20.913024
72 6.947861 34.110705
73 -23.854458 6.947861
74 14.343222 -23.854458
75 44.540903 14.343222
76 -42.261417 44.540903
77 -14.063736 -42.261417
78 87.133944 -14.063736
79 67.331625 87.133944
80 30.529305 67.331625
81 -36.273014 30.529305
82 -32.075333 -36.273014
83 27.122347 -32.075333
84 24.959503 27.122347
85 -5.842816 24.959503
86 27.354864 -5.842816
87 -41.447455 27.354864
88 -10.249774 -41.447455
89 -23.052094 -10.249774
90 1.145587 -23.052094
91 87.343267 1.145587
92 -10.459052 87.343267
93 2.738628 -10.459052
94 21.936309 2.738628
95 17.133989 21.936309
96 11.971146 17.133989
97 3.168826 11.971146
98 -41.633493 3.168826
99 -89.435813 -41.633493
100 -97.238132 -89.435813
101 -57.040452 -97.238132
102 5.157229 -57.040452
103 104.354910 5.157229
104 2.552590 104.354910
105 119.750271 2.552590
106 78.947951 119.750271
107 12.145632 78.947951
108 142.982788 12.145632
109 88.180469 142.982788
110 6.378149 88.180469
111 -2.424170 6.378149
112 -78.226490 -2.424170
113 -70.028809 -78.226490
114 6.168871 -70.028809
115 -3.633448 6.168871
116 27.564232 -3.633448
117 21.761913 27.564232
> 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/7z3dt1322055808.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/83l5r1322055808.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/9d99y1322055808.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/109lqs1322055808.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/11232q1322055808.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/123sh61322055808.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/13tppe1322055808.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/146yi51322055808.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/152vz51322055808.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/16h4xu1322055808.tab")
+ }
>
> try(system("convert tmp/1or3i1322055808.ps tmp/1or3i1322055808.png",intern=TRUE))
character(0)
> try(system("convert tmp/2k0br1322055808.ps tmp/2k0br1322055808.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ihd91322055808.ps tmp/3ihd91322055808.png",intern=TRUE))
character(0)
> try(system("convert tmp/4nzqr1322055808.ps tmp/4nzqr1322055808.png",intern=TRUE))
character(0)
> try(system("convert tmp/5217j1322055808.ps tmp/5217j1322055808.png",intern=TRUE))
character(0)
> try(system("convert tmp/6lh831322055808.ps tmp/6lh831322055808.png",intern=TRUE))
character(0)
> try(system("convert tmp/7z3dt1322055808.ps tmp/7z3dt1322055808.png",intern=TRUE))
character(0)
> try(system("convert tmp/83l5r1322055808.ps tmp/83l5r1322055808.png",intern=TRUE))
character(0)
> try(system("convert tmp/9d99y1322055808.ps tmp/9d99y1322055808.png",intern=TRUE))
character(0)
> try(system("convert tmp/109lqs1322055808.ps tmp/109lqs1322055808.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.064 0.553 4.663