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(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 = '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
month Year robberies
1 1 1966 41
2 2 1966 39
3 3 1966 50
4 4 1966 40
5 5 1966 43
6 6 1966 38
7 7 1966 44
8 8 1966 35
9 9 1966 39
10 10 1966 35
11 11 1966 29
12 12 1966 49
13 1 1967 50
14 2 1967 59
15 3 1967 63
16 4 1967 32
17 5 1967 39
18 6 1967 47
19 7 1967 53
20 8 1967 60
21 9 1967 57
22 10 1967 52
23 11 1967 70
24 12 1967 90
25 1 1968 74
26 2 1968 62
27 3 1968 55
28 4 1968 84
29 5 1968 94
30 6 1968 70
31 7 1968 108
32 8 1968 139
33 9 1968 120
34 10 1968 97
35 11 1968 126
36 12 1968 149
37 1 1969 158
38 2 1969 124
39 3 1969 140
40 4 1969 109
41 5 1969 114
42 6 1969 77
43 7 1969 120
44 8 1969 133
45 9 1969 110
46 10 1969 92
47 11 1969 97
48 12 1969 78
49 1 1970 99
50 2 1970 107
51 3 1970 112
52 4 1970 90
53 5 1970 98
54 6 1970 125
55 7 1970 155
56 8 1970 190
57 9 1970 236
58 10 1970 189
59 11 1970 174
60 12 1970 178
61 1 1971 136
62 2 1971 161
63 3 1971 171
64 4 1971 149
65 5 1971 184
66 6 1971 155
67 7 1971 276
68 8 1971 224
69 9 1971 213
70 10 1971 279
71 11 1971 268
72 12 1971 287
73 1 1972 238
74 2 1972 213
75 3 1972 257
76 4 1972 293
77 5 1972 212
78 6 1972 246
79 7 1972 353
80 8 1972 339
81 9 1972 308
82 10 1972 247
83 11 1972 257
84 12 1972 322
85 1 1973 298
86 2 1973 273
87 3 1973 312
88 4 1973 249
89 5 1973 286
90 6 1973 279
91 7 1973 309
92 8 1973 401
93 9 1973 309
94 10 1973 328
95 11 1973 353
96 12 1973 354
97 1 1974 327
98 2 1974 324
99 3 1974 285
100 4 1974 243
101 5 1974 241
102 6 1974 287
103 7 1974 355
104 8 1974 460
105 9 1974 364
106 10 1974 487
107 11 1974 452
108 12 1974 391
109 1 1975 500
110 2 1975 451
111 3 1975 375
112 4 1975 372
113 5 1975 302
114 6 1975 316
115 7 1975 398
116 8 1975 394
117 9 1975 431
118 10 1975 431
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Year robberies
2547.65328 -1.29266 0.02984
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8.5613 -2.2737 0.3269 2.1349 7.2736
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.548e+03 5.376e+02 4.739 6.20e-06 ***
Year -1.293e+00 2.734e-01 -4.728 6.48e-06 ***
robberies 2.984e-02 6.081e-03 4.907 3.08e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.146 on 115 degrees of freedom
Multiple R-squared: 0.1744, Adjusted R-squared: 0.1601
F-statistic: 12.15 on 2 and 115 DF, p-value: 1.633e-05
> 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.3964920 0.79298400 0.60350800
[2,] 0.4454164 0.89083281 0.55458359
[3,] 0.4537599 0.90751982 0.54624009
[4,] 0.5200594 0.95988118 0.47994059
[5,] 0.5164996 0.96700080 0.48350040
[6,] 0.4275999 0.85519979 0.57240010
[7,] 0.8608077 0.27838457 0.13919228
[8,] 0.8260021 0.34799583 0.17399791
[9,] 0.8004451 0.39910979 0.19955489
[10,] 0.7788150 0.44237007 0.22118504
[11,] 0.7211957 0.55760862 0.27880431
[12,] 0.6627741 0.67445187 0.33722593
[13,] 0.6353418 0.72931640 0.36465820
[14,] 0.6500142 0.69997170 0.34998585
[15,] 0.6984356 0.60312873 0.30156437
[16,] 0.7385296 0.52294074 0.26147037
[17,] 0.7843617 0.43127652 0.21563826
[18,] 0.8341541 0.33169186 0.16584593
[19,] 0.8334142 0.33317157 0.16658579
[20,] 0.8820432 0.23591363 0.11795681
[21,] 0.8713138 0.25737234 0.12868617
[22,] 0.8458975 0.30820501 0.15410250
[23,] 0.8241311 0.35173782 0.17586891
[24,] 0.7938121 0.41237589 0.20618794
[25,] 0.7625670 0.47486594 0.23743297
[26,] 0.7155447 0.56891060 0.28445530
[27,] 0.6706203 0.65875939 0.32937970
[28,] 0.6212256 0.75754874 0.37877437
[29,] 0.6263392 0.74732169 0.37366085
[30,] 0.5930890 0.81382196 0.40691098
[31,] 0.5482036 0.90359276 0.45179638
[32,] 0.7932628 0.41347430 0.20673715
[33,] 0.8211669 0.35766610 0.17883305
[34,] 0.8419494 0.31610113 0.15805056
[35,] 0.8242849 0.35143018 0.17571509
[36,] 0.8013447 0.39731056 0.19865528
[37,] 0.8057641 0.38847172 0.19423586
[38,] 0.7789896 0.44202086 0.22101043
[39,] 0.7494134 0.50117319 0.25058660
[40,] 0.7577240 0.48455194 0.24227597
[41,] 0.8123792 0.37524162 0.18762081
[42,] 0.8689721 0.26205577 0.13102789
[43,] 0.9418843 0.11623135 0.05811568
[44,] 0.9487507 0.10249868 0.05124934
[45,] 0.9470150 0.10596996 0.05298498
[46,] 0.9388025 0.12239504 0.06119752
[47,] 0.9218635 0.15627297 0.07813649
[48,] 0.9013833 0.19723331 0.09861665
[49,] 0.8768728 0.24625432 0.12312716
[50,] 0.8479024 0.30419529 0.15209764
[51,] 0.8141980 0.37160401 0.18580200
[52,] 0.7776015 0.44479709 0.22239855
[53,] 0.7533540 0.49329190 0.24664595
[54,] 0.7646141 0.47077183 0.23538592
[55,] 0.8053482 0.38930353 0.19465176
[56,] 0.8134107 0.37317863 0.18658931
[57,] 0.8201108 0.35977847 0.17988924
[58,] 0.8131971 0.37360570 0.18680285
[59,] 0.7798439 0.44031221 0.22015611
[60,] 0.7460583 0.50788338 0.25394169
[61,] 0.7035397 0.59292055 0.29646027
[62,] 0.6889750 0.62204992 0.31102496
[63,] 0.6396711 0.72065774 0.36032887
[64,] 0.6028831 0.79423380 0.39711690
[65,] 0.5525950 0.89481002 0.44740501
[66,] 0.5221503 0.95569944 0.47784972
[67,] 0.5089945 0.98201110 0.49100555
[68,] 0.6148032 0.77039360 0.38519680
[69,] 0.6334621 0.73307582 0.36653791
[70,] 0.6713926 0.65721489 0.32860744
[71,] 0.7139674 0.57206511 0.28603255
[72,] 0.6721051 0.65578985 0.32789492
[73,] 0.6268452 0.74630959 0.37315480
[74,] 0.6266212 0.74675764 0.37337882
[75,] 0.5906405 0.81871908 0.40935954
[76,] 0.5348407 0.93031869 0.46515935
[77,] 0.5213173 0.95736538 0.47868269
[78,] 0.5401430 0.91971390 0.45985695
[79,] 0.5387427 0.92251468 0.46125734
[80,] 0.6685245 0.66295090 0.33147545
[81,] 0.7143335 0.57133308 0.28566654
[82,] 0.7714362 0.45712769 0.22856385
[83,] 0.7508748 0.49825049 0.24912524
[84,] 0.7318276 0.53634479 0.26817240
[85,] 0.6922778 0.61544447 0.30772224
[86,] 0.6458128 0.70837436 0.35418718
[87,] 0.6260459 0.74790815 0.37395407
[88,] 0.5674451 0.86510982 0.43255491
[89,] 0.5095369 0.98092612 0.49046306
[90,] 0.4578290 0.91565804 0.54217098
[91,] 0.4528356 0.90567122 0.54716439
[92,] 0.5841191 0.83176173 0.41588086
[93,] 0.6740222 0.65195566 0.32597783
[94,] 0.6947735 0.61045304 0.30522652
[95,] 0.6796182 0.64076355 0.32038178
[96,] 0.6653889 0.66922218 0.33461109
[97,] 0.6685497 0.66290055 0.33145027
[98,] 0.6688318 0.66233641 0.33116820
[99,] 0.6245118 0.75097645 0.37548822
[100,] 0.5845416 0.83091683 0.41545842
[101,] 0.4917580 0.98351590 0.50824205
[102,] 0.3906362 0.78127250 0.60936375
[103,] 0.3009976 0.60199527 0.69900236
[104,] 0.4647881 0.92957625 0.53521188
[105,] 0.7907118 0.41857637 0.20928818
[106,] 0.8918322 0.21633569 0.10816785
[107,] 0.9750005 0.04999901 0.02499950
> postscript(file="/var/wessaorg/rcomp/tmp/1owrs1321954189.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/29k1y1321954189.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/3ymbh1321954189.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/4jbh51321954189.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/55ovd1321954189.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
-6.50040646 -5.44073365 -4.76893412 -3.47057005 -2.56007927 -1.41089724
7 8 9 10 11 12
-0.58991568 0.67861198 1.55926635 2.67861198 3.85763042 4.26090229
13 14 15 16 17 18
-5.47627076 -4.74479842 -3.86414404 -1.93921544 -1.14807029 -0.38676154
19 20 21 22 23 24
0.43422002 1.22536518 2.31487440 3.46405643 3.92700111 4.33027299
25 26 27 28 29 30
-4.89968115 -3.54164428 -2.33278943 -2.19804522 -1.49640928 0.21966447
31 32 33 34 35 36
0.08588103 0.16095243 1.72784415 3.41408150 3.54882571 3.86258837
37 38 39 40 41 42
-6.11327593 -4.09883812 -3.57622062 -1.65129202 -0.80047405 1.30347299
43 44 45 46 47 48
1.02050751 1.63263423 3.31887157 4.85592689 5.70674486 7.27363658
49 50 51 52 53 54
-3.06026460 -2.29895585 -1.44813788 0.20826306 0.96957181 1.16398884
55 56 57 58 59 60
1.26889664 1.22462242 0.85214773 3.25445883 4.70200492 5.58265930
61 62 63 64 65 66
-2.87154828 -2.61745844 -1.91582250 -0.25942156 -0.30369578 1.56156000
67 68 69 70 71 72
-1.04864517 1.50284796 2.83104843 1.86184561 3.19004608 3.62315436
73 74 75 76 77 78
-4.62219837 -2.87628821 -3.18909009 -3.26320072 0.15354820 0.13911038
79 80 81 82 83 84
-2.05338511 -0.63567542 1.28925318 4.10927397 4.81090991 3.87154349
85 86 87 88 89 90
-5.11971940 -3.37380924 -3.53742909 -0.65773548 -0.76168252 0.44717233
91 92 93 94 95 96
0.55208013 -1.19286926 2.55208013 2.98518841 3.23927825 4.20944185
97 98 99 100 101 102
-4.69231182 -3.60280260 -1.43918275 0.81394631 1.87361913 1.50114443
103 104 105 106 107 108
0.47226880 -1.66055388 2.20374114 -0.46613685 1.57813738 4.39815817
109 110 111 112 113 114
-8.56134677 -6.09936286 -2.83179597 -1.74228675 1.34626170 1.92855201
115 116 117 118
0.48196668 1.60131231 1.49736527 2.49736527
> postscript(file="/var/wessaorg/rcomp/tmp/65j4e1321954189.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 -6.50040646 NA
1 -5.44073365 -6.50040646
2 -4.76893412 -5.44073365
3 -3.47057005 -4.76893412
4 -2.56007927 -3.47057005
5 -1.41089724 -2.56007927
6 -0.58991568 -1.41089724
7 0.67861198 -0.58991568
8 1.55926635 0.67861198
9 2.67861198 1.55926635
10 3.85763042 2.67861198
11 4.26090229 3.85763042
12 -5.47627076 4.26090229
13 -4.74479842 -5.47627076
14 -3.86414404 -4.74479842
15 -1.93921544 -3.86414404
16 -1.14807029 -1.93921544
17 -0.38676154 -1.14807029
18 0.43422002 -0.38676154
19 1.22536518 0.43422002
20 2.31487440 1.22536518
21 3.46405643 2.31487440
22 3.92700111 3.46405643
23 4.33027299 3.92700111
24 -4.89968115 4.33027299
25 -3.54164428 -4.89968115
26 -2.33278943 -3.54164428
27 -2.19804522 -2.33278943
28 -1.49640928 -2.19804522
29 0.21966447 -1.49640928
30 0.08588103 0.21966447
31 0.16095243 0.08588103
32 1.72784415 0.16095243
33 3.41408150 1.72784415
34 3.54882571 3.41408150
35 3.86258837 3.54882571
36 -6.11327593 3.86258837
37 -4.09883812 -6.11327593
38 -3.57622062 -4.09883812
39 -1.65129202 -3.57622062
40 -0.80047405 -1.65129202
41 1.30347299 -0.80047405
42 1.02050751 1.30347299
43 1.63263423 1.02050751
44 3.31887157 1.63263423
45 4.85592689 3.31887157
46 5.70674486 4.85592689
47 7.27363658 5.70674486
48 -3.06026460 7.27363658
49 -2.29895585 -3.06026460
50 -1.44813788 -2.29895585
51 0.20826306 -1.44813788
52 0.96957181 0.20826306
53 1.16398884 0.96957181
54 1.26889664 1.16398884
55 1.22462242 1.26889664
56 0.85214773 1.22462242
57 3.25445883 0.85214773
58 4.70200492 3.25445883
59 5.58265930 4.70200492
60 -2.87154828 5.58265930
61 -2.61745844 -2.87154828
62 -1.91582250 -2.61745844
63 -0.25942156 -1.91582250
64 -0.30369578 -0.25942156
65 1.56156000 -0.30369578
66 -1.04864517 1.56156000
67 1.50284796 -1.04864517
68 2.83104843 1.50284796
69 1.86184561 2.83104843
70 3.19004608 1.86184561
71 3.62315436 3.19004608
72 -4.62219837 3.62315436
73 -2.87628821 -4.62219837
74 -3.18909009 -2.87628821
75 -3.26320072 -3.18909009
76 0.15354820 -3.26320072
77 0.13911038 0.15354820
78 -2.05338511 0.13911038
79 -0.63567542 -2.05338511
80 1.28925318 -0.63567542
81 4.10927397 1.28925318
82 4.81090991 4.10927397
83 3.87154349 4.81090991
84 -5.11971940 3.87154349
85 -3.37380924 -5.11971940
86 -3.53742909 -3.37380924
87 -0.65773548 -3.53742909
88 -0.76168252 -0.65773548
89 0.44717233 -0.76168252
90 0.55208013 0.44717233
91 -1.19286926 0.55208013
92 2.55208013 -1.19286926
93 2.98518841 2.55208013
94 3.23927825 2.98518841
95 4.20944185 3.23927825
96 -4.69231182 4.20944185
97 -3.60280260 -4.69231182
98 -1.43918275 -3.60280260
99 0.81394631 -1.43918275
100 1.87361913 0.81394631
101 1.50114443 1.87361913
102 0.47226880 1.50114443
103 -1.66055388 0.47226880
104 2.20374114 -1.66055388
105 -0.46613685 2.20374114
106 1.57813738 -0.46613685
107 4.39815817 1.57813738
108 -8.56134677 4.39815817
109 -6.09936286 -8.56134677
110 -2.83179597 -6.09936286
111 -1.74228675 -2.83179597
112 1.34626170 -1.74228675
113 1.92855201 1.34626170
114 0.48196668 1.92855201
115 1.60131231 0.48196668
116 1.49736527 1.60131231
117 2.49736527 1.49736527
118 NA 2.49736527
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -5.44073365 -6.50040646
[2,] -4.76893412 -5.44073365
[3,] -3.47057005 -4.76893412
[4,] -2.56007927 -3.47057005
[5,] -1.41089724 -2.56007927
[6,] -0.58991568 -1.41089724
[7,] 0.67861198 -0.58991568
[8,] 1.55926635 0.67861198
[9,] 2.67861198 1.55926635
[10,] 3.85763042 2.67861198
[11,] 4.26090229 3.85763042
[12,] -5.47627076 4.26090229
[13,] -4.74479842 -5.47627076
[14,] -3.86414404 -4.74479842
[15,] -1.93921544 -3.86414404
[16,] -1.14807029 -1.93921544
[17,] -0.38676154 -1.14807029
[18,] 0.43422002 -0.38676154
[19,] 1.22536518 0.43422002
[20,] 2.31487440 1.22536518
[21,] 3.46405643 2.31487440
[22,] 3.92700111 3.46405643
[23,] 4.33027299 3.92700111
[24,] -4.89968115 4.33027299
[25,] -3.54164428 -4.89968115
[26,] -2.33278943 -3.54164428
[27,] -2.19804522 -2.33278943
[28,] -1.49640928 -2.19804522
[29,] 0.21966447 -1.49640928
[30,] 0.08588103 0.21966447
[31,] 0.16095243 0.08588103
[32,] 1.72784415 0.16095243
[33,] 3.41408150 1.72784415
[34,] 3.54882571 3.41408150
[35,] 3.86258837 3.54882571
[36,] -6.11327593 3.86258837
[37,] -4.09883812 -6.11327593
[38,] -3.57622062 -4.09883812
[39,] -1.65129202 -3.57622062
[40,] -0.80047405 -1.65129202
[41,] 1.30347299 -0.80047405
[42,] 1.02050751 1.30347299
[43,] 1.63263423 1.02050751
[44,] 3.31887157 1.63263423
[45,] 4.85592689 3.31887157
[46,] 5.70674486 4.85592689
[47,] 7.27363658 5.70674486
[48,] -3.06026460 7.27363658
[49,] -2.29895585 -3.06026460
[50,] -1.44813788 -2.29895585
[51,] 0.20826306 -1.44813788
[52,] 0.96957181 0.20826306
[53,] 1.16398884 0.96957181
[54,] 1.26889664 1.16398884
[55,] 1.22462242 1.26889664
[56,] 0.85214773 1.22462242
[57,] 3.25445883 0.85214773
[58,] 4.70200492 3.25445883
[59,] 5.58265930 4.70200492
[60,] -2.87154828 5.58265930
[61,] -2.61745844 -2.87154828
[62,] -1.91582250 -2.61745844
[63,] -0.25942156 -1.91582250
[64,] -0.30369578 -0.25942156
[65,] 1.56156000 -0.30369578
[66,] -1.04864517 1.56156000
[67,] 1.50284796 -1.04864517
[68,] 2.83104843 1.50284796
[69,] 1.86184561 2.83104843
[70,] 3.19004608 1.86184561
[71,] 3.62315436 3.19004608
[72,] -4.62219837 3.62315436
[73,] -2.87628821 -4.62219837
[74,] -3.18909009 -2.87628821
[75,] -3.26320072 -3.18909009
[76,] 0.15354820 -3.26320072
[77,] 0.13911038 0.15354820
[78,] -2.05338511 0.13911038
[79,] -0.63567542 -2.05338511
[80,] 1.28925318 -0.63567542
[81,] 4.10927397 1.28925318
[82,] 4.81090991 4.10927397
[83,] 3.87154349 4.81090991
[84,] -5.11971940 3.87154349
[85,] -3.37380924 -5.11971940
[86,] -3.53742909 -3.37380924
[87,] -0.65773548 -3.53742909
[88,] -0.76168252 -0.65773548
[89,] 0.44717233 -0.76168252
[90,] 0.55208013 0.44717233
[91,] -1.19286926 0.55208013
[92,] 2.55208013 -1.19286926
[93,] 2.98518841 2.55208013
[94,] 3.23927825 2.98518841
[95,] 4.20944185 3.23927825
[96,] -4.69231182 4.20944185
[97,] -3.60280260 -4.69231182
[98,] -1.43918275 -3.60280260
[99,] 0.81394631 -1.43918275
[100,] 1.87361913 0.81394631
[101,] 1.50114443 1.87361913
[102,] 0.47226880 1.50114443
[103,] -1.66055388 0.47226880
[104,] 2.20374114 -1.66055388
[105,] -0.46613685 2.20374114
[106,] 1.57813738 -0.46613685
[107,] 4.39815817 1.57813738
[108,] -8.56134677 4.39815817
[109,] -6.09936286 -8.56134677
[110,] -2.83179597 -6.09936286
[111,] -1.74228675 -2.83179597
[112,] 1.34626170 -1.74228675
[113,] 1.92855201 1.34626170
[114,] 0.48196668 1.92855201
[115,] 1.60131231 0.48196668
[116,] 1.49736527 1.60131231
[117,] 2.49736527 1.49736527
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -5.44073365 -6.50040646
2 -4.76893412 -5.44073365
3 -3.47057005 -4.76893412
4 -2.56007927 -3.47057005
5 -1.41089724 -2.56007927
6 -0.58991568 -1.41089724
7 0.67861198 -0.58991568
8 1.55926635 0.67861198
9 2.67861198 1.55926635
10 3.85763042 2.67861198
11 4.26090229 3.85763042
12 -5.47627076 4.26090229
13 -4.74479842 -5.47627076
14 -3.86414404 -4.74479842
15 -1.93921544 -3.86414404
16 -1.14807029 -1.93921544
17 -0.38676154 -1.14807029
18 0.43422002 -0.38676154
19 1.22536518 0.43422002
20 2.31487440 1.22536518
21 3.46405643 2.31487440
22 3.92700111 3.46405643
23 4.33027299 3.92700111
24 -4.89968115 4.33027299
25 -3.54164428 -4.89968115
26 -2.33278943 -3.54164428
27 -2.19804522 -2.33278943
28 -1.49640928 -2.19804522
29 0.21966447 -1.49640928
30 0.08588103 0.21966447
31 0.16095243 0.08588103
32 1.72784415 0.16095243
33 3.41408150 1.72784415
34 3.54882571 3.41408150
35 3.86258837 3.54882571
36 -6.11327593 3.86258837
37 -4.09883812 -6.11327593
38 -3.57622062 -4.09883812
39 -1.65129202 -3.57622062
40 -0.80047405 -1.65129202
41 1.30347299 -0.80047405
42 1.02050751 1.30347299
43 1.63263423 1.02050751
44 3.31887157 1.63263423
45 4.85592689 3.31887157
46 5.70674486 4.85592689
47 7.27363658 5.70674486
48 -3.06026460 7.27363658
49 -2.29895585 -3.06026460
50 -1.44813788 -2.29895585
51 0.20826306 -1.44813788
52 0.96957181 0.20826306
53 1.16398884 0.96957181
54 1.26889664 1.16398884
55 1.22462242 1.26889664
56 0.85214773 1.22462242
57 3.25445883 0.85214773
58 4.70200492 3.25445883
59 5.58265930 4.70200492
60 -2.87154828 5.58265930
61 -2.61745844 -2.87154828
62 -1.91582250 -2.61745844
63 -0.25942156 -1.91582250
64 -0.30369578 -0.25942156
65 1.56156000 -0.30369578
66 -1.04864517 1.56156000
67 1.50284796 -1.04864517
68 2.83104843 1.50284796
69 1.86184561 2.83104843
70 3.19004608 1.86184561
71 3.62315436 3.19004608
72 -4.62219837 3.62315436
73 -2.87628821 -4.62219837
74 -3.18909009 -2.87628821
75 -3.26320072 -3.18909009
76 0.15354820 -3.26320072
77 0.13911038 0.15354820
78 -2.05338511 0.13911038
79 -0.63567542 -2.05338511
80 1.28925318 -0.63567542
81 4.10927397 1.28925318
82 4.81090991 4.10927397
83 3.87154349 4.81090991
84 -5.11971940 3.87154349
85 -3.37380924 -5.11971940
86 -3.53742909 -3.37380924
87 -0.65773548 -3.53742909
88 -0.76168252 -0.65773548
89 0.44717233 -0.76168252
90 0.55208013 0.44717233
91 -1.19286926 0.55208013
92 2.55208013 -1.19286926
93 2.98518841 2.55208013
94 3.23927825 2.98518841
95 4.20944185 3.23927825
96 -4.69231182 4.20944185
97 -3.60280260 -4.69231182
98 -1.43918275 -3.60280260
99 0.81394631 -1.43918275
100 1.87361913 0.81394631
101 1.50114443 1.87361913
102 0.47226880 1.50114443
103 -1.66055388 0.47226880
104 2.20374114 -1.66055388
105 -0.46613685 2.20374114
106 1.57813738 -0.46613685
107 4.39815817 1.57813738
108 -8.56134677 4.39815817
109 -6.09936286 -8.56134677
110 -2.83179597 -6.09936286
111 -1.74228675 -2.83179597
112 1.34626170 -1.74228675
113 1.92855201 1.34626170
114 0.48196668 1.92855201
115 1.60131231 0.48196668
116 1.49736527 1.60131231
117 2.49736527 1.49736527
> 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/7lae11321954189.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/8ys9n1321954189.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/9p0fy1321954189.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/10khwj1321954189.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/11470n1321954189.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/12jito1321954189.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/136p2c1321954189.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/1486qu1321954189.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/15g5we1321954189.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/16x51m1321954189.tab")
+ }
>
> try(system("convert tmp/1owrs1321954189.ps tmp/1owrs1321954189.png",intern=TRUE))
character(0)
> try(system("convert tmp/29k1y1321954189.ps tmp/29k1y1321954189.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ymbh1321954189.ps tmp/3ymbh1321954189.png",intern=TRUE))
character(0)
> try(system("convert tmp/4jbh51321954189.ps tmp/4jbh51321954189.png",intern=TRUE))
character(0)
> try(system("convert tmp/55ovd1321954189.ps tmp/55ovd1321954189.png",intern=TRUE))
character(0)
> try(system("convert tmp/65j4e1321954189.ps tmp/65j4e1321954189.png",intern=TRUE))
character(0)
> try(system("convert tmp/7lae11321954189.ps tmp/7lae11321954189.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ys9n1321954189.ps tmp/8ys9n1321954189.png",intern=TRUE))
character(0)
> try(system("convert tmp/9p0fy1321954189.ps tmp/9p0fy1321954189.png",intern=TRUE))
character(0)
> try(system("convert tmp/10khwj1321954189.ps tmp/10khwj1321954189.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.995 0.634 4.762