R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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(4.031636
+ ,0.5215052
+ ,3.702076
+ ,0.4248284
+ ,3.056176
+ ,0.4250311
+ ,3.280707
+ ,0.4771938
+ ,2.984728
+ ,0.8280212
+ ,3.693712
+ ,0.6156186
+ ,3.226317
+ ,0.366627
+ ,2.190349
+ ,0.4308883
+ ,2.599515
+ ,0.2810287
+ ,3.080288
+ ,0.4646245
+ ,2.929672
+ ,0.2693951
+ ,2.922548
+ ,0.5779049
+ ,3.234943
+ ,0.5661151
+ ,2.983081
+ ,0.5077584
+ ,3.284389
+ ,0.7507175
+ ,3.806511
+ ,0.6808395
+ ,3.784579
+ ,0.7661091
+ ,2.645654
+ ,0.4561473
+ ,3.092081
+ ,0.4977496
+ ,3.204859
+ ,0.4193273
+ ,3.107225
+ ,0.6095514
+ ,3.466909
+ ,0.457337
+ ,2.984404
+ ,0.5705478
+ ,3.218072
+ ,0.3478996
+ ,2.82731
+ ,0.3874993
+ ,3.182049
+ ,0.5824285
+ ,2.236319
+ ,0.2391033
+ ,2.033218
+ ,0.2367445
+ ,1.644804
+ ,0.2626158
+ ,1.627971
+ ,0.4240934
+ ,1.677559
+ ,0.365275
+ ,2.330828
+ ,0.3750758
+ ,2.493615
+ ,0.4090056
+ ,2.257172
+ ,0.3891676
+ ,2.655517
+ ,0.240261
+ ,2.298655
+ ,0.1589496
+ ,2.600402
+ ,0.4393373
+ ,3.04523
+ ,0.5094681
+ ,2.790583
+ ,0.3743465
+ ,3.227052
+ ,0.4339828
+ ,2.967479
+ ,0.4130557
+ ,2.938817
+ ,0.3288928
+ ,3.277961
+ ,0.5186648
+ ,3.423985
+ ,0.5486504
+ ,3.072646
+ ,0.5469111
+ ,2.754253
+ ,0.4963494
+ ,2.910431
+ ,0.5308929
+ ,3.174369
+ ,0.5957761
+ ,3.068387
+ ,0.5570584
+ ,3.089543
+ ,0.5731325
+ ,2.906654
+ ,0.5005416
+ ,2.931161
+ ,0.5431269
+ ,3.02566
+ ,0.5593657
+ ,2.939551
+ ,0.6911693
+ ,2.691019
+ ,0.4403485
+ ,3.19812
+ ,0.5676662
+ ,3.07639
+ ,0.5969114
+ ,2.863873
+ ,0.4735537
+ ,3.013802
+ ,0.5923935
+ ,3.053364
+ ,0.5975556
+ ,2.864753
+ ,0.6334127
+ ,3.057062
+ ,0.6057115
+ ,2.959365
+ ,0.7046107
+ ,3.252258
+ ,0.4805263
+ ,3.602988
+ ,0.702686
+ ,3.497704
+ ,0.7009017
+ ,3.296867
+ ,0.6030854
+ ,3.602417
+ ,0.6980919
+ ,3.3001
+ ,0.597656
+ ,3.40193
+ ,0.8023421
+ ,3.502591
+ ,0.6017109
+ ,3.402348
+ ,0.5993127
+ ,3.498551
+ ,0.6025625
+ ,3.199823
+ ,0.7016625
+ ,2.700064
+ ,0.4995714
+ ,2.801034
+ ,0.4980918
+ ,2.898628
+ ,0.497569
+ ,2.800854
+ ,0.600183
+ ,2.399942
+ ,0.3339542
+ ,2.402724
+ ,0.274437
+ ,2.202331
+ ,0.3209428
+ ,2.102594
+ ,0.5406671
+ ,1.798293
+ ,0.4050209
+ ,1.202484
+ ,0.2885961
+ ,1.400201
+ ,0.3275942
+ ,1.200832
+ ,0.3132606
+ ,1.298083
+ ,0.2575562
+ ,1.099742
+ ,0.2138386
+ ,1.001377
+ ,0.1861856
+ ,0.8361743
+ ,0.1592713)
+ ,dim=c(2
+ ,90)
+ ,dimnames=list(c('firearmsuicide'
+ ,'firearmhomicide')
+ ,1:90))
> y <- array(NA,dim=c(2,90),dimnames=list(c('firearmsuicide','firearmhomicide'),1:90))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '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.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
firearmsuicide firearmhomicide t
1 4.0316360 0.5215052 1
2 3.7020760 0.4248284 2
3 3.0561760 0.4250311 3
4 3.2807070 0.4771938 4
5 2.9847280 0.8280212 5
6 3.6937120 0.6156186 6
7 3.2263170 0.3666270 7
8 2.1903490 0.4308883 8
9 2.5995150 0.2810287 9
10 3.0802880 0.4646245 10
11 2.9296720 0.2693951 11
12 2.9225480 0.5779049 12
13 3.2349430 0.5661151 13
14 2.9830810 0.5077584 14
15 3.2843890 0.7507175 15
16 3.8065110 0.6808395 16
17 3.7845790 0.7661091 17
18 2.6456540 0.4561473 18
19 3.0920810 0.4977496 19
20 3.2048590 0.4193273 20
21 3.1072250 0.6095514 21
22 3.4669090 0.4573370 22
23 2.9844040 0.5705478 23
24 3.2180720 0.3478996 24
25 2.8273100 0.3874993 25
26 3.1820490 0.5824285 26
27 2.2363190 0.2391033 27
28 2.0332180 0.2367445 28
29 1.6448040 0.2626158 29
30 1.6279710 0.4240934 30
31 1.6775590 0.3652750 31
32 2.3308280 0.3750758 32
33 2.4936150 0.4090056 33
34 2.2571720 0.3891676 34
35 2.6555170 0.2402610 35
36 2.2986550 0.1589496 36
37 2.6004020 0.4393373 37
38 3.0452300 0.5094681 38
39 2.7905830 0.3743465 39
40 3.2270520 0.4339828 40
41 2.9674790 0.4130557 41
42 2.9388170 0.3288928 42
43 3.2779610 0.5186648 43
44 3.4239850 0.5486504 44
45 3.0726460 0.5469111 45
46 2.7542530 0.4963494 46
47 2.9104310 0.5308929 47
48 3.1743690 0.5957761 48
49 3.0683870 0.5570584 49
50 3.0895430 0.5731325 50
51 2.9066540 0.5005416 51
52 2.9311610 0.5431269 52
53 3.0256600 0.5593657 53
54 2.9395510 0.6911693 54
55 2.6910190 0.4403485 55
56 3.1981200 0.5676662 56
57 3.0763900 0.5969114 57
58 2.8638730 0.4735537 58
59 3.0138020 0.5923935 59
60 3.0533640 0.5975556 60
61 2.8647530 0.6334127 61
62 3.0570620 0.6057115 62
63 2.9593650 0.7046107 63
64 3.2522580 0.4805263 64
65 3.6029880 0.7026860 65
66 3.4977040 0.7009017 66
67 3.2968670 0.6030854 67
68 3.6024170 0.6980919 68
69 3.3001000 0.5976560 69
70 3.4019300 0.8023421 70
71 3.5025910 0.6017109 71
72 3.4023480 0.5993127 72
73 3.4985510 0.6025625 73
74 3.1998230 0.7016625 74
75 2.7000640 0.4995714 75
76 2.8010340 0.4980918 76
77 2.8986280 0.4975690 77
78 2.8008540 0.6001830 78
79 2.3999420 0.3339542 79
80 2.4027240 0.2744370 80
81 2.2023310 0.3209428 81
82 2.1025940 0.5406671 82
83 1.7982930 0.4050209 83
84 1.2024840 0.2885961 84
85 1.4002010 0.3275942 85
86 1.2008320 0.3132606 86
87 1.2980830 0.2575562 87
88 1.0997420 0.2138386 88
89 1.0013770 0.1861856 89
90 0.8361743 0.1592713 90
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) firearmhomicide t
1.752194 3.064703 -0.009476
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.25772 -0.26337 0.04479 0.33627 0.69066
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.752194 0.171696 10.205 < 2e-16 ***
firearmhomicide 3.064703 0.294062 10.422 < 2e-16 ***
t -0.009476 0.001709 -5.543 3.13e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4207 on 87 degrees of freedom
Multiple R-squared: 0.6267, Adjusted R-squared: 0.6181
F-statistic: 73.03 on 2 and 87 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.7580024 0.4839951872 0.2419975936
[2,] 0.6262328 0.7475344713 0.3737672356
[3,] 0.7964581 0.4070838809 0.2035419405
[4,] 0.7008912 0.5982176412 0.2991088206
[5,] 0.7428511 0.5142978042 0.2571489021
[6,] 0.7261389 0.5477222791 0.2738611395
[7,] 0.6878066 0.6243868746 0.3121934373
[8,] 0.7135270 0.5729460475 0.2864730238
[9,] 0.6537057 0.6925887000 0.3462943500
[10,] 0.6546153 0.6907693240 0.3453846620
[11,] 0.7776189 0.4447621009 0.2223810505
[12,] 0.7840619 0.4318762173 0.2159381086
[13,] 0.7502898 0.4994204681 0.2497102341
[14,] 0.6937751 0.6124497701 0.3062248850
[15,] 0.6736252 0.6527496170 0.3263748085
[16,] 0.6232876 0.7534248847 0.3767124423
[17,] 0.6571296 0.6857408876 0.3428704438
[18,] 0.6129744 0.7740511049 0.3870255524
[19,] 0.6269422 0.7461155387 0.3730577693
[20,] 0.5645151 0.8709697060 0.4354848530
[21,] 0.5008806 0.9982388711 0.4991194355
[22,] 0.4997602 0.9995203469 0.5002398265
[23,] 0.5140869 0.9718261255 0.4859130628
[24,] 0.6607856 0.6784288898 0.3392144449
[25,] 0.9120585 0.1758829515 0.0879414758
[26,] 0.9726316 0.0547368050 0.0273684025
[27,] 0.9694851 0.0610297844 0.0305148922
[28,] 0.9666541 0.0666918260 0.0333459130
[29,] 0.9714511 0.0570977996 0.0285488998
[30,] 0.9771694 0.0456611161 0.0228305581
[31,] 0.9739310 0.0521380884 0.0260690442
[32,] 0.9717291 0.0565418827 0.0282709414
[33,] 0.9722294 0.0555412352 0.0277706176
[34,] 0.9685365 0.0629269909 0.0314634954
[35,] 0.9774382 0.0451236824 0.0225618412
[36,] 0.9747795 0.0504410222 0.0252205111
[37,] 0.9807421 0.0385158065 0.0192579033
[38,] 0.9794150 0.0411699715 0.0205849858
[39,] 0.9795241 0.0409517270 0.0204758635
[40,] 0.9712448 0.0575103855 0.0287551928
[41,] 0.9613441 0.0773117539 0.0386558770
[42,] 0.9483211 0.1033578346 0.0516789173
[43,] 0.9326272 0.1347456621 0.0673728310
[44,] 0.9118557 0.1762886969 0.0881443485
[45,] 0.8873848 0.2252303599 0.1126151800
[46,] 0.8553046 0.2893907541 0.1446953770
[47,] 0.8225911 0.3548178968 0.1774089484
[48,] 0.7840375 0.4319249179 0.2159624589
[49,] 0.8488979 0.3022042939 0.1511021469
[50,] 0.8134076 0.3731847659 0.1865923829
[51,] 0.7757954 0.4484092271 0.2242046135
[52,] 0.7482560 0.5034880968 0.2517440484
[53,] 0.7006283 0.5987434741 0.2993717370
[54,] 0.6852099 0.6295802472 0.3147901236
[55,] 0.6764124 0.6471752747 0.3235876374
[56,] 0.7980789 0.4038421901 0.2019210951
[57,] 0.8501600 0.2996800025 0.1498400013
[58,] 0.9814045 0.0371910521 0.0185955260
[59,] 0.9862104 0.0275792447 0.0137896224
[60,] 0.9880867 0.0238266437 0.0119133219
[61,] 0.9919209 0.0161582157 0.0080791079
[62,] 0.9974204 0.0051592879 0.0025796440
[63,] 0.9971441 0.0057118557 0.0028559278
[64,] 0.9990818 0.0018363824 0.0009181912
[65,] 0.9995473 0.0009054117 0.0004527058
[66,] 0.9991958 0.0016084194 0.0008042097
[67,] 0.9984387 0.0031226798 0.0015613399
[68,] 0.9977243 0.0045514216 0.0022757108
[69,] 0.9951941 0.0096118772 0.0048059386
[70,] 0.9969230 0.0061539532 0.0030769766
[71,] 0.9946708 0.0106583487 0.0053291744
[72,] 0.9894992 0.0210015891 0.0105007946
[73,] 0.9842921 0.0314157806 0.0157078903
[74,] 0.9671384 0.0657231953 0.0328615977
[75,] 0.9585829 0.0828341688 0.0414170844
[76,] 0.9973435 0.0053129820 0.0026564910
[77,] 0.9921754 0.0156492222 0.0078246111
[78,] 0.9901842 0.0196316865 0.0098158433
[79,] 0.9951220 0.0097560345 0.0048780172
> postscript(file="/var/www/html/rcomp/tmp/175v81293203072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/20xct1293203072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/30xct1293203072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/40xct1293203072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5aotw1293203072.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 = 90
Frequency = 1
1 2 3 4 5 6
0.690659278 0.666861273 0.029816326 0.103960386 -1.257724304 0.111686955
7 8 9 10 11 12
0.416853649 -0.806579913 0.071337595 -0.001079826 0.456100666 -0.487038127
13 14 15 16 17 18
-0.129034617 -0.192574366 -0.626387700 0.119365919 -0.154415854 -0.333923574
19 20 21 22 23 24
-0.005519020 0.357076345 -0.314061850 0.521590422 -0.298395844 0.627099140
25 26 27 28 29 30
0.124452070 -0.108732862 0.007203346 -0.179192363 -0.637417958 -1.139655654
31 32 33 34 35 36
-0.900330430 -0.267621707 -0.199343216 -0.365512359 0.498663486 0.400473086
37 38 39 40 41 42
-0.147608809 0.091765352 0.260702260 0.523879951 0.337918576 0.576667178
43 44 45 46 47 48
0.343692535 0.407295831 0.070763538 -0.083196575 -0.023407892 0.051158607
49 50 51 52 53 54
0.073311146 0.054681064 0.103737917 0.007209868 0.061418029 -0.419153656
55 56 57 58 59 60
0.110481994 0.236868263 0.034986665 0.210000707 0.005197225 0.038415188
61 62 63 64 65 66
-0.250610923 0.036070310 -0.355247146 0.633874366 0.313227026 0.222887645
67 68 69 70 71 72
0.331304869 0.355164386 0.370129907 -0.145866030 0.579146378 0.495729418
73 74 75 76 77 78
0.591449013 -0.001514834 0.127551734 0.242532538 0.351205033 -0.051574183
79 80 81 82 83 84
0.372902419 0.567563258 0.234120039 -0.529530522 -0.408639871 -0.638165111
85 86 87 88 89 90
-0.550489456 -0.696453954 -0.419009216 -0.473892466 -0.478032952 -0.551275035
> postscript(file="/var/www/html/rcomp/tmp/6aotw1293203072.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 = 90
Frequency = 1
lag(myerror, k = 1) myerror
0 0.690659278 NA
1 0.666861273 0.690659278
2 0.029816326 0.666861273
3 0.103960386 0.029816326
4 -1.257724304 0.103960386
5 0.111686955 -1.257724304
6 0.416853649 0.111686955
7 -0.806579913 0.416853649
8 0.071337595 -0.806579913
9 -0.001079826 0.071337595
10 0.456100666 -0.001079826
11 -0.487038127 0.456100666
12 -0.129034617 -0.487038127
13 -0.192574366 -0.129034617
14 -0.626387700 -0.192574366
15 0.119365919 -0.626387700
16 -0.154415854 0.119365919
17 -0.333923574 -0.154415854
18 -0.005519020 -0.333923574
19 0.357076345 -0.005519020
20 -0.314061850 0.357076345
21 0.521590422 -0.314061850
22 -0.298395844 0.521590422
23 0.627099140 -0.298395844
24 0.124452070 0.627099140
25 -0.108732862 0.124452070
26 0.007203346 -0.108732862
27 -0.179192363 0.007203346
28 -0.637417958 -0.179192363
29 -1.139655654 -0.637417958
30 -0.900330430 -1.139655654
31 -0.267621707 -0.900330430
32 -0.199343216 -0.267621707
33 -0.365512359 -0.199343216
34 0.498663486 -0.365512359
35 0.400473086 0.498663486
36 -0.147608809 0.400473086
37 0.091765352 -0.147608809
38 0.260702260 0.091765352
39 0.523879951 0.260702260
40 0.337918576 0.523879951
41 0.576667178 0.337918576
42 0.343692535 0.576667178
43 0.407295831 0.343692535
44 0.070763538 0.407295831
45 -0.083196575 0.070763538
46 -0.023407892 -0.083196575
47 0.051158607 -0.023407892
48 0.073311146 0.051158607
49 0.054681064 0.073311146
50 0.103737917 0.054681064
51 0.007209868 0.103737917
52 0.061418029 0.007209868
53 -0.419153656 0.061418029
54 0.110481994 -0.419153656
55 0.236868263 0.110481994
56 0.034986665 0.236868263
57 0.210000707 0.034986665
58 0.005197225 0.210000707
59 0.038415188 0.005197225
60 -0.250610923 0.038415188
61 0.036070310 -0.250610923
62 -0.355247146 0.036070310
63 0.633874366 -0.355247146
64 0.313227026 0.633874366
65 0.222887645 0.313227026
66 0.331304869 0.222887645
67 0.355164386 0.331304869
68 0.370129907 0.355164386
69 -0.145866030 0.370129907
70 0.579146378 -0.145866030
71 0.495729418 0.579146378
72 0.591449013 0.495729418
73 -0.001514834 0.591449013
74 0.127551734 -0.001514834
75 0.242532538 0.127551734
76 0.351205033 0.242532538
77 -0.051574183 0.351205033
78 0.372902419 -0.051574183
79 0.567563258 0.372902419
80 0.234120039 0.567563258
81 -0.529530522 0.234120039
82 -0.408639871 -0.529530522
83 -0.638165111 -0.408639871
84 -0.550489456 -0.638165111
85 -0.696453954 -0.550489456
86 -0.419009216 -0.696453954
87 -0.473892466 -0.419009216
88 -0.478032952 -0.473892466
89 -0.551275035 -0.478032952
90 NA -0.551275035
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.666861273 0.690659278
[2,] 0.029816326 0.666861273
[3,] 0.103960386 0.029816326
[4,] -1.257724304 0.103960386
[5,] 0.111686955 -1.257724304
[6,] 0.416853649 0.111686955
[7,] -0.806579913 0.416853649
[8,] 0.071337595 -0.806579913
[9,] -0.001079826 0.071337595
[10,] 0.456100666 -0.001079826
[11,] -0.487038127 0.456100666
[12,] -0.129034617 -0.487038127
[13,] -0.192574366 -0.129034617
[14,] -0.626387700 -0.192574366
[15,] 0.119365919 -0.626387700
[16,] -0.154415854 0.119365919
[17,] -0.333923574 -0.154415854
[18,] -0.005519020 -0.333923574
[19,] 0.357076345 -0.005519020
[20,] -0.314061850 0.357076345
[21,] 0.521590422 -0.314061850
[22,] -0.298395844 0.521590422
[23,] 0.627099140 -0.298395844
[24,] 0.124452070 0.627099140
[25,] -0.108732862 0.124452070
[26,] 0.007203346 -0.108732862
[27,] -0.179192363 0.007203346
[28,] -0.637417958 -0.179192363
[29,] -1.139655654 -0.637417958
[30,] -0.900330430 -1.139655654
[31,] -0.267621707 -0.900330430
[32,] -0.199343216 -0.267621707
[33,] -0.365512359 -0.199343216
[34,] 0.498663486 -0.365512359
[35,] 0.400473086 0.498663486
[36,] -0.147608809 0.400473086
[37,] 0.091765352 -0.147608809
[38,] 0.260702260 0.091765352
[39,] 0.523879951 0.260702260
[40,] 0.337918576 0.523879951
[41,] 0.576667178 0.337918576
[42,] 0.343692535 0.576667178
[43,] 0.407295831 0.343692535
[44,] 0.070763538 0.407295831
[45,] -0.083196575 0.070763538
[46,] -0.023407892 -0.083196575
[47,] 0.051158607 -0.023407892
[48,] 0.073311146 0.051158607
[49,] 0.054681064 0.073311146
[50,] 0.103737917 0.054681064
[51,] 0.007209868 0.103737917
[52,] 0.061418029 0.007209868
[53,] -0.419153656 0.061418029
[54,] 0.110481994 -0.419153656
[55,] 0.236868263 0.110481994
[56,] 0.034986665 0.236868263
[57,] 0.210000707 0.034986665
[58,] 0.005197225 0.210000707
[59,] 0.038415188 0.005197225
[60,] -0.250610923 0.038415188
[61,] 0.036070310 -0.250610923
[62,] -0.355247146 0.036070310
[63,] 0.633874366 -0.355247146
[64,] 0.313227026 0.633874366
[65,] 0.222887645 0.313227026
[66,] 0.331304869 0.222887645
[67,] 0.355164386 0.331304869
[68,] 0.370129907 0.355164386
[69,] -0.145866030 0.370129907
[70,] 0.579146378 -0.145866030
[71,] 0.495729418 0.579146378
[72,] 0.591449013 0.495729418
[73,] -0.001514834 0.591449013
[74,] 0.127551734 -0.001514834
[75,] 0.242532538 0.127551734
[76,] 0.351205033 0.242532538
[77,] -0.051574183 0.351205033
[78,] 0.372902419 -0.051574183
[79,] 0.567563258 0.372902419
[80,] 0.234120039 0.567563258
[81,] -0.529530522 0.234120039
[82,] -0.408639871 -0.529530522
[83,] -0.638165111 -0.408639871
[84,] -0.550489456 -0.638165111
[85,] -0.696453954 -0.550489456
[86,] -0.419009216 -0.696453954
[87,] -0.473892466 -0.419009216
[88,] -0.478032952 -0.473892466
[89,] -0.551275035 -0.478032952
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.666861273 0.690659278
2 0.029816326 0.666861273
3 0.103960386 0.029816326
4 -1.257724304 0.103960386
5 0.111686955 -1.257724304
6 0.416853649 0.111686955
7 -0.806579913 0.416853649
8 0.071337595 -0.806579913
9 -0.001079826 0.071337595
10 0.456100666 -0.001079826
11 -0.487038127 0.456100666
12 -0.129034617 -0.487038127
13 -0.192574366 -0.129034617
14 -0.626387700 -0.192574366
15 0.119365919 -0.626387700
16 -0.154415854 0.119365919
17 -0.333923574 -0.154415854
18 -0.005519020 -0.333923574
19 0.357076345 -0.005519020
20 -0.314061850 0.357076345
21 0.521590422 -0.314061850
22 -0.298395844 0.521590422
23 0.627099140 -0.298395844
24 0.124452070 0.627099140
25 -0.108732862 0.124452070
26 0.007203346 -0.108732862
27 -0.179192363 0.007203346
28 -0.637417958 -0.179192363
29 -1.139655654 -0.637417958
30 -0.900330430 -1.139655654
31 -0.267621707 -0.900330430
32 -0.199343216 -0.267621707
33 -0.365512359 -0.199343216
34 0.498663486 -0.365512359
35 0.400473086 0.498663486
36 -0.147608809 0.400473086
37 0.091765352 -0.147608809
38 0.260702260 0.091765352
39 0.523879951 0.260702260
40 0.337918576 0.523879951
41 0.576667178 0.337918576
42 0.343692535 0.576667178
43 0.407295831 0.343692535
44 0.070763538 0.407295831
45 -0.083196575 0.070763538
46 -0.023407892 -0.083196575
47 0.051158607 -0.023407892
48 0.073311146 0.051158607
49 0.054681064 0.073311146
50 0.103737917 0.054681064
51 0.007209868 0.103737917
52 0.061418029 0.007209868
53 -0.419153656 0.061418029
54 0.110481994 -0.419153656
55 0.236868263 0.110481994
56 0.034986665 0.236868263
57 0.210000707 0.034986665
58 0.005197225 0.210000707
59 0.038415188 0.005197225
60 -0.250610923 0.038415188
61 0.036070310 -0.250610923
62 -0.355247146 0.036070310
63 0.633874366 -0.355247146
64 0.313227026 0.633874366
65 0.222887645 0.313227026
66 0.331304869 0.222887645
67 0.355164386 0.331304869
68 0.370129907 0.355164386
69 -0.145866030 0.370129907
70 0.579146378 -0.145866030
71 0.495729418 0.579146378
72 0.591449013 0.495729418
73 -0.001514834 0.591449013
74 0.127551734 -0.001514834
75 0.242532538 0.127551734
76 0.351205033 0.242532538
77 -0.051574183 0.351205033
78 0.372902419 -0.051574183
79 0.567563258 0.372902419
80 0.234120039 0.567563258
81 -0.529530522 0.234120039
82 -0.408639871 -0.529530522
83 -0.638165111 -0.408639871
84 -0.550489456 -0.638165111
85 -0.696453954 -0.550489456
86 -0.419009216 -0.696453954
87 -0.473892466 -0.419009216
88 -0.478032952 -0.473892466
89 -0.551275035 -0.478032952
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/73fsh1293203072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/83fsh1293203072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9eoa21293203072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10eoa21293203072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/114afe1293203072.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/12kqpw1293203072.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/139qm71293203072.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14ki3s1293203072.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15niky1293203072.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/169j041293203072.tab")
+ }
>
> try(system("convert tmp/175v81293203072.ps tmp/175v81293203072.png",intern=TRUE))
character(0)
> try(system("convert tmp/20xct1293203072.ps tmp/20xct1293203072.png",intern=TRUE))
character(0)
> try(system("convert tmp/30xct1293203072.ps tmp/30xct1293203072.png",intern=TRUE))
character(0)
> try(system("convert tmp/40xct1293203072.ps tmp/40xct1293203072.png",intern=TRUE))
character(0)
> try(system("convert tmp/5aotw1293203072.ps tmp/5aotw1293203072.png",intern=TRUE))
character(0)
> try(system("convert tmp/6aotw1293203072.ps tmp/6aotw1293203072.png",intern=TRUE))
character(0)
> try(system("convert tmp/73fsh1293203072.ps tmp/73fsh1293203072.png",intern=TRUE))
character(0)
> try(system("convert tmp/83fsh1293203072.ps tmp/83fsh1293203072.png",intern=TRUE))
character(0)
> try(system("convert tmp/9eoa21293203072.ps tmp/9eoa21293203072.png",intern=TRUE))
character(0)
> try(system("convert tmp/10eoa21293203072.ps tmp/10eoa21293203072.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.821 1.679 7.276