R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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('month','robberies'),1:118))
> y <- array(NA,dim=c(2,118),dimnames=list(c('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'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal Dummies'
> par1 <- '2'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, as.Date.numeric
> 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 month
1 41 1
2 39 2
3 50 3
4 40 4
5 43 5
6 38 6
7 44 7
8 35 8
9 39 9
10 35 10
11 29 11
12 49 12
13 50 1
14 59 2
15 63 3
16 32 4
17 39 5
18 47 6
19 53 7
20 60 8
21 57 9
22 52 10
23 70 11
24 90 12
25 74 1
26 62 2
27 55 3
28 84 4
29 94 5
30 70 6
31 108 7
32 139 8
33 120 9
34 97 10
35 126 11
36 149 12
37 158 1
38 124 2
39 140 3
40 109 4
41 114 5
42 77 6
43 120 7
44 133 8
45 110 9
46 92 10
47 97 11
48 78 12
49 99 1
50 107 2
51 112 3
52 90 4
53 98 5
54 125 6
55 155 7
56 190 8
57 236 9
58 189 10
59 174 11
60 178 12
61 136 1
62 161 2
63 171 3
64 149 4
65 184 5
66 155 6
67 276 7
68 224 8
69 213 9
70 279 10
71 268 11
72 287 12
73 238 1
74 213 2
75 257 3
76 293 4
77 212 5
78 246 6
79 353 7
80 339 8
81 308 9
82 247 10
83 257 11
84 322 12
85 298 1
86 273 2
87 312 3
88 249 4
89 286 5
90 279 6
91 309 7
92 401 8
93 309 9
94 328 10
95 353 11
96 354 12
97 327 1
98 324 2
99 285 3
100 243 4
101 241 5
102 287 6
103 355 7
104 460 8
105 364 9
106 487 10
107 452 11
108 391 12
109 500 1
110 451 2
111 375 3
112 372 4
113 302 5
114 316 6
115 398 7
116 394 8
117 431 9
118 431 10
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) month
168.006 4.409
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-187.50 -116.87 -29.75 101.26 327.59
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 168.006 24.999 6.721 7.12e-10 ***
month 4.409 3.439 1.282 0.202
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 127.7 on 116 degrees of freedom
Multiple R-squared: 0.01397, Adjusted R-squared: 0.005468
F-statistic: 1.643 on 1 and 116 DF, p-value: 0.2024
> 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,] 8.455396e-05 1.691079e-04 0.999915446
[2,] 4.241213e-06 8.482425e-06 0.999995759
[3,] 1.673869e-07 3.347737e-07 0.999999833
[4,] 1.192142e-08 2.384285e-08 0.999999988
[5,] 4.173297e-10 8.346594e-10 1.000000000
[6,] 1.710538e-11 3.421077e-11 1.000000000
[7,] 1.693836e-12 3.387673e-12 1.000000000
[8,] 2.278228e-12 4.556456e-12 1.000000000
[9,] 2.014821e-13 4.029641e-13 1.000000000
[10,] 1.149839e-13 2.299679e-13 1.000000000
[11,] 8.979604e-14 1.795921e-13 1.000000000
[12,] 2.254970e-14 4.509940e-14 1.000000000
[13,] 2.039996e-15 4.079993e-15 1.000000000
[14,] 1.903107e-16 3.806215e-16 1.000000000
[15,] 3.841006e-17 7.682011e-17 1.000000000
[16,] 2.722247e-17 5.444493e-17 1.000000000
[17,] 8.620147e-18 1.724029e-17 1.000000000
[18,] 1.387939e-18 2.775878e-18 1.000000000
[19,] 2.957636e-18 5.915272e-18 1.000000000
[20,] 7.513041e-17 1.502608e-16 1.000000000
[21,] 8.264036e-17 1.652807e-16 1.000000000
[22,] 2.015267e-17 4.030534e-17 1.000000000
[23,] 3.637077e-18 7.274154e-18 1.000000000
[24,] 7.164778e-18 1.432956e-17 1.000000000
[25,] 3.202968e-17 6.405936e-17 1.000000000
[26,] 1.119415e-17 2.238831e-17 1.000000000
[27,] 1.442825e-16 2.885649e-16 1.000000000
[28,] 2.404437e-14 4.808874e-14 1.000000000
[29,] 8.790752e-14 1.758150e-13 1.000000000
[30,] 6.182472e-14 1.236494e-13 1.000000000
[31,] 1.501105e-13 3.002210e-13 1.000000000
[32,] 7.967917e-13 1.593583e-12 1.000000000
[33,] 2.561965e-11 5.123930e-11 1.000000000
[34,] 3.767670e-11 7.535340e-11 1.000000000
[35,] 8.155118e-11 1.631024e-10 1.000000000
[36,] 6.039848e-11 1.207970e-10 1.000000000
[37,] 4.969885e-11 9.939769e-11 1.000000000
[38,] 3.185196e-11 6.370392e-11 1.000000000
[39,] 3.181359e-11 6.362718e-11 1.000000000
[40,] 4.328868e-11 8.657737e-11 1.000000000
[41,] 3.792991e-11 7.585982e-11 1.000000000
[42,] 3.382254e-11 6.764507e-11 1.000000000
[43,] 3.562191e-11 7.124382e-11 1.000000000
[44,] 5.920640e-11 1.184128e-10 1.000000000
[45,] 4.782219e-11 9.564437e-11 1.000000000
[46,] 4.614844e-11 9.229687e-11 1.000000000
[47,] 5.230820e-11 1.046164e-10 1.000000000
[48,] 6.495004e-11 1.299001e-10 1.000000000
[49,] 9.617750e-11 1.923550e-10 1.000000000
[50,] 1.988092e-10 3.976185e-10 1.000000000
[51,] 8.285993e-10 1.657199e-09 0.999999999
[52,] 9.405682e-09 1.881136e-08 0.999999991
[53,] 3.526395e-07 7.052791e-07 0.999999647
[54,] 1.355373e-06 2.710746e-06 0.999998645
[55,] 3.897734e-06 7.795468e-06 0.999996102
[56,] 1.232895e-05 2.465791e-05 0.999987671
[57,] 2.170990e-05 4.341979e-05 0.999978290
[58,] 4.628654e-05 9.257308e-05 0.999953713
[59,] 1.060572e-04 2.121144e-04 0.999893943
[60,] 2.451735e-04 4.903469e-04 0.999754827
[61,] 6.288725e-04 1.257745e-03 0.999371127
[62,] 1.701785e-03 3.403571e-03 0.998298215
[63,] 9.138028e-03 1.827606e-02 0.990861972
[64,] 1.982673e-02 3.965346e-02 0.980173271
[65,] 3.905790e-02 7.811580e-02 0.960942101
[66,] 8.367858e-02 1.673572e-01 0.916321417
[67,] 1.388934e-01 2.777867e-01 0.861106637
[68,] 2.113521e-01 4.227043e-01 0.788647873
[69,] 2.758234e-01 5.516468e-01 0.724176606
[70,] 3.395804e-01 6.791607e-01 0.660419647
[71,] 4.096302e-01 8.192605e-01 0.590369757
[72,] 4.928319e-01 9.856639e-01 0.507168070
[73,] 5.771402e-01 8.457196e-01 0.422859791
[74,] 6.420470e-01 7.159059e-01 0.357952962
[75,] 7.403841e-01 5.192318e-01 0.259615920
[76,] 7.921803e-01 4.156393e-01 0.207819662
[77,] 8.158007e-01 3.683986e-01 0.184199298
[78,] 8.574014e-01 2.851971e-01 0.142598565
[79,] 8.968227e-01 2.063545e-01 0.103177265
[80,] 9.119241e-01 1.761517e-01 0.088075871
[81,] 9.175249e-01 1.649503e-01 0.082475132
[82,] 9.199071e-01 1.601858e-01 0.080092887
[83,] 9.211406e-01 1.577188e-01 0.078859392
[84,] 9.334597e-01 1.330806e-01 0.066540300
[85,] 9.367879e-01 1.264242e-01 0.063212125
[86,] 9.436177e-01 1.127646e-01 0.056382286
[87,] 9.446596e-01 1.106808e-01 0.055340376
[88,] 9.513566e-01 9.728671e-02 0.048643356
[89,] 9.520818e-01 9.583641e-02 0.047918203
[90,] 9.505385e-01 9.892290e-02 0.049461450
[91,] 9.461902e-01 1.076196e-01 0.053809777
[92,] 9.430763e-01 1.138473e-01 0.056923667
[93,] 9.320392e-01 1.359217e-01 0.067960846
[94,] 9.172682e-01 1.654636e-01 0.082731819
[95,] 9.110637e-01 1.778726e-01 0.088936319
[96,] 9.421158e-01 1.157683e-01 0.057884153
[97,] 9.773945e-01 4.521096e-02 0.022605480
[98,] 9.877428e-01 2.451430e-02 0.012257150
[99,] 9.850499e-01 2.990025e-02 0.014950126
[100,] 9.845164e-01 3.096711e-02 0.015483556
[101,] 9.784069e-01 4.318611e-02 0.021593057
[102,] 9.834544e-01 3.309110e-02 0.016545551
[103,] 9.808019e-01 3.839624e-02 0.019198120
[104,] 9.643181e-01 7.136376e-02 0.035681882
[105,] 9.847105e-01 3.057901e-02 0.015289504
[106,] 9.955608e-01 8.878312e-03 0.004439156
[107,] 9.941874e-01 1.162517e-02 0.005812584
[108,] 9.979035e-01 4.192902e-03 0.002096451
[109,] 9.884789e-01 2.304222e-02 0.011521108
> postscript(file="/var/wessaorg/rcomp/tmp/1g98c1353426653.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/21kya1353426653.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/31u3u1353426653.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/4r1ts1353426653.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/55i901353426653.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
-131.414627 -137.823194 -131.231760 -145.640327 -147.048893 -156.457460
7 8 9 10 11 12
-154.866026 -168.274593 -168.683159 -177.091726 -187.500292 -171.908858
13 14 15 16 17 18
-122.414627 -117.823194 -118.231760 -153.640327 -151.048893 -147.457460
19 20 21 22 23 24
-145.866026 -143.274593 -150.683159 -160.091726 -146.500292 -130.908858
25 26 27 28 29 30
-98.414627 -114.823194 -126.231760 -101.640327 -96.048893 -124.457460
31 32 33 34 35 36
-90.866026 -64.274593 -87.683159 -115.091726 -90.500292 -71.908858
37 38 39 40 41 42
-14.414627 -52.823194 -41.231760 -76.640327 -76.048893 -117.457460
43 44 45 46 47 48
-78.866026 -70.274593 -97.683159 -120.091726 -119.500292 -142.908858
49 50 51 52 53 54
-73.414627 -69.823194 -69.231760 -95.640327 -92.048893 -69.457460
55 56 57 58 59 60
-43.866026 -13.274593 28.316841 -23.091726 -42.500292 -42.908858
61 62 63 64 65 66
-36.414627 -15.823194 -10.231760 -36.640327 -6.048893 -39.457460
67 68 69 70 71 72
77.133974 20.725407 5.316841 66.908274 51.499708 66.091142
73 74 75 76 77 78
65.585373 36.176806 75.768240 107.359673 21.951107 51.542540
79 80 81 82 83 84
154.133974 135.725407 100.316841 34.908274 40.499708 101.091142
85 86 87 88 89 90
125.585373 96.176806 130.768240 63.359673 95.951107 84.542540
91 92 93 94 95 96
110.133974 197.725407 101.316841 115.908274 136.499708 133.091142
97 98 99 100 101 102
154.585373 147.176806 103.768240 57.359673 50.951107 92.542540
103 104 105 106 107 108
156.133974 256.725407 156.316841 274.908274 235.499708 170.091142
109 110 111 112 113 114
327.585373 274.176806 193.768240 186.359673 111.951107 121.542540
115 116 117 118
199.133974 190.725407 223.316841 218.908274
> postscript(file="/var/wessaorg/rcomp/tmp/6vtbz1353426653.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 -131.414627 NA
1 -137.823194 -131.414627
2 -131.231760 -137.823194
3 -145.640327 -131.231760
4 -147.048893 -145.640327
5 -156.457460 -147.048893
6 -154.866026 -156.457460
7 -168.274593 -154.866026
8 -168.683159 -168.274593
9 -177.091726 -168.683159
10 -187.500292 -177.091726
11 -171.908858 -187.500292
12 -122.414627 -171.908858
13 -117.823194 -122.414627
14 -118.231760 -117.823194
15 -153.640327 -118.231760
16 -151.048893 -153.640327
17 -147.457460 -151.048893
18 -145.866026 -147.457460
19 -143.274593 -145.866026
20 -150.683159 -143.274593
21 -160.091726 -150.683159
22 -146.500292 -160.091726
23 -130.908858 -146.500292
24 -98.414627 -130.908858
25 -114.823194 -98.414627
26 -126.231760 -114.823194
27 -101.640327 -126.231760
28 -96.048893 -101.640327
29 -124.457460 -96.048893
30 -90.866026 -124.457460
31 -64.274593 -90.866026
32 -87.683159 -64.274593
33 -115.091726 -87.683159
34 -90.500292 -115.091726
35 -71.908858 -90.500292
36 -14.414627 -71.908858
37 -52.823194 -14.414627
38 -41.231760 -52.823194
39 -76.640327 -41.231760
40 -76.048893 -76.640327
41 -117.457460 -76.048893
42 -78.866026 -117.457460
43 -70.274593 -78.866026
44 -97.683159 -70.274593
45 -120.091726 -97.683159
46 -119.500292 -120.091726
47 -142.908858 -119.500292
48 -73.414627 -142.908858
49 -69.823194 -73.414627
50 -69.231760 -69.823194
51 -95.640327 -69.231760
52 -92.048893 -95.640327
53 -69.457460 -92.048893
54 -43.866026 -69.457460
55 -13.274593 -43.866026
56 28.316841 -13.274593
57 -23.091726 28.316841
58 -42.500292 -23.091726
59 -42.908858 -42.500292
60 -36.414627 -42.908858
61 -15.823194 -36.414627
62 -10.231760 -15.823194
63 -36.640327 -10.231760
64 -6.048893 -36.640327
65 -39.457460 -6.048893
66 77.133974 -39.457460
67 20.725407 77.133974
68 5.316841 20.725407
69 66.908274 5.316841
70 51.499708 66.908274
71 66.091142 51.499708
72 65.585373 66.091142
73 36.176806 65.585373
74 75.768240 36.176806
75 107.359673 75.768240
76 21.951107 107.359673
77 51.542540 21.951107
78 154.133974 51.542540
79 135.725407 154.133974
80 100.316841 135.725407
81 34.908274 100.316841
82 40.499708 34.908274
83 101.091142 40.499708
84 125.585373 101.091142
85 96.176806 125.585373
86 130.768240 96.176806
87 63.359673 130.768240
88 95.951107 63.359673
89 84.542540 95.951107
90 110.133974 84.542540
91 197.725407 110.133974
92 101.316841 197.725407
93 115.908274 101.316841
94 136.499708 115.908274
95 133.091142 136.499708
96 154.585373 133.091142
97 147.176806 154.585373
98 103.768240 147.176806
99 57.359673 103.768240
100 50.951107 57.359673
101 92.542540 50.951107
102 156.133974 92.542540
103 256.725407 156.133974
104 156.316841 256.725407
105 274.908274 156.316841
106 235.499708 274.908274
107 170.091142 235.499708
108 327.585373 170.091142
109 274.176806 327.585373
110 193.768240 274.176806
111 186.359673 193.768240
112 111.951107 186.359673
113 121.542540 111.951107
114 199.133974 121.542540
115 190.725407 199.133974
116 223.316841 190.725407
117 218.908274 223.316841
118 NA 218.908274
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -137.823194 -131.414627
[2,] -131.231760 -137.823194
[3,] -145.640327 -131.231760
[4,] -147.048893 -145.640327
[5,] -156.457460 -147.048893
[6,] -154.866026 -156.457460
[7,] -168.274593 -154.866026
[8,] -168.683159 -168.274593
[9,] -177.091726 -168.683159
[10,] -187.500292 -177.091726
[11,] -171.908858 -187.500292
[12,] -122.414627 -171.908858
[13,] -117.823194 -122.414627
[14,] -118.231760 -117.823194
[15,] -153.640327 -118.231760
[16,] -151.048893 -153.640327
[17,] -147.457460 -151.048893
[18,] -145.866026 -147.457460
[19,] -143.274593 -145.866026
[20,] -150.683159 -143.274593
[21,] -160.091726 -150.683159
[22,] -146.500292 -160.091726
[23,] -130.908858 -146.500292
[24,] -98.414627 -130.908858
[25,] -114.823194 -98.414627
[26,] -126.231760 -114.823194
[27,] -101.640327 -126.231760
[28,] -96.048893 -101.640327
[29,] -124.457460 -96.048893
[30,] -90.866026 -124.457460
[31,] -64.274593 -90.866026
[32,] -87.683159 -64.274593
[33,] -115.091726 -87.683159
[34,] -90.500292 -115.091726
[35,] -71.908858 -90.500292
[36,] -14.414627 -71.908858
[37,] -52.823194 -14.414627
[38,] -41.231760 -52.823194
[39,] -76.640327 -41.231760
[40,] -76.048893 -76.640327
[41,] -117.457460 -76.048893
[42,] -78.866026 -117.457460
[43,] -70.274593 -78.866026
[44,] -97.683159 -70.274593
[45,] -120.091726 -97.683159
[46,] -119.500292 -120.091726
[47,] -142.908858 -119.500292
[48,] -73.414627 -142.908858
[49,] -69.823194 -73.414627
[50,] -69.231760 -69.823194
[51,] -95.640327 -69.231760
[52,] -92.048893 -95.640327
[53,] -69.457460 -92.048893
[54,] -43.866026 -69.457460
[55,] -13.274593 -43.866026
[56,] 28.316841 -13.274593
[57,] -23.091726 28.316841
[58,] -42.500292 -23.091726
[59,] -42.908858 -42.500292
[60,] -36.414627 -42.908858
[61,] -15.823194 -36.414627
[62,] -10.231760 -15.823194
[63,] -36.640327 -10.231760
[64,] -6.048893 -36.640327
[65,] -39.457460 -6.048893
[66,] 77.133974 -39.457460
[67,] 20.725407 77.133974
[68,] 5.316841 20.725407
[69,] 66.908274 5.316841
[70,] 51.499708 66.908274
[71,] 66.091142 51.499708
[72,] 65.585373 66.091142
[73,] 36.176806 65.585373
[74,] 75.768240 36.176806
[75,] 107.359673 75.768240
[76,] 21.951107 107.359673
[77,] 51.542540 21.951107
[78,] 154.133974 51.542540
[79,] 135.725407 154.133974
[80,] 100.316841 135.725407
[81,] 34.908274 100.316841
[82,] 40.499708 34.908274
[83,] 101.091142 40.499708
[84,] 125.585373 101.091142
[85,] 96.176806 125.585373
[86,] 130.768240 96.176806
[87,] 63.359673 130.768240
[88,] 95.951107 63.359673
[89,] 84.542540 95.951107
[90,] 110.133974 84.542540
[91,] 197.725407 110.133974
[92,] 101.316841 197.725407
[93,] 115.908274 101.316841
[94,] 136.499708 115.908274
[95,] 133.091142 136.499708
[96,] 154.585373 133.091142
[97,] 147.176806 154.585373
[98,] 103.768240 147.176806
[99,] 57.359673 103.768240
[100,] 50.951107 57.359673
[101,] 92.542540 50.951107
[102,] 156.133974 92.542540
[103,] 256.725407 156.133974
[104,] 156.316841 256.725407
[105,] 274.908274 156.316841
[106,] 235.499708 274.908274
[107,] 170.091142 235.499708
[108,] 327.585373 170.091142
[109,] 274.176806 327.585373
[110,] 193.768240 274.176806
[111,] 186.359673 193.768240
[112,] 111.951107 186.359673
[113,] 121.542540 111.951107
[114,] 199.133974 121.542540
[115,] 190.725407 199.133974
[116,] 223.316841 190.725407
[117,] 218.908274 223.316841
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -137.823194 -131.414627
2 -131.231760 -137.823194
3 -145.640327 -131.231760
4 -147.048893 -145.640327
5 -156.457460 -147.048893
6 -154.866026 -156.457460
7 -168.274593 -154.866026
8 -168.683159 -168.274593
9 -177.091726 -168.683159
10 -187.500292 -177.091726
11 -171.908858 -187.500292
12 -122.414627 -171.908858
13 -117.823194 -122.414627
14 -118.231760 -117.823194
15 -153.640327 -118.231760
16 -151.048893 -153.640327
17 -147.457460 -151.048893
18 -145.866026 -147.457460
19 -143.274593 -145.866026
20 -150.683159 -143.274593
21 -160.091726 -150.683159
22 -146.500292 -160.091726
23 -130.908858 -146.500292
24 -98.414627 -130.908858
25 -114.823194 -98.414627
26 -126.231760 -114.823194
27 -101.640327 -126.231760
28 -96.048893 -101.640327
29 -124.457460 -96.048893
30 -90.866026 -124.457460
31 -64.274593 -90.866026
32 -87.683159 -64.274593
33 -115.091726 -87.683159
34 -90.500292 -115.091726
35 -71.908858 -90.500292
36 -14.414627 -71.908858
37 -52.823194 -14.414627
38 -41.231760 -52.823194
39 -76.640327 -41.231760
40 -76.048893 -76.640327
41 -117.457460 -76.048893
42 -78.866026 -117.457460
43 -70.274593 -78.866026
44 -97.683159 -70.274593
45 -120.091726 -97.683159
46 -119.500292 -120.091726
47 -142.908858 -119.500292
48 -73.414627 -142.908858
49 -69.823194 -73.414627
50 -69.231760 -69.823194
51 -95.640327 -69.231760
52 -92.048893 -95.640327
53 -69.457460 -92.048893
54 -43.866026 -69.457460
55 -13.274593 -43.866026
56 28.316841 -13.274593
57 -23.091726 28.316841
58 -42.500292 -23.091726
59 -42.908858 -42.500292
60 -36.414627 -42.908858
61 -15.823194 -36.414627
62 -10.231760 -15.823194
63 -36.640327 -10.231760
64 -6.048893 -36.640327
65 -39.457460 -6.048893
66 77.133974 -39.457460
67 20.725407 77.133974
68 5.316841 20.725407
69 66.908274 5.316841
70 51.499708 66.908274
71 66.091142 51.499708
72 65.585373 66.091142
73 36.176806 65.585373
74 75.768240 36.176806
75 107.359673 75.768240
76 21.951107 107.359673
77 51.542540 21.951107
78 154.133974 51.542540
79 135.725407 154.133974
80 100.316841 135.725407
81 34.908274 100.316841
82 40.499708 34.908274
83 101.091142 40.499708
84 125.585373 101.091142
85 96.176806 125.585373
86 130.768240 96.176806
87 63.359673 130.768240
88 95.951107 63.359673
89 84.542540 95.951107
90 110.133974 84.542540
91 197.725407 110.133974
92 101.316841 197.725407
93 115.908274 101.316841
94 136.499708 115.908274
95 133.091142 136.499708
96 154.585373 133.091142
97 147.176806 154.585373
98 103.768240 147.176806
99 57.359673 103.768240
100 50.951107 57.359673
101 92.542540 50.951107
102 156.133974 92.542540
103 256.725407 156.133974
104 156.316841 256.725407
105 274.908274 156.316841
106 235.499708 274.908274
107 170.091142 235.499708
108 327.585373 170.091142
109 274.176806 327.585373
110 193.768240 274.176806
111 186.359673 193.768240
112 111.951107 186.359673
113 121.542540 111.951107
114 199.133974 121.542540
115 190.725407 199.133974
116 223.316841 190.725407
117 218.908274 223.316841
> 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/7k6sz1353426653.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/8ar551353426653.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/91ukn1353426653.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/10w2fk1353426653.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/114lko1353426654.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/12yfou1353426654.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/13svxa1353426654.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/14m3o21353426654.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/152hd91353426654.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/16a2sx1353426654.tab")
+ }
>
> try(system("convert tmp/1g98c1353426653.ps tmp/1g98c1353426653.png",intern=TRUE))
character(0)
> try(system("convert tmp/21kya1353426653.ps tmp/21kya1353426653.png",intern=TRUE))
character(0)
> try(system("convert tmp/31u3u1353426653.ps tmp/31u3u1353426653.png",intern=TRUE))
character(0)
> try(system("convert tmp/4r1ts1353426653.ps tmp/4r1ts1353426653.png",intern=TRUE))
character(0)
> try(system("convert tmp/55i901353426653.ps tmp/55i901353426653.png",intern=TRUE))
character(0)
> try(system("convert tmp/6vtbz1353426653.ps tmp/6vtbz1353426653.png",intern=TRUE))
character(0)
> try(system("convert tmp/7k6sz1353426653.ps tmp/7k6sz1353426653.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ar551353426653.ps tmp/8ar551353426653.png",intern=TRUE))
character(0)
> try(system("convert tmp/91ukn1353426653.ps tmp/91ukn1353426653.png",intern=TRUE))
character(0)
> try(system("convert tmp/10w2fk1353426653.ps tmp/10w2fk1353426653.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
8.754 1.031 9.797