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(41,39,50,40,43,38,44,35,39,35,29,49,50,59,63,32,39,47,53,60,57,52,70,90,74,62,55,84,94,70,108,139,120,97,126,149,158,124,140,109,114,77,120,133,110,92,97,78,99,107,112,90,98,125,155,190,236,189,174,178,136,161,171,149,184,155,276,224,213,279,268,287,238,213,257,293,212,246,353,339,308,247,257,322,298,273,312,249,286,279,309,401,309,328,353,354,327,324,285,243,241,287,355,460,364,487,452,391,500,451,375,372,302,316,398,394,431,431),dim=c(1,118),dimnames=list(c('robberies'),1:118))
> y <- array(NA,dim=c(1,118),dimnames=list(c('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 = 'Include Monthly Dummies'
> par1 = '1'
> par3 <- 'No Linear Trend'
> par2 <- 'Include Monthly Dummies'
> par1 <- '1'
> #'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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 41 1 0 0 0 0 0 0 0 0 0 0
2 39 0 1 0 0 0 0 0 0 0 0 0
3 50 0 0 1 0 0 0 0 0 0 0 0
4 40 0 0 0 1 0 0 0 0 0 0 0
5 43 0 0 0 0 1 0 0 0 0 0 0
6 38 0 0 0 0 0 1 0 0 0 0 0
7 44 0 0 0 0 0 0 1 0 0 0 0
8 35 0 0 0 0 0 0 0 1 0 0 0
9 39 0 0 0 0 0 0 0 0 1 0 0
10 35 0 0 0 0 0 0 0 0 0 1 0
11 29 0 0 0 0 0 0 0 0 0 0 1
12 49 0 0 0 0 0 0 0 0 0 0 0
13 50 1 0 0 0 0 0 0 0 0 0 0
14 59 0 1 0 0 0 0 0 0 0 0 0
15 63 0 0 1 0 0 0 0 0 0 0 0
16 32 0 0 0 1 0 0 0 0 0 0 0
17 39 0 0 0 0 1 0 0 0 0 0 0
18 47 0 0 0 0 0 1 0 0 0 0 0
19 53 0 0 0 0 0 0 1 0 0 0 0
20 60 0 0 0 0 0 0 0 1 0 0 0
21 57 0 0 0 0 0 0 0 0 1 0 0
22 52 0 0 0 0 0 0 0 0 0 1 0
23 70 0 0 0 0 0 0 0 0 0 0 1
24 90 0 0 0 0 0 0 0 0 0 0 0
25 74 1 0 0 0 0 0 0 0 0 0 0
26 62 0 1 0 0 0 0 0 0 0 0 0
27 55 0 0 1 0 0 0 0 0 0 0 0
28 84 0 0 0 1 0 0 0 0 0 0 0
29 94 0 0 0 0 1 0 0 0 0 0 0
30 70 0 0 0 0 0 1 0 0 0 0 0
31 108 0 0 0 0 0 0 1 0 0 0 0
32 139 0 0 0 0 0 0 0 1 0 0 0
33 120 0 0 0 0 0 0 0 0 1 0 0
34 97 0 0 0 0 0 0 0 0 0 1 0
35 126 0 0 0 0 0 0 0 0 0 0 1
36 149 0 0 0 0 0 0 0 0 0 0 0
37 158 1 0 0 0 0 0 0 0 0 0 0
38 124 0 1 0 0 0 0 0 0 0 0 0
39 140 0 0 1 0 0 0 0 0 0 0 0
40 109 0 0 0 1 0 0 0 0 0 0 0
41 114 0 0 0 0 1 0 0 0 0 0 0
42 77 0 0 0 0 0 1 0 0 0 0 0
43 120 0 0 0 0 0 0 1 0 0 0 0
44 133 0 0 0 0 0 0 0 1 0 0 0
45 110 0 0 0 0 0 0 0 0 1 0 0
46 92 0 0 0 0 0 0 0 0 0 1 0
47 97 0 0 0 0 0 0 0 0 0 0 1
48 78 0 0 0 0 0 0 0 0 0 0 0
49 99 1 0 0 0 0 0 0 0 0 0 0
50 107 0 1 0 0 0 0 0 0 0 0 0
51 112 0 0 1 0 0 0 0 0 0 0 0
52 90 0 0 0 1 0 0 0 0 0 0 0
53 98 0 0 0 0 1 0 0 0 0 0 0
54 125 0 0 0 0 0 1 0 0 0 0 0
55 155 0 0 0 0 0 0 1 0 0 0 0
56 190 0 0 0 0 0 0 0 1 0 0 0
57 236 0 0 0 0 0 0 0 0 1 0 0
58 189 0 0 0 0 0 0 0 0 0 1 0
59 174 0 0 0 0 0 0 0 0 0 0 1
60 178 0 0 0 0 0 0 0 0 0 0 0
61 136 1 0 0 0 0 0 0 0 0 0 0
62 161 0 1 0 0 0 0 0 0 0 0 0
63 171 0 0 1 0 0 0 0 0 0 0 0
64 149 0 0 0 1 0 0 0 0 0 0 0
65 184 0 0 0 0 1 0 0 0 0 0 0
66 155 0 0 0 0 0 1 0 0 0 0 0
67 276 0 0 0 0 0 0 1 0 0 0 0
68 224 0 0 0 0 0 0 0 1 0 0 0
69 213 0 0 0 0 0 0 0 0 1 0 0
70 279 0 0 0 0 0 0 0 0 0 1 0
71 268 0 0 0 0 0 0 0 0 0 0 1
72 287 0 0 0 0 0 0 0 0 0 0 0
73 238 1 0 0 0 0 0 0 0 0 0 0
74 213 0 1 0 0 0 0 0 0 0 0 0
75 257 0 0 1 0 0 0 0 0 0 0 0
76 293 0 0 0 1 0 0 0 0 0 0 0
77 212 0 0 0 0 1 0 0 0 0 0 0
78 246 0 0 0 0 0 1 0 0 0 0 0
79 353 0 0 0 0 0 0 1 0 0 0 0
80 339 0 0 0 0 0 0 0 1 0 0 0
81 308 0 0 0 0 0 0 0 0 1 0 0
82 247 0 0 0 0 0 0 0 0 0 1 0
83 257 0 0 0 0 0 0 0 0 0 0 1
84 322 0 0 0 0 0 0 0 0 0 0 0
85 298 1 0 0 0 0 0 0 0 0 0 0
86 273 0 1 0 0 0 0 0 0 0 0 0
87 312 0 0 1 0 0 0 0 0 0 0 0
88 249 0 0 0 1 0 0 0 0 0 0 0
89 286 0 0 0 0 1 0 0 0 0 0 0
90 279 0 0 0 0 0 1 0 0 0 0 0
91 309 0 0 0 0 0 0 1 0 0 0 0
92 401 0 0 0 0 0 0 0 1 0 0 0
93 309 0 0 0 0 0 0 0 0 1 0 0
94 328 0 0 0 0 0 0 0 0 0 1 0
95 353 0 0 0 0 0 0 0 0 0 0 1
96 354 0 0 0 0 0 0 0 0 0 0 0
97 327 1 0 0 0 0 0 0 0 0 0 0
98 324 0 1 0 0 0 0 0 0 0 0 0
99 285 0 0 1 0 0 0 0 0 0 0 0
100 243 0 0 0 1 0 0 0 0 0 0 0
101 241 0 0 0 0 1 0 0 0 0 0 0
102 287 0 0 0 0 0 1 0 0 0 0 0
103 355 0 0 0 0 0 0 1 0 0 0 0
104 460 0 0 0 0 0 0 0 1 0 0 0
105 364 0 0 0 0 0 0 0 0 1 0 0
106 487 0 0 0 0 0 0 0 0 0 1 0
107 452 0 0 0 0 0 0 0 0 0 0 1
108 391 0 0 0 0 0 0 0 0 0 0 0
109 500 1 0 0 0 0 0 0 0 0 0 0
110 451 0 1 0 0 0 0 0 0 0 0 0
111 375 0 0 1 0 0 0 0 0 0 0 0
112 372 0 0 0 1 0 0 0 0 0 0 0
113 302 0 0 0 0 1 0 0 0 0 0 0
114 316 0 0 0 0 0 1 0 0 0 0 0
115 398 0 0 0 0 0 0 1 0 0 0 0
116 394 0 0 0 0 0 0 0 1 0 0 0
117 431 0 0 0 0 0 0 0 0 1 0 0
118 431 0 0 0 0 0 0 0 0 0 1 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
210.889 -18.789 -29.589 -28.889 -44.789 -49.589
M6 M7 M8 M9 M10 M11
-46.889 6.211 26.611 7.811 12.811 -8.000
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-202.50 -115.03 -30.89 103.97 307.90
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 210.889 43.987 4.794 5.35e-06 ***
M1 -18.789 60.631 -0.310 0.757
M2 -29.589 60.631 -0.488 0.627
M3 -28.889 60.631 -0.476 0.635
M4 -44.789 60.631 -0.739 0.462
M5 -49.589 60.631 -0.818 0.415
M6 -46.889 60.631 -0.773 0.441
M7 6.211 60.631 0.102 0.919
M8 26.611 60.631 0.439 0.662
M9 7.811 60.631 0.129 0.898
M10 12.811 60.631 0.211 0.833
M11 -8.000 62.207 -0.129 0.898
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 132 on 106 degrees of freedom
Multiple R-squared: 0.03775, Adjusted R-squared: -0.06211
F-statistic: 0.378 on 11 and 106 DF, p-value: 0.962
> 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,] 7.842632e-04 1.568526e-03 0.99921574
[2,] 6.533532e-05 1.306706e-04 0.99993466
[3,] 4.550915e-06 9.101831e-06 0.99999545
[4,] 3.850224e-07 7.700448e-07 0.99999961
[5,] 3.444106e-08 6.888213e-08 0.99999997
[6,] 1.892971e-08 3.785941e-08 0.99999998
[7,] 3.700644e-09 7.401289e-09 1.00000000
[8,] 6.891523e-10 1.378305e-09 1.00000000
[9,] 1.475201e-09 2.950401e-09 1.00000000
[10,] 1.471745e-09 2.943490e-09 1.00000000
[11,] 6.416233e-10 1.283247e-09 1.00000000
[12,] 1.197043e-10 2.394086e-10 1.00000000
[13,] 1.812023e-11 3.624047e-11 1.00000000
[14,] 4.559606e-11 9.119212e-11 1.00000000
[15,] 1.007317e-10 2.014634e-10 1.00000000
[16,] 3.524054e-11 7.048107e-11 1.00000000
[17,] 9.431298e-11 1.886260e-10 1.00000000
[18,] 1.650756e-09 3.301512e-09 1.00000000
[19,] 3.473085e-09 6.946169e-09 1.00000000
[20,] 3.295411e-09 6.590821e-09 1.00000000
[21,] 6.306849e-09 1.261370e-08 0.99999999
[22,] 1.105640e-08 2.211279e-08 0.99999999
[23,] 5.341014e-08 1.068203e-07 0.99999995
[24,] 6.145262e-08 1.229052e-07 0.99999994
[25,] 9.617381e-08 1.923476e-07 0.99999990
[26,] 7.329352e-08 1.465870e-07 0.99999993
[27,] 5.166997e-08 1.033399e-07 0.99999995
[28,] 2.809932e-08 5.619864e-08 0.99999997
[29,] 2.600156e-08 5.200312e-08 0.99999997
[30,] 2.882089e-08 5.764177e-08 0.99999997
[31,] 2.516004e-08 5.032007e-08 0.99999997
[32,] 2.746317e-08 5.492633e-08 0.99999997
[33,] 2.320119e-08 4.640237e-08 0.99999998
[34,] 2.502309e-08 5.004619e-08 0.99999997
[35,] 2.346341e-08 4.692683e-08 0.99999998
[36,] 2.353668e-08 4.707335e-08 0.99999998
[37,] 2.278465e-08 4.556931e-08 0.99999998
[38,] 2.092325e-08 4.184649e-08 0.99999998
[39,] 1.677722e-08 3.355445e-08 0.99999998
[40,] 2.642384e-08 5.284768e-08 0.99999997
[41,] 7.953141e-08 1.590628e-07 0.99999992
[42,] 4.327134e-07 8.654268e-07 0.99999957
[43,] 6.410178e-06 1.282036e-05 0.99999359
[44,] 3.379425e-05 6.758850e-05 0.99996621
[45,] 8.357842e-05 1.671568e-04 0.99991642
[46,] 1.853818e-04 3.707637e-04 0.99981462
[47,] 4.713233e-04 9.426466e-04 0.99952868
[48,] 9.969558e-04 1.993912e-03 0.99900304
[49,] 1.865413e-03 3.730826e-03 0.99813459
[50,] 3.258953e-03 6.517906e-03 0.99674105
[51,] 4.867568e-03 9.735136e-03 0.99513243
[52,] 7.741659e-03 1.548332e-02 0.99225834
[53,] 2.551154e-02 5.102309e-02 0.97448846
[54,] 6.394115e-02 1.278823e-01 0.93605885
[55,] 1.050789e-01 2.101578e-01 0.89492111
[56,] 2.064381e-01 4.128761e-01 0.79356194
[57,] 2.962088e-01 5.924177e-01 0.70379117
[58,] 3.858867e-01 7.717735e-01 0.61411327
[59,] 5.162870e-01 9.674260e-01 0.48371301
[60,] 6.241997e-01 7.516005e-01 0.37580026
[61,] 6.725803e-01 6.548395e-01 0.32741975
[62,] 7.329494e-01 5.341012e-01 0.26705061
[63,] 7.340439e-01 5.319122e-01 0.26595612
[64,] 7.492106e-01 5.015788e-01 0.25078941
[65,] 7.930128e-01 4.139744e-01 0.20698719
[66,] 8.380306e-01 3.239387e-01 0.16196935
[67,] 8.489806e-01 3.020389e-01 0.15101944
[68,] 9.234464e-01 1.531072e-01 0.07655361
[69,] 9.551166e-01 8.976683e-02 0.04488341
[70,] 9.548985e-01 9.020305e-02 0.04510152
[71,] 9.705915e-01 5.881694e-02 0.02940847
[72,] 9.801703e-01 3.965930e-02 0.01982965
[73,] 9.762255e-01 4.754908e-02 0.02377454
[74,] 9.716193e-01 5.676148e-02 0.02838074
[75,] 9.631305e-01 7.373903e-02 0.03686952
[76,] 9.520162e-01 9.596753e-02 0.04798377
[77,] 9.450158e-01 1.099685e-01 0.05498424
[78,] 9.363281e-01 1.273438e-01 0.06367192
[79,] 9.340106e-01 1.319788e-01 0.06598939
[80,] 9.571669e-01 8.566618e-02 0.04283309
[81,] 9.573463e-01 8.530738e-02 0.04265369
[82,] 9.379780e-01 1.240441e-01 0.06202203
[83,] 9.739444e-01 5.211115e-02 0.02605557
[84,] 9.827690e-01 3.446206e-02 0.01723103
[85,] 9.801173e-01 3.976537e-02 0.01988268
[86,] 9.915689e-01 1.686222e-02 0.00843111
[87,] 9.855706e-01 2.885885e-02 0.01442942
[88,] 9.626405e-01 7.471896e-02 0.03735948
[89,] 9.161617e-01 1.676767e-01 0.08383835
> postscript(file="/var/wessaorg/rcomp/tmp/1pjgs1356033713.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/2g49k1356033713.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/3yrch1356033713.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/4nwif1356033713.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/5vker1356033713.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
-151.10000 -142.30000 -132.00000 -126.10000 -118.30000 -126.00000 -173.10000
8 9 10 11 12 13 14
-202.50000 -179.70000 -188.70000 -173.88889 -161.88889 -142.10000 -122.30000
15 16 17 18 19 20 21
-119.00000 -134.10000 -122.30000 -117.00000 -164.10000 -177.50000 -161.70000
22 23 24 25 26 27 28
-171.70000 -132.88889 -120.88889 -118.10000 -119.30000 -127.00000 -82.10000
29 30 31 32 33 34 35
-67.30000 -94.00000 -109.10000 -98.50000 -98.70000 -126.70000 -76.88889
36 37 38 39 40 41 42
-61.88889 -34.10000 -57.30000 -42.00000 -57.10000 -47.30000 -87.00000
43 44 45 46 47 48 49
-97.10000 -104.50000 -108.70000 -131.70000 -105.88889 -132.88889 -93.10000
50 51 52 53 54 55 56
-74.30000 -70.00000 -76.10000 -63.30000 -39.00000 -62.10000 -47.50000
57 58 59 60 61 62 63
17.30000 -34.70000 -28.88889 -32.88889 -56.10000 -20.30000 -11.00000
64 65 66 67 68 69 70
-17.10000 22.70000 -9.00000 58.90000 -13.50000 -5.70000 55.30000
71 72 73 74 75 76 77
65.11111 76.11111 45.90000 31.70000 75.00000 126.90000 50.70000
78 79 80 81 82 83 84
82.00000 135.90000 101.50000 89.30000 23.30000 54.11111 111.11111
85 86 87 88 89 90 91
105.90000 91.70000 130.00000 82.90000 124.70000 115.00000 91.90000
92 93 94 95 96 97 98
163.50000 90.30000 104.30000 150.11111 143.11111 134.90000 142.70000
99 100 101 102 103 104 105
103.00000 76.90000 79.70000 123.00000 137.90000 222.50000 145.30000
106 107 108 109 110 111 112
263.30000 249.11111 180.11111 307.90000 269.70000 193.00000 205.90000
113 114 115 116 117 118
140.70000 152.00000 180.90000 156.50000 212.30000 207.30000
> postscript(file="/var/wessaorg/rcomp/tmp/63yur1356033713.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 -151.10000 NA
1 -142.30000 -151.10000
2 -132.00000 -142.30000
3 -126.10000 -132.00000
4 -118.30000 -126.10000
5 -126.00000 -118.30000
6 -173.10000 -126.00000
7 -202.50000 -173.10000
8 -179.70000 -202.50000
9 -188.70000 -179.70000
10 -173.88889 -188.70000
11 -161.88889 -173.88889
12 -142.10000 -161.88889
13 -122.30000 -142.10000
14 -119.00000 -122.30000
15 -134.10000 -119.00000
16 -122.30000 -134.10000
17 -117.00000 -122.30000
18 -164.10000 -117.00000
19 -177.50000 -164.10000
20 -161.70000 -177.50000
21 -171.70000 -161.70000
22 -132.88889 -171.70000
23 -120.88889 -132.88889
24 -118.10000 -120.88889
25 -119.30000 -118.10000
26 -127.00000 -119.30000
27 -82.10000 -127.00000
28 -67.30000 -82.10000
29 -94.00000 -67.30000
30 -109.10000 -94.00000
31 -98.50000 -109.10000
32 -98.70000 -98.50000
33 -126.70000 -98.70000
34 -76.88889 -126.70000
35 -61.88889 -76.88889
36 -34.10000 -61.88889
37 -57.30000 -34.10000
38 -42.00000 -57.30000
39 -57.10000 -42.00000
40 -47.30000 -57.10000
41 -87.00000 -47.30000
42 -97.10000 -87.00000
43 -104.50000 -97.10000
44 -108.70000 -104.50000
45 -131.70000 -108.70000
46 -105.88889 -131.70000
47 -132.88889 -105.88889
48 -93.10000 -132.88889
49 -74.30000 -93.10000
50 -70.00000 -74.30000
51 -76.10000 -70.00000
52 -63.30000 -76.10000
53 -39.00000 -63.30000
54 -62.10000 -39.00000
55 -47.50000 -62.10000
56 17.30000 -47.50000
57 -34.70000 17.30000
58 -28.88889 -34.70000
59 -32.88889 -28.88889
60 -56.10000 -32.88889
61 -20.30000 -56.10000
62 -11.00000 -20.30000
63 -17.10000 -11.00000
64 22.70000 -17.10000
65 -9.00000 22.70000
66 58.90000 -9.00000
67 -13.50000 58.90000
68 -5.70000 -13.50000
69 55.30000 -5.70000
70 65.11111 55.30000
71 76.11111 65.11111
72 45.90000 76.11111
73 31.70000 45.90000
74 75.00000 31.70000
75 126.90000 75.00000
76 50.70000 126.90000
77 82.00000 50.70000
78 135.90000 82.00000
79 101.50000 135.90000
80 89.30000 101.50000
81 23.30000 89.30000
82 54.11111 23.30000
83 111.11111 54.11111
84 105.90000 111.11111
85 91.70000 105.90000
86 130.00000 91.70000
87 82.90000 130.00000
88 124.70000 82.90000
89 115.00000 124.70000
90 91.90000 115.00000
91 163.50000 91.90000
92 90.30000 163.50000
93 104.30000 90.30000
94 150.11111 104.30000
95 143.11111 150.11111
96 134.90000 143.11111
97 142.70000 134.90000
98 103.00000 142.70000
99 76.90000 103.00000
100 79.70000 76.90000
101 123.00000 79.70000
102 137.90000 123.00000
103 222.50000 137.90000
104 145.30000 222.50000
105 263.30000 145.30000
106 249.11111 263.30000
107 180.11111 249.11111
108 307.90000 180.11111
109 269.70000 307.90000
110 193.00000 269.70000
111 205.90000 193.00000
112 140.70000 205.90000
113 152.00000 140.70000
114 180.90000 152.00000
115 156.50000 180.90000
116 212.30000 156.50000
117 207.30000 212.30000
118 NA 207.30000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -142.30000 -151.10000
[2,] -132.00000 -142.30000
[3,] -126.10000 -132.00000
[4,] -118.30000 -126.10000
[5,] -126.00000 -118.30000
[6,] -173.10000 -126.00000
[7,] -202.50000 -173.10000
[8,] -179.70000 -202.50000
[9,] -188.70000 -179.70000
[10,] -173.88889 -188.70000
[11,] -161.88889 -173.88889
[12,] -142.10000 -161.88889
[13,] -122.30000 -142.10000
[14,] -119.00000 -122.30000
[15,] -134.10000 -119.00000
[16,] -122.30000 -134.10000
[17,] -117.00000 -122.30000
[18,] -164.10000 -117.00000
[19,] -177.50000 -164.10000
[20,] -161.70000 -177.50000
[21,] -171.70000 -161.70000
[22,] -132.88889 -171.70000
[23,] -120.88889 -132.88889
[24,] -118.10000 -120.88889
[25,] -119.30000 -118.10000
[26,] -127.00000 -119.30000
[27,] -82.10000 -127.00000
[28,] -67.30000 -82.10000
[29,] -94.00000 -67.30000
[30,] -109.10000 -94.00000
[31,] -98.50000 -109.10000
[32,] -98.70000 -98.50000
[33,] -126.70000 -98.70000
[34,] -76.88889 -126.70000
[35,] -61.88889 -76.88889
[36,] -34.10000 -61.88889
[37,] -57.30000 -34.10000
[38,] -42.00000 -57.30000
[39,] -57.10000 -42.00000
[40,] -47.30000 -57.10000
[41,] -87.00000 -47.30000
[42,] -97.10000 -87.00000
[43,] -104.50000 -97.10000
[44,] -108.70000 -104.50000
[45,] -131.70000 -108.70000
[46,] -105.88889 -131.70000
[47,] -132.88889 -105.88889
[48,] -93.10000 -132.88889
[49,] -74.30000 -93.10000
[50,] -70.00000 -74.30000
[51,] -76.10000 -70.00000
[52,] -63.30000 -76.10000
[53,] -39.00000 -63.30000
[54,] -62.10000 -39.00000
[55,] -47.50000 -62.10000
[56,] 17.30000 -47.50000
[57,] -34.70000 17.30000
[58,] -28.88889 -34.70000
[59,] -32.88889 -28.88889
[60,] -56.10000 -32.88889
[61,] -20.30000 -56.10000
[62,] -11.00000 -20.30000
[63,] -17.10000 -11.00000
[64,] 22.70000 -17.10000
[65,] -9.00000 22.70000
[66,] 58.90000 -9.00000
[67,] -13.50000 58.90000
[68,] -5.70000 -13.50000
[69,] 55.30000 -5.70000
[70,] 65.11111 55.30000
[71,] 76.11111 65.11111
[72,] 45.90000 76.11111
[73,] 31.70000 45.90000
[74,] 75.00000 31.70000
[75,] 126.90000 75.00000
[76,] 50.70000 126.90000
[77,] 82.00000 50.70000
[78,] 135.90000 82.00000
[79,] 101.50000 135.90000
[80,] 89.30000 101.50000
[81,] 23.30000 89.30000
[82,] 54.11111 23.30000
[83,] 111.11111 54.11111
[84,] 105.90000 111.11111
[85,] 91.70000 105.90000
[86,] 130.00000 91.70000
[87,] 82.90000 130.00000
[88,] 124.70000 82.90000
[89,] 115.00000 124.70000
[90,] 91.90000 115.00000
[91,] 163.50000 91.90000
[92,] 90.30000 163.50000
[93,] 104.30000 90.30000
[94,] 150.11111 104.30000
[95,] 143.11111 150.11111
[96,] 134.90000 143.11111
[97,] 142.70000 134.90000
[98,] 103.00000 142.70000
[99,] 76.90000 103.00000
[100,] 79.70000 76.90000
[101,] 123.00000 79.70000
[102,] 137.90000 123.00000
[103,] 222.50000 137.90000
[104,] 145.30000 222.50000
[105,] 263.30000 145.30000
[106,] 249.11111 263.30000
[107,] 180.11111 249.11111
[108,] 307.90000 180.11111
[109,] 269.70000 307.90000
[110,] 193.00000 269.70000
[111,] 205.90000 193.00000
[112,] 140.70000 205.90000
[113,] 152.00000 140.70000
[114,] 180.90000 152.00000
[115,] 156.50000 180.90000
[116,] 212.30000 156.50000
[117,] 207.30000 212.30000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -142.30000 -151.10000
2 -132.00000 -142.30000
3 -126.10000 -132.00000
4 -118.30000 -126.10000
5 -126.00000 -118.30000
6 -173.10000 -126.00000
7 -202.50000 -173.10000
8 -179.70000 -202.50000
9 -188.70000 -179.70000
10 -173.88889 -188.70000
11 -161.88889 -173.88889
12 -142.10000 -161.88889
13 -122.30000 -142.10000
14 -119.00000 -122.30000
15 -134.10000 -119.00000
16 -122.30000 -134.10000
17 -117.00000 -122.30000
18 -164.10000 -117.00000
19 -177.50000 -164.10000
20 -161.70000 -177.50000
21 -171.70000 -161.70000
22 -132.88889 -171.70000
23 -120.88889 -132.88889
24 -118.10000 -120.88889
25 -119.30000 -118.10000
26 -127.00000 -119.30000
27 -82.10000 -127.00000
28 -67.30000 -82.10000
29 -94.00000 -67.30000
30 -109.10000 -94.00000
31 -98.50000 -109.10000
32 -98.70000 -98.50000
33 -126.70000 -98.70000
34 -76.88889 -126.70000
35 -61.88889 -76.88889
36 -34.10000 -61.88889
37 -57.30000 -34.10000
38 -42.00000 -57.30000
39 -57.10000 -42.00000
40 -47.30000 -57.10000
41 -87.00000 -47.30000
42 -97.10000 -87.00000
43 -104.50000 -97.10000
44 -108.70000 -104.50000
45 -131.70000 -108.70000
46 -105.88889 -131.70000
47 -132.88889 -105.88889
48 -93.10000 -132.88889
49 -74.30000 -93.10000
50 -70.00000 -74.30000
51 -76.10000 -70.00000
52 -63.30000 -76.10000
53 -39.00000 -63.30000
54 -62.10000 -39.00000
55 -47.50000 -62.10000
56 17.30000 -47.50000
57 -34.70000 17.30000
58 -28.88889 -34.70000
59 -32.88889 -28.88889
60 -56.10000 -32.88889
61 -20.30000 -56.10000
62 -11.00000 -20.30000
63 -17.10000 -11.00000
64 22.70000 -17.10000
65 -9.00000 22.70000
66 58.90000 -9.00000
67 -13.50000 58.90000
68 -5.70000 -13.50000
69 55.30000 -5.70000
70 65.11111 55.30000
71 76.11111 65.11111
72 45.90000 76.11111
73 31.70000 45.90000
74 75.00000 31.70000
75 126.90000 75.00000
76 50.70000 126.90000
77 82.00000 50.70000
78 135.90000 82.00000
79 101.50000 135.90000
80 89.30000 101.50000
81 23.30000 89.30000
82 54.11111 23.30000
83 111.11111 54.11111
84 105.90000 111.11111
85 91.70000 105.90000
86 130.00000 91.70000
87 82.90000 130.00000
88 124.70000 82.90000
89 115.00000 124.70000
90 91.90000 115.00000
91 163.50000 91.90000
92 90.30000 163.50000
93 104.30000 90.30000
94 150.11111 104.30000
95 143.11111 150.11111
96 134.90000 143.11111
97 142.70000 134.90000
98 103.00000 142.70000
99 76.90000 103.00000
100 79.70000 76.90000
101 123.00000 79.70000
102 137.90000 123.00000
103 222.50000 137.90000
104 145.30000 222.50000
105 263.30000 145.30000
106 249.11111 263.30000
107 180.11111 249.11111
108 307.90000 180.11111
109 269.70000 307.90000
110 193.00000 269.70000
111 205.90000 193.00000
112 140.70000 205.90000
113 152.00000 140.70000
114 180.90000 152.00000
115 156.50000 180.90000
116 212.30000 156.50000
117 207.30000 212.30000
> 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/7ihj21356033713.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/8xycq1356033713.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/91gti1356033713.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/10su6r1356033713.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/11irdk1356033713.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/12ym1a1356033713.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/13xyzo1356033713.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/1482691356033713.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/15irbz1356033713.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/16zzjr1356033713.tab")
+ }
>
> try(system("convert tmp/1pjgs1356033713.ps tmp/1pjgs1356033713.png",intern=TRUE))
character(0)
> try(system("convert tmp/2g49k1356033713.ps tmp/2g49k1356033713.png",intern=TRUE))
character(0)
> try(system("convert tmp/3yrch1356033713.ps tmp/3yrch1356033713.png",intern=TRUE))
character(0)
> try(system("convert tmp/4nwif1356033713.ps tmp/4nwif1356033713.png",intern=TRUE))
character(0)
> try(system("convert tmp/5vker1356033713.ps tmp/5vker1356033713.png",intern=TRUE))
character(0)
> try(system("convert tmp/63yur1356033713.ps tmp/63yur1356033713.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ihj21356033713.ps tmp/7ihj21356033713.png",intern=TRUE))
character(0)
> try(system("convert tmp/8xycq1356033713.ps tmp/8xycq1356033713.png",intern=TRUE))
character(0)
> try(system("convert tmp/91gti1356033713.ps tmp/91gti1356033713.png",intern=TRUE))
character(0)
> try(system("convert tmp/10su6r1356033713.ps tmp/10su6r1356033713.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
7.064 1.043 8.123