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(1966
+ ,1
+ ,41
+ ,1966
+ ,2
+ ,39
+ ,1966
+ ,3
+ ,50
+ ,1966
+ ,4
+ ,40
+ ,1966
+ ,5
+ ,43
+ ,1966
+ ,6
+ ,38
+ ,1966
+ ,7
+ ,44
+ ,1966
+ ,8
+ ,35
+ ,1966
+ ,9
+ ,39
+ ,1966
+ ,10
+ ,35
+ ,1966
+ ,11
+ ,29
+ ,1966
+ ,12
+ ,49
+ ,1967
+ ,1
+ ,50
+ ,1967
+ ,2
+ ,59
+ ,1967
+ ,3
+ ,63
+ ,1967
+ ,4
+ ,32
+ ,1967
+ ,5
+ ,39
+ ,1967
+ ,6
+ ,47
+ ,1967
+ ,7
+ ,53
+ ,1967
+ ,8
+ ,60
+ ,1967
+ ,9
+ ,57
+ ,1967
+ ,10
+ ,52
+ ,1967
+ ,11
+ ,70
+ ,1967
+ ,12
+ ,90
+ ,1968
+ ,1
+ ,74
+ ,1968
+ ,2
+ ,62
+ ,1968
+ ,3
+ ,55
+ ,1968
+ ,4
+ ,84
+ ,1968
+ ,5
+ ,94
+ ,1968
+ ,6
+ ,70
+ ,1968
+ ,7
+ ,108
+ ,1968
+ ,8
+ ,139
+ ,1968
+ ,9
+ ,120
+ ,1968
+ ,10
+ ,97
+ ,1968
+ ,11
+ ,126
+ ,1968
+ ,12
+ ,149
+ ,1969
+ ,1
+ ,158
+ ,1969
+ ,2
+ ,124
+ ,1969
+ ,3
+ ,140
+ ,1969
+ ,4
+ ,109
+ ,1969
+ ,5
+ ,114
+ ,1969
+ ,6
+ ,77
+ ,1969
+ ,7
+ ,120
+ ,1969
+ ,8
+ ,133
+ ,1969
+ ,9
+ ,110
+ ,1969
+ ,10
+ ,92
+ ,1969
+ ,11
+ ,97
+ ,1969
+ ,12
+ ,78
+ ,1970
+ ,1
+ ,99
+ ,1970
+ ,2
+ ,107
+ ,1970
+ ,3
+ ,112
+ ,1970
+ ,4
+ ,90
+ ,1970
+ ,5
+ ,98
+ ,1970
+ ,6
+ ,125
+ ,1970
+ ,7
+ ,155
+ ,1970
+ ,8
+ ,190
+ ,1970
+ ,9
+ ,236
+ ,1970
+ ,10
+ ,189
+ ,1970
+ ,11
+ ,174
+ ,1970
+ ,12
+ ,178
+ ,1971
+ ,1
+ ,136
+ ,1971
+ ,2
+ ,161
+ ,1971
+ ,3
+ ,171
+ ,1971
+ ,4
+ ,149
+ ,1971
+ ,5
+ ,184
+ ,1971
+ ,6
+ ,155
+ ,1971
+ ,7
+ ,276
+ ,1971
+ ,8
+ ,224
+ ,1971
+ ,9
+ ,213
+ ,1971
+ ,10
+ ,279
+ ,1971
+ ,11
+ ,268
+ ,1971
+ ,12
+ ,287
+ ,1972
+ ,1
+ ,238
+ ,1972
+ ,2
+ ,213
+ ,1972
+ ,3
+ ,257
+ ,1972
+ ,4
+ ,293
+ ,1972
+ ,5
+ ,212
+ ,1972
+ ,6
+ ,246
+ ,1972
+ ,7
+ ,353
+ ,1972
+ ,8
+ ,339
+ ,1972
+ ,9
+ ,308
+ ,1972
+ ,10
+ ,247
+ ,1972
+ ,11
+ ,257
+ ,1972
+ ,12
+ ,322
+ ,1973
+ ,1
+ ,298
+ ,1973
+ ,2
+ ,273
+ ,1973
+ ,3
+ ,312
+ ,1973
+ ,4
+ ,249
+ ,1973
+ ,5
+ ,286
+ ,1973
+ ,6
+ ,279
+ ,1973
+ ,7
+ ,309
+ ,1973
+ ,8
+ ,401
+ ,1973
+ ,9
+ ,309
+ ,1973
+ ,10
+ ,328
+ ,1973
+ ,11
+ ,353
+ ,1973
+ ,12
+ ,354
+ ,1974
+ ,1
+ ,327
+ ,1974
+ ,2
+ ,324
+ ,1974
+ ,3
+ ,285
+ ,1974
+ ,4
+ ,243
+ ,1974
+ ,5
+ ,241
+ ,1974
+ ,6
+ ,287
+ ,1974
+ ,7
+ ,355
+ ,1974
+ ,8
+ ,460
+ ,1974
+ ,9
+ ,364
+ ,1974
+ ,10
+ ,487
+ ,1974
+ ,11
+ ,452
+ ,1974
+ ,12
+ ,391
+ ,1975
+ ,1
+ ,500
+ ,1975
+ ,2
+ ,451
+ ,1975
+ ,3
+ ,375
+ ,1975
+ ,4
+ ,372
+ ,1975
+ ,5
+ ,302
+ ,1975
+ ,6
+ ,316
+ ,1975
+ ,7
+ ,398
+ ,1975
+ ,8
+ ,394
+ ,1975
+ ,9
+ ,431
+ ,1975
+ ,10
+ ,431)
+ ,dim=c(3
+ ,118)
+ ,dimnames=list(c('Year'
+ ,'month'
+ ,'robberies')
+ ,1:118))
> y <- array(NA,dim=c(3,118),dimnames=list(c('Year','month','robberies'),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 = '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
robberies Year month
1 41 1966 1
2 39 1966 2
3 50 1966 3
4 40 1966 4
5 43 1966 5
6 38 1966 6
7 44 1966 7
8 35 1966 8
9 39 1966 9
10 35 1966 10
11 29 1966 11
12 49 1966 12
13 50 1967 1
14 59 1967 2
15 63 1967 3
16 32 1967 4
17 39 1967 5
18 47 1967 6
19 53 1967 7
20 60 1967 8
21 57 1967 9
22 52 1967 10
23 70 1967 11
24 90 1967 12
25 74 1968 1
26 62 1968 2
27 55 1968 3
28 84 1968 4
29 94 1968 5
30 70 1968 6
31 108 1968 7
32 139 1968 8
33 120 1968 9
34 97 1968 10
35 126 1968 11
36 149 1968 12
37 158 1969 1
38 124 1969 2
39 140 1969 3
40 109 1969 4
41 114 1969 5
42 77 1969 6
43 120 1969 7
44 133 1969 8
45 110 1969 9
46 92 1969 10
47 97 1969 11
48 78 1969 12
49 99 1970 1
50 107 1970 2
51 112 1970 3
52 90 1970 4
53 98 1970 5
54 125 1970 6
55 155 1970 7
56 190 1970 8
57 236 1970 9
58 189 1970 10
59 174 1970 11
60 178 1970 12
61 136 1971 1
62 161 1971 2
63 171 1971 3
64 149 1971 4
65 184 1971 5
66 155 1971 6
67 276 1971 7
68 224 1971 8
69 213 1971 9
70 279 1971 10
71 268 1971 11
72 287 1971 12
73 238 1972 1
74 213 1972 2
75 257 1972 3
76 293 1972 4
77 212 1972 5
78 246 1972 6
79 353 1972 7
80 339 1972 8
81 308 1972 9
82 247 1972 10
83 257 1972 11
84 322 1972 12
85 298 1973 1
86 273 1973 2
87 312 1973 3
88 249 1973 4
89 286 1973 5
90 279 1973 6
91 309 1973 7
92 401 1973 8
93 309 1973 9
94 328 1973 10
95 353 1973 11
96 354 1973 12
97 327 1974 1
98 324 1974 2
99 285 1974 3
100 243 1974 4
101 241 1974 5
102 287 1974 6
103 355 1974 7
104 460 1974 8
105 364 1974 9
106 487 1974 10
107 452 1974 11
108 391 1974 12
109 500 1975 1
110 451 1975 2
111 375 1975 3
112 372 1975 4
113 302 1975 5
114 316 1975 6
115 398 1975 7
116 394 1975 8
117 431 1975 9
118 431 1975 10
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Year month
-82575.792 41.988 5.802
> (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) -82575.792 2808.959 -29.397 < 2e-16 ***
Year 41.988 1.425 29.457 < 2e-16 ***
month 5.802 1.182 4.907 3.08e-06 ***
---
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/www/rcomp/tmp/1luyv1321954473.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/2x1dj1321954473.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/3vh691321954473.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/4uycv1321954473.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/5eqot1321954473.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/www/rcomp/tmp/6uiac1321954473.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/www/rcomp/tmp/76zka1321954473.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/8kir51321954473.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/9cp1a1321954473.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/10qftm1321954473.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/11x0vr1321954473.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/12y7mc1321954473.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/13yb6t1321954473.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/14tg601321954473.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/15fqb01321954473.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/16826l1321954473.tab")
+ }
>
> try(system("convert tmp/1luyv1321954473.ps tmp/1luyv1321954473.png",intern=TRUE))
character(0)
> try(system("convert tmp/2x1dj1321954473.ps tmp/2x1dj1321954473.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vh691321954473.ps tmp/3vh691321954473.png",intern=TRUE))
character(0)
> try(system("convert tmp/4uycv1321954473.ps tmp/4uycv1321954473.png",intern=TRUE))
character(0)
> try(system("convert tmp/5eqot1321954473.ps tmp/5eqot1321954473.png",intern=TRUE))
character(0)
> try(system("convert tmp/6uiac1321954473.ps tmp/6uiac1321954473.png",intern=TRUE))
character(0)
> try(system("convert tmp/76zka1321954473.ps tmp/76zka1321954473.png",intern=TRUE))
character(0)
> try(system("convert tmp/8kir51321954473.ps tmp/8kir51321954473.png",intern=TRUE))
character(0)
> try(system("convert tmp/9cp1a1321954473.ps tmp/9cp1a1321954473.png",intern=TRUE))
character(0)
> try(system("convert tmp/10qftm1321954473.ps tmp/10qftm1321954473.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.820 0.624 5.670