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(37,0,30,1,47,0,35,0,30,1,43,0,82,0,40,0,47,0,19,1,52,0,136,0,80,0,42,0,54,0,66,0,81,0,63,0,137,0,72,0,107,0,58,0,36,0,52,0,79,0,77,0,54,0,84,0,48,0,96,0,83,0,66,0,61,0,53,0,30,1,74,0,69,0,59,0,42,0,65,0,70,0,100,0,63,0,105,0,82,0,81,0,75,0,102,0,121,0,98,0,76,0,77,0,63,0,37,0,35,0,23,1,40,0,29,1,37,0,51,0,20,1,28,1,13,1,22,1,25,1,13,1,16,1,13,1,16,1,17,1,9,1,17,1,25,1,14,1,8,1,7,1,10,1,7,1,10,1,3,1),dim=c(2,80),dimnames=list(c('Soldiers','Dummy'),1:80))
> y <- array(NA,dim=c(2,80),dimnames=list(c('Soldiers','Dummy'),1:80))
> 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 = '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.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
Soldiers Dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 37 0 1 0 0 0 0 0 0 0 0 0 0 1
2 30 1 0 1 0 0 0 0 0 0 0 0 0 2
3 47 0 0 0 1 0 0 0 0 0 0 0 0 3
4 35 0 0 0 0 1 0 0 0 0 0 0 0 4
5 30 1 0 0 0 0 1 0 0 0 0 0 0 5
6 43 0 0 0 0 0 0 1 0 0 0 0 0 6
7 82 0 0 0 0 0 0 0 1 0 0 0 0 7
8 40 0 0 0 0 0 0 0 0 1 0 0 0 8
9 47 0 0 0 0 0 0 0 0 0 1 0 0 9
10 19 1 0 0 0 0 0 0 0 0 0 1 0 10
11 52 0 0 0 0 0 0 0 0 0 0 0 1 11
12 136 0 0 0 0 0 0 0 0 0 0 0 0 12
13 80 0 1 0 0 0 0 0 0 0 0 0 0 13
14 42 0 0 1 0 0 0 0 0 0 0 0 0 14
15 54 0 0 0 1 0 0 0 0 0 0 0 0 15
16 66 0 0 0 0 1 0 0 0 0 0 0 0 16
17 81 0 0 0 0 0 1 0 0 0 0 0 0 17
18 63 0 0 0 0 0 0 1 0 0 0 0 0 18
19 137 0 0 0 0 0 0 0 1 0 0 0 0 19
20 72 0 0 0 0 0 0 0 0 1 0 0 0 20
21 107 0 0 0 0 0 0 0 0 0 1 0 0 21
22 58 0 0 0 0 0 0 0 0 0 0 1 0 22
23 36 0 0 0 0 0 0 0 0 0 0 0 1 23
24 52 0 0 0 0 0 0 0 0 0 0 0 0 24
25 79 0 1 0 0 0 0 0 0 0 0 0 0 25
26 77 0 0 1 0 0 0 0 0 0 0 0 0 26
27 54 0 0 0 1 0 0 0 0 0 0 0 0 27
28 84 0 0 0 0 1 0 0 0 0 0 0 0 28
29 48 0 0 0 0 0 1 0 0 0 0 0 0 29
30 96 0 0 0 0 0 0 1 0 0 0 0 0 30
31 83 0 0 0 0 0 0 0 1 0 0 0 0 31
32 66 0 0 0 0 0 0 0 0 1 0 0 0 32
33 61 0 0 0 0 0 0 0 0 0 1 0 0 33
34 53 0 0 0 0 0 0 0 0 0 0 1 0 34
35 30 1 0 0 0 0 0 0 0 0 0 0 1 35
36 74 0 0 0 0 0 0 0 0 0 0 0 0 36
37 69 0 1 0 0 0 0 0 0 0 0 0 0 37
38 59 0 0 1 0 0 0 0 0 0 0 0 0 38
39 42 0 0 0 1 0 0 0 0 0 0 0 0 39
40 65 0 0 0 0 1 0 0 0 0 0 0 0 40
41 70 0 0 0 0 0 1 0 0 0 0 0 0 41
42 100 0 0 0 0 0 0 1 0 0 0 0 0 42
43 63 0 0 0 0 0 0 0 1 0 0 0 0 43
44 105 0 0 0 0 0 0 0 0 1 0 0 0 44
45 82 0 0 0 0 0 0 0 0 0 1 0 0 45
46 81 0 0 0 0 0 0 0 0 0 0 1 0 46
47 75 0 0 0 0 0 0 0 0 0 0 0 1 47
48 102 0 0 0 0 0 0 0 0 0 0 0 0 48
49 121 0 1 0 0 0 0 0 0 0 0 0 0 49
50 98 0 0 1 0 0 0 0 0 0 0 0 0 50
51 76 0 0 0 1 0 0 0 0 0 0 0 0 51
52 77 0 0 0 0 1 0 0 0 0 0 0 0 52
53 63 0 0 0 0 0 1 0 0 0 0 0 0 53
54 37 0 0 0 0 0 0 1 0 0 0 0 0 54
55 35 0 0 0 0 0 0 0 1 0 0 0 0 55
56 23 1 0 0 0 0 0 0 0 1 0 0 0 56
57 40 0 0 0 0 0 0 0 0 0 1 0 0 57
58 29 1 0 0 0 0 0 0 0 0 0 1 0 58
59 37 0 0 0 0 0 0 0 0 0 0 0 1 59
60 51 0 0 0 0 0 0 0 0 0 0 0 0 60
61 20 1 1 0 0 0 0 0 0 0 0 0 0 61
62 28 1 0 1 0 0 0 0 0 0 0 0 0 62
63 13 1 0 0 1 0 0 0 0 0 0 0 0 63
64 22 1 0 0 0 1 0 0 0 0 0 0 0 64
65 25 1 0 0 0 0 1 0 0 0 0 0 0 65
66 13 1 0 0 0 0 0 1 0 0 0 0 0 66
67 16 1 0 0 0 0 0 0 1 0 0 0 0 67
68 13 1 0 0 0 0 0 0 0 1 0 0 0 68
69 16 1 0 0 0 0 0 0 0 0 1 0 0 69
70 17 1 0 0 0 0 0 0 0 0 0 1 0 70
71 9 1 0 0 0 0 0 0 0 0 0 0 1 71
72 17 1 0 0 0 0 0 0 0 0 0 0 0 72
73 25 1 1 0 0 0 0 0 0 0 0 0 0 73
74 14 1 0 1 0 0 0 0 0 0 0 0 0 74
75 8 1 0 0 1 0 0 0 0 0 0 0 0 75
76 7 1 0 0 0 1 0 0 0 0 0 0 0 76
77 10 1 0 0 0 0 1 0 0 0 0 0 0 77
78 7 1 0 0 0 0 0 1 0 0 0 0 0 78
79 10 1 0 0 0 0 0 0 1 0 0 0 0 79
80 3 1 0 0 0 0 0 0 0 1 0 0 0 80
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dummy M1 M2 M3 M4
82.198 -48.842 -4.859 -9.690 -24.332 -15.426
M5 M6 M7 M8 M9 M10
-12.543 -14.900 -5.279 -13.110 -13.314 -12.984
M11 t
-24.075 -0.049
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-40.290 -11.688 -1.201 7.429 61.012
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.1983 10.2343 8.032 2.32e-11 ***
Dummy -48.8422 6.7786 -7.205 6.98e-10 ***
M1 -4.8590 12.2056 -0.398 0.6918
M2 -9.6897 12.3218 -0.786 0.4345
M3 -24.3324 12.1826 -1.997 0.0499 *
M4 -15.4263 12.1733 -1.267 0.2095
M5 -12.5427 12.2745 -1.022 0.3106
M6 -14.8997 12.1593 -1.225 0.2248
M7 -5.2793 12.1545 -0.434 0.6654
M8 -13.1100 12.2404 -1.071 0.2881
M9 -13.3137 12.5969 -1.057 0.2944
M10 -12.9839 12.8239 -1.012 0.3150
M11 -24.0753 12.6492 -1.903 0.0614 .
t -0.0490 0.1350 -0.363 0.7178
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.81 on 66 degrees of freedom
Multiple R-squared: 0.61, Adjusted R-squared: 0.5332
F-statistic: 7.942 on 13 and 66 DF, p-value: 3.049e-09
> 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.3809357 7.618713e-01 6.190643e-01
[2,] 0.2347706 4.695413e-01 7.652294e-01
[3,] 0.3373336 6.746672e-01 6.626664e-01
[4,] 0.2190770 4.381540e-01 7.809230e-01
[5,] 0.2777198 5.554397e-01 7.222802e-01
[6,] 0.1968565 3.937130e-01 8.031435e-01
[7,] 0.4920700 9.841401e-01 5.079300e-01
[8,] 0.9894918 2.101635e-02 1.050817e-02
[9,] 0.9826008 3.479834e-02 1.739917e-02
[10,] 0.9716402 5.671965e-02 2.835982e-02
[11,] 0.9639241 7.215188e-02 3.607594e-02
[12,] 0.9462029 1.075942e-01 5.379710e-02
[13,] 0.9653236 6.935271e-02 3.467636e-02
[14,] 0.9597078 8.058438e-02 4.029219e-02
[15,] 0.9722473 5.550532e-02 2.775266e-02
[16,] 0.9616945 7.661100e-02 3.830550e-02
[17,] 0.9581811 8.363788e-02 4.181894e-02
[18,] 0.9580674 8.386513e-02 4.193256e-02
[19,] 0.9370663 1.258674e-01 6.293369e-02
[20,] 0.9324645 1.350710e-01 6.753550e-02
[21,] 0.9349536 1.300928e-01 6.504638e-02
[22,] 0.9468360 1.063280e-01 5.316398e-02
[23,] 0.9720287 5.594268e-02 2.797134e-02
[24,] 0.9702621 5.947574e-02 2.973787e-02
[25,] 0.9631453 7.370948e-02 3.685474e-02
[26,] 0.9685837 6.283260e-02 3.141630e-02
[27,] 0.9756204 4.875927e-02 2.437963e-02
[28,] 0.9849371 3.012580e-02 1.506290e-02
[29,] 0.9788162 4.236758e-02 2.118379e-02
[30,] 0.9685318 6.293634e-02 3.146817e-02
[31,] 0.9588357 8.232865e-02 4.116433e-02
[32,] 0.9704716 5.905686e-02 2.952843e-02
[33,] 0.9984400 3.120052e-03 1.560026e-03
[34,] 0.9997901 4.198816e-04 2.099408e-04
[35,] 0.9999582 8.361809e-05 4.180905e-05
[36,] 0.9999993 1.426673e-06 7.133363e-07
[37,] 0.9999998 3.390404e-07 1.695202e-07
[38,] 0.9999996 7.913940e-07 3.956970e-07
[39,] 0.9999998 4.422845e-07 2.211422e-07
[40,] 0.9999988 2.308626e-06 1.154313e-06
[41,] 0.9999964 7.185183e-06 3.592591e-06
[42,] 0.9999818 3.648058e-05 1.824029e-05
[43,] 0.9999211 1.577595e-04 7.887975e-05
[44,] 0.9996401 7.198795e-04 3.599398e-04
[45,] 0.9998325 3.350529e-04 1.675265e-04
[46,] 0.9990626 1.874714e-03 9.373572e-04
[47,] 0.9953465 9.306946e-03 4.653473e-03
> postscript(file="/var/www/html/rcomp/tmp/13l4w1291029019.ps",horizontal=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/2wclz1291029019.ps",horizontal=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/3wclz1291029019.ps",horizontal=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/4wclz1291029019.ps",horizontal=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/5wclz1291029019.ps",horizontal=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 = 80
Frequency = 1
1 2 3 4 5 6
-40.29025670 6.43167013 -10.71882813 -31.57597098 9.43167013 -24.00454241
7 8 9 10 11 12
5.42402902 -28.69629212 -21.44362478 -0.88212613 -5.58399944 54.38970856
13 14 15 16 17 18
3.29771001 -29.82261112 -3.13086142 0.01199573 12.17738888 -3.41657570
19 20 21 22 23 24
61.01199573 3.89167459 39.14434194 -10.13640739 -20.99603272 -29.02232473
25 26 27 28 29 30
2.88567673 5.76535559 -2.54289470 18.59996244 -20.23464441 30.17139101
31 32 33 34 35 36
7.59996244 -1.52035870 -6.26769135 -14.54844067 22.43418195 -6.43435802
37 38 39 40 41 42
-6.52635656 -11.64667770 -13.95492799 0.18792915 2.35332230 34.75935772
43 44 45 46 47 48
-11.81207085 38.06760801 15.32027536 14.03952604 19.17990070 22.15360870
49 50 51 52 53 54
46.06161015 27.94128901 20.63303872 12.77589586 -4.05871099 -27.65267556
55 56 57 58 59 60
-39.22410414 5.49782269 -26.09175793 11.46974072 -18.23213259 -28.25842459
61 62 63 64 65 66
-5.50817517 7.37150369 7.06325340 7.20611054 7.37150369 -2.22246088
67 68 69 70 71 72
-8.79388946 -3.91421059 -0.66154325 0.05770743 3.19808209 -12.82820991
73 74 75 76 77 78
0.07979154 -6.04052960 2.65122011 -7.20592274 -7.04052960 -7.63449417
79 80
-14.20592274 -13.32624388
> postscript(file="/var/www/html/rcomp/tmp/674k11291029019.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 80
Frequency = 1
lag(myerror, k = 1) myerror
0 -40.29025670 NA
1 6.43167013 -40.29025670
2 -10.71882813 6.43167013
3 -31.57597098 -10.71882813
4 9.43167013 -31.57597098
5 -24.00454241 9.43167013
6 5.42402902 -24.00454241
7 -28.69629212 5.42402902
8 -21.44362478 -28.69629212
9 -0.88212613 -21.44362478
10 -5.58399944 -0.88212613
11 54.38970856 -5.58399944
12 3.29771001 54.38970856
13 -29.82261112 3.29771001
14 -3.13086142 -29.82261112
15 0.01199573 -3.13086142
16 12.17738888 0.01199573
17 -3.41657570 12.17738888
18 61.01199573 -3.41657570
19 3.89167459 61.01199573
20 39.14434194 3.89167459
21 -10.13640739 39.14434194
22 -20.99603272 -10.13640739
23 -29.02232473 -20.99603272
24 2.88567673 -29.02232473
25 5.76535559 2.88567673
26 -2.54289470 5.76535559
27 18.59996244 -2.54289470
28 -20.23464441 18.59996244
29 30.17139101 -20.23464441
30 7.59996244 30.17139101
31 -1.52035870 7.59996244
32 -6.26769135 -1.52035870
33 -14.54844067 -6.26769135
34 22.43418195 -14.54844067
35 -6.43435802 22.43418195
36 -6.52635656 -6.43435802
37 -11.64667770 -6.52635656
38 -13.95492799 -11.64667770
39 0.18792915 -13.95492799
40 2.35332230 0.18792915
41 34.75935772 2.35332230
42 -11.81207085 34.75935772
43 38.06760801 -11.81207085
44 15.32027536 38.06760801
45 14.03952604 15.32027536
46 19.17990070 14.03952604
47 22.15360870 19.17990070
48 46.06161015 22.15360870
49 27.94128901 46.06161015
50 20.63303872 27.94128901
51 12.77589586 20.63303872
52 -4.05871099 12.77589586
53 -27.65267556 -4.05871099
54 -39.22410414 -27.65267556
55 5.49782269 -39.22410414
56 -26.09175793 5.49782269
57 11.46974072 -26.09175793
58 -18.23213259 11.46974072
59 -28.25842459 -18.23213259
60 -5.50817517 -28.25842459
61 7.37150369 -5.50817517
62 7.06325340 7.37150369
63 7.20611054 7.06325340
64 7.37150369 7.20611054
65 -2.22246088 7.37150369
66 -8.79388946 -2.22246088
67 -3.91421059 -8.79388946
68 -0.66154325 -3.91421059
69 0.05770743 -0.66154325
70 3.19808209 0.05770743
71 -12.82820991 3.19808209
72 0.07979154 -12.82820991
73 -6.04052960 0.07979154
74 2.65122011 -6.04052960
75 -7.20592274 2.65122011
76 -7.04052960 -7.20592274
77 -7.63449417 -7.04052960
78 -14.20592274 -7.63449417
79 -13.32624388 -14.20592274
80 NA -13.32624388
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 6.43167013 -40.29025670
[2,] -10.71882813 6.43167013
[3,] -31.57597098 -10.71882813
[4,] 9.43167013 -31.57597098
[5,] -24.00454241 9.43167013
[6,] 5.42402902 -24.00454241
[7,] -28.69629212 5.42402902
[8,] -21.44362478 -28.69629212
[9,] -0.88212613 -21.44362478
[10,] -5.58399944 -0.88212613
[11,] 54.38970856 -5.58399944
[12,] 3.29771001 54.38970856
[13,] -29.82261112 3.29771001
[14,] -3.13086142 -29.82261112
[15,] 0.01199573 -3.13086142
[16,] 12.17738888 0.01199573
[17,] -3.41657570 12.17738888
[18,] 61.01199573 -3.41657570
[19,] 3.89167459 61.01199573
[20,] 39.14434194 3.89167459
[21,] -10.13640739 39.14434194
[22,] -20.99603272 -10.13640739
[23,] -29.02232473 -20.99603272
[24,] 2.88567673 -29.02232473
[25,] 5.76535559 2.88567673
[26,] -2.54289470 5.76535559
[27,] 18.59996244 -2.54289470
[28,] -20.23464441 18.59996244
[29,] 30.17139101 -20.23464441
[30,] 7.59996244 30.17139101
[31,] -1.52035870 7.59996244
[32,] -6.26769135 -1.52035870
[33,] -14.54844067 -6.26769135
[34,] 22.43418195 -14.54844067
[35,] -6.43435802 22.43418195
[36,] -6.52635656 -6.43435802
[37,] -11.64667770 -6.52635656
[38,] -13.95492799 -11.64667770
[39,] 0.18792915 -13.95492799
[40,] 2.35332230 0.18792915
[41,] 34.75935772 2.35332230
[42,] -11.81207085 34.75935772
[43,] 38.06760801 -11.81207085
[44,] 15.32027536 38.06760801
[45,] 14.03952604 15.32027536
[46,] 19.17990070 14.03952604
[47,] 22.15360870 19.17990070
[48,] 46.06161015 22.15360870
[49,] 27.94128901 46.06161015
[50,] 20.63303872 27.94128901
[51,] 12.77589586 20.63303872
[52,] -4.05871099 12.77589586
[53,] -27.65267556 -4.05871099
[54,] -39.22410414 -27.65267556
[55,] 5.49782269 -39.22410414
[56,] -26.09175793 5.49782269
[57,] 11.46974072 -26.09175793
[58,] -18.23213259 11.46974072
[59,] -28.25842459 -18.23213259
[60,] -5.50817517 -28.25842459
[61,] 7.37150369 -5.50817517
[62,] 7.06325340 7.37150369
[63,] 7.20611054 7.06325340
[64,] 7.37150369 7.20611054
[65,] -2.22246088 7.37150369
[66,] -8.79388946 -2.22246088
[67,] -3.91421059 -8.79388946
[68,] -0.66154325 -3.91421059
[69,] 0.05770743 -0.66154325
[70,] 3.19808209 0.05770743
[71,] -12.82820991 3.19808209
[72,] 0.07979154 -12.82820991
[73,] -6.04052960 0.07979154
[74,] 2.65122011 -6.04052960
[75,] -7.20592274 2.65122011
[76,] -7.04052960 -7.20592274
[77,] -7.63449417 -7.04052960
[78,] -14.20592274 -7.63449417
[79,] -13.32624388 -14.20592274
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 6.43167013 -40.29025670
2 -10.71882813 6.43167013
3 -31.57597098 -10.71882813
4 9.43167013 -31.57597098
5 -24.00454241 9.43167013
6 5.42402902 -24.00454241
7 -28.69629212 5.42402902
8 -21.44362478 -28.69629212
9 -0.88212613 -21.44362478
10 -5.58399944 -0.88212613
11 54.38970856 -5.58399944
12 3.29771001 54.38970856
13 -29.82261112 3.29771001
14 -3.13086142 -29.82261112
15 0.01199573 -3.13086142
16 12.17738888 0.01199573
17 -3.41657570 12.17738888
18 61.01199573 -3.41657570
19 3.89167459 61.01199573
20 39.14434194 3.89167459
21 -10.13640739 39.14434194
22 -20.99603272 -10.13640739
23 -29.02232473 -20.99603272
24 2.88567673 -29.02232473
25 5.76535559 2.88567673
26 -2.54289470 5.76535559
27 18.59996244 -2.54289470
28 -20.23464441 18.59996244
29 30.17139101 -20.23464441
30 7.59996244 30.17139101
31 -1.52035870 7.59996244
32 -6.26769135 -1.52035870
33 -14.54844067 -6.26769135
34 22.43418195 -14.54844067
35 -6.43435802 22.43418195
36 -6.52635656 -6.43435802
37 -11.64667770 -6.52635656
38 -13.95492799 -11.64667770
39 0.18792915 -13.95492799
40 2.35332230 0.18792915
41 34.75935772 2.35332230
42 -11.81207085 34.75935772
43 38.06760801 -11.81207085
44 15.32027536 38.06760801
45 14.03952604 15.32027536
46 19.17990070 14.03952604
47 22.15360870 19.17990070
48 46.06161015 22.15360870
49 27.94128901 46.06161015
50 20.63303872 27.94128901
51 12.77589586 20.63303872
52 -4.05871099 12.77589586
53 -27.65267556 -4.05871099
54 -39.22410414 -27.65267556
55 5.49782269 -39.22410414
56 -26.09175793 5.49782269
57 11.46974072 -26.09175793
58 -18.23213259 11.46974072
59 -28.25842459 -18.23213259
60 -5.50817517 -28.25842459
61 7.37150369 -5.50817517
62 7.06325340 7.37150369
63 7.20611054 7.06325340
64 7.37150369 7.20611054
65 -2.22246088 7.37150369
66 -8.79388946 -2.22246088
67 -3.91421059 -8.79388946
68 -0.66154325 -3.91421059
69 0.05770743 -0.66154325
70 3.19808209 0.05770743
71 -12.82820991 3.19808209
72 0.07979154 -12.82820991
73 -6.04052960 0.07979154
74 2.65122011 -6.04052960
75 -7.20592274 2.65122011
76 -7.04052960 -7.20592274
77 -7.63449417 -7.04052960
78 -14.20592274 -7.63449417
79 -13.32624388 -14.20592274
> 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/7hvjm1291029019.ps",horizontal=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/8hvjm1291029019.ps",horizontal=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/9am0p1291029019.ps",horizontal=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/10am0p1291029019.ps",horizontal=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/11v5hd1291029019.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/12zng11291029019.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/13vxva1291029019.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/14gxcy1291029019.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/152gsm1291029019.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/16g88u1291029019.tab")
+ }
>
> try(system("convert tmp/13l4w1291029019.ps tmp/13l4w1291029019.png",intern=TRUE))
character(0)
> try(system("convert tmp/2wclz1291029019.ps tmp/2wclz1291029019.png",intern=TRUE))
character(0)
> try(system("convert tmp/3wclz1291029019.ps tmp/3wclz1291029019.png",intern=TRUE))
character(0)
> try(system("convert tmp/4wclz1291029019.ps tmp/4wclz1291029019.png",intern=TRUE))
character(0)
> try(system("convert tmp/5wclz1291029019.ps tmp/5wclz1291029019.png",intern=TRUE))
character(0)
> try(system("convert tmp/674k11291029019.ps tmp/674k11291029019.png",intern=TRUE))
character(0)
> try(system("convert tmp/7hvjm1291029019.ps tmp/7hvjm1291029019.png",intern=TRUE))
character(0)
> try(system("convert tmp/8hvjm1291029019.ps tmp/8hvjm1291029019.png",intern=TRUE))
character(0)
> try(system("convert tmp/9am0p1291029019.ps tmp/9am0p1291029019.png",intern=TRUE))
character(0)
> try(system("convert tmp/10am0p1291029019.ps tmp/10am0p1291029019.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.756 1.691 13.075