R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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.
Natural language support but running in an English locale
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(9084
+ ,0
+ ,9700
+ ,9081
+ ,9743
+ ,0
+ ,9081
+ ,9084
+ ,8587
+ ,0
+ ,9084
+ ,9743
+ ,9731
+ ,0
+ ,9743
+ ,8587
+ ,9563
+ ,0
+ ,8587
+ ,9731
+ ,9998
+ ,0
+ ,9731
+ ,9563
+ ,9437
+ ,0
+ ,9563
+ ,9998
+ ,10038
+ ,0
+ ,9998
+ ,9437
+ ,9918
+ ,0
+ ,9437
+ ,10038
+ ,9252
+ ,0
+ ,10038
+ ,9918
+ ,9737
+ ,0
+ ,9918
+ ,9252
+ ,9035
+ ,0
+ ,9252
+ ,9737
+ ,9133
+ ,0
+ ,9737
+ ,9035
+ ,9487
+ ,0
+ ,9035
+ ,9133
+ ,8700
+ ,0
+ ,9133
+ ,9487
+ ,9627
+ ,0
+ ,9487
+ ,8700
+ ,8947
+ ,0
+ ,8700
+ ,9627
+ ,9283
+ ,0
+ ,9627
+ ,8947
+ ,8829
+ ,0
+ ,8947
+ ,9283
+ ,9947
+ ,0
+ ,9283
+ ,8829
+ ,9628
+ ,0
+ ,8829
+ ,9947
+ ,9318
+ ,0
+ ,9947
+ ,9628
+ ,9605
+ ,0
+ ,9628
+ ,9318
+ ,8640
+ ,0
+ ,9318
+ ,9605
+ ,9214
+ ,0
+ ,9605
+ ,8640
+ ,9567
+ ,0
+ ,8640
+ ,9214
+ ,8547
+ ,0
+ ,9214
+ ,9567
+ ,9185
+ ,0
+ ,9567
+ ,8547
+ ,9470
+ ,0
+ ,8547
+ ,9185
+ ,9123
+ ,0
+ ,9185
+ ,9470
+ ,9278
+ ,0
+ ,9470
+ ,9123
+ ,10170
+ ,0
+ ,9123
+ ,9278
+ ,9434
+ ,0
+ ,9278
+ ,10170
+ ,9655
+ ,0
+ ,10170
+ ,9434
+ ,9429
+ ,0
+ ,9434
+ ,9655
+ ,8739
+ ,0
+ ,9655
+ ,9429
+ ,9552
+ ,0
+ ,9429
+ ,8739
+ ,9687
+ ,0
+ ,8739
+ ,9552
+ ,9019
+ ,0
+ ,9552
+ ,9687
+ ,9672
+ ,0
+ ,9687
+ ,9019
+ ,9206
+ ,0
+ ,9019
+ ,9672
+ ,9069
+ ,0
+ ,9672
+ ,9206
+ ,9788
+ ,0
+ ,9206
+ ,9069
+ ,10312
+ ,0
+ ,9069
+ ,9788
+ ,10105
+ ,0
+ ,9788
+ ,10312
+ ,9863
+ ,0
+ ,10312
+ ,10105
+ ,9656
+ ,0
+ ,10105
+ ,9863
+ ,9295
+ ,0
+ ,9863
+ ,9656
+ ,9946
+ ,0
+ ,9656
+ ,9295
+ ,9701
+ ,0
+ ,9295
+ ,9946
+ ,9049
+ ,0
+ ,9946
+ ,9701
+ ,10190
+ ,0
+ ,9701
+ ,9049
+ ,9706
+ ,0
+ ,9049
+ ,10190
+ ,9765
+ ,0
+ ,10190
+ ,9706
+ ,9893
+ ,0
+ ,9706
+ ,9765
+ ,9994
+ ,0
+ ,9765
+ ,9893
+ ,10433
+ ,1
+ ,9893
+ ,9994
+ ,10073
+ ,1
+ ,9994
+ ,10433
+ ,10112
+ ,1
+ ,10433
+ ,10073
+ ,9266
+ ,1
+ ,10073
+ ,10112
+ ,9820
+ ,1
+ ,10112
+ ,9266
+ ,10097
+ ,1
+ ,9266
+ ,9820
+ ,9115
+ ,1
+ ,9820
+ ,10097
+ ,10411
+ ,1
+ ,10097
+ ,9115
+ ,9678
+ ,1
+ ,9115
+ ,10411
+ ,10408
+ ,1
+ ,10411
+ ,9678
+ ,10153
+ ,1
+ ,9678
+ ,10408
+ ,10368
+ ,1
+ ,10408
+ ,10153
+ ,10581
+ ,1
+ ,10153
+ ,10368
+ ,10597
+ ,1
+ ,10368
+ ,10581
+ ,10680
+ ,1
+ ,10581
+ ,10597
+ ,9738
+ ,1
+ ,10597
+ ,10680
+ ,9556
+ ,1
+ ,10680
+ ,9738)
+ ,dim=c(4
+ ,73)
+ ,dimnames=list(c('Births'
+ ,'x'
+ ,'y-1'
+ ,'y-2')
+ ,1:73))
> y <- array(NA,dim=c(4,73),dimnames=list(c('Births','x','y-1','y-2'),1:73))
> 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
Births x y-1 y-2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 9084 0 9700 9081 1 0 0 0 0 0 0 0 0 0 0 1
2 9743 0 9081 9084 0 1 0 0 0 0 0 0 0 0 0 2
3 8587 0 9084 9743 0 0 1 0 0 0 0 0 0 0 0 3
4 9731 0 9743 8587 0 0 0 1 0 0 0 0 0 0 0 4
5 9563 0 8587 9731 0 0 0 0 1 0 0 0 0 0 0 5
6 9998 0 9731 9563 0 0 0 0 0 1 0 0 0 0 0 6
7 9437 0 9563 9998 0 0 0 0 0 0 1 0 0 0 0 7
8 10038 0 9998 9437 0 0 0 0 0 0 0 1 0 0 0 8
9 9918 0 9437 10038 0 0 0 0 0 0 0 0 1 0 0 9
10 9252 0 10038 9918 0 0 0 0 0 0 0 0 0 1 0 10
11 9737 0 9918 9252 0 0 0 0 0 0 0 0 0 0 1 11
12 9035 0 9252 9737 0 0 0 0 0 0 0 0 0 0 0 12
13 9133 0 9737 9035 1 0 0 0 0 0 0 0 0 0 0 13
14 9487 0 9035 9133 0 1 0 0 0 0 0 0 0 0 0 14
15 8700 0 9133 9487 0 0 1 0 0 0 0 0 0 0 0 15
16 9627 0 9487 8700 0 0 0 1 0 0 0 0 0 0 0 16
17 8947 0 8700 9627 0 0 0 0 1 0 0 0 0 0 0 17
18 9283 0 9627 8947 0 0 0 0 0 1 0 0 0 0 0 18
19 8829 0 8947 9283 0 0 0 0 0 0 1 0 0 0 0 19
20 9947 0 9283 8829 0 0 0 0 0 0 0 1 0 0 0 20
21 9628 0 8829 9947 0 0 0 0 0 0 0 0 1 0 0 21
22 9318 0 9947 9628 0 0 0 0 0 0 0 0 0 1 0 22
23 9605 0 9628 9318 0 0 0 0 0 0 0 0 0 0 1 23
24 8640 0 9318 9605 0 0 0 0 0 0 0 0 0 0 0 24
25 9214 0 9605 8640 1 0 0 0 0 0 0 0 0 0 0 25
26 9567 0 8640 9214 0 1 0 0 0 0 0 0 0 0 0 26
27 8547 0 9214 9567 0 0 1 0 0 0 0 0 0 0 0 27
28 9185 0 9567 8547 0 0 0 1 0 0 0 0 0 0 0 28
29 9470 0 8547 9185 0 0 0 0 1 0 0 0 0 0 0 29
30 9123 0 9185 9470 0 0 0 0 0 1 0 0 0 0 0 30
31 9278 0 9470 9123 0 0 0 0 0 0 1 0 0 0 0 31
32 10170 0 9123 9278 0 0 0 0 0 0 0 1 0 0 0 32
33 9434 0 9278 10170 0 0 0 0 0 0 0 0 1 0 0 33
34 9655 0 10170 9434 0 0 0 0 0 0 0 0 0 1 0 34
35 9429 0 9434 9655 0 0 0 0 0 0 0 0 0 0 1 35
36 8739 0 9655 9429 0 0 0 0 0 0 0 0 0 0 0 36
37 9552 0 9429 8739 1 0 0 0 0 0 0 0 0 0 0 37
38 9687 0 8739 9552 0 1 0 0 0 0 0 0 0 0 0 38
39 9019 0 9552 9687 0 0 1 0 0 0 0 0 0 0 0 39
40 9672 0 9687 9019 0 0 0 1 0 0 0 0 0 0 0 40
41 9206 0 9019 9672 0 0 0 0 1 0 0 0 0 0 0 41
42 9069 0 9672 9206 0 0 0 0 0 1 0 0 0 0 0 42
43 9788 0 9206 9069 0 0 0 0 0 0 1 0 0 0 0 43
44 10312 0 9069 9788 0 0 0 0 0 0 0 1 0 0 0 44
45 10105 0 9788 10312 0 0 0 0 0 0 0 0 1 0 0 45
46 9863 0 10312 10105 0 0 0 0 0 0 0 0 0 1 0 46
47 9656 0 10105 9863 0 0 0 0 0 0 0 0 0 0 1 47
48 9295 0 9863 9656 0 0 0 0 0 0 0 0 0 0 0 48
49 9946 0 9656 9295 1 0 0 0 0 0 0 0 0 0 0 49
50 9701 0 9295 9946 0 1 0 0 0 0 0 0 0 0 0 50
51 9049 0 9946 9701 0 0 1 0 0 0 0 0 0 0 0 51
52 10190 0 9701 9049 0 0 0 1 0 0 0 0 0 0 0 52
53 9706 0 9049 10190 0 0 0 0 1 0 0 0 0 0 0 53
54 9765 0 10190 9706 0 0 0 0 0 1 0 0 0 0 0 54
55 9893 0 9706 9765 0 0 0 0 0 0 1 0 0 0 0 55
56 9994 0 9765 9893 0 0 0 0 0 0 0 1 0 0 0 56
57 10433 1 9893 9994 0 0 0 0 0 0 0 0 1 0 0 57
58 10073 1 9994 10433 0 0 0 0 0 0 0 0 0 1 0 58
59 10112 1 10433 10073 0 0 0 0 0 0 0 0 0 0 1 59
60 9266 1 10073 10112 0 0 0 0 0 0 0 0 0 0 0 60
61 9820 1 10112 9266 1 0 0 0 0 0 0 0 0 0 0 61
62 10097 1 9266 9820 0 1 0 0 0 0 0 0 0 0 0 62
63 9115 1 9820 10097 0 0 1 0 0 0 0 0 0 0 0 63
64 10411 1 10097 9115 0 0 0 1 0 0 0 0 0 0 0 64
65 9678 1 9115 10411 0 0 0 0 1 0 0 0 0 0 0 65
66 10408 1 10411 9678 0 0 0 0 0 1 0 0 0 0 0 66
67 10153 1 9678 10408 0 0 0 0 0 0 1 0 0 0 0 67
68 10368 1 10408 10153 0 0 0 0 0 0 0 1 0 0 0 68
69 10581 1 10153 10368 0 0 0 0 0 0 0 0 1 0 0 69
70 10597 1 10368 10581 0 0 0 0 0 0 0 0 0 1 0 70
71 10680 1 10581 10597 0 0 0 0 0 0 0 0 0 0 1 71
72 9738 1 10597 10680 0 0 0 0 0 0 0 0 0 0 0 72
73 9556 1 10680 9738 1 0 0 0 0 0 0 0 0 0 0 73
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x `y-1` `y-2` M1 M2
5632.4265 233.3273 0.1849 0.1419 485.5181 884.0771
M3 M4 M5 M6 M7 M8
-117.3717 921.6543 567.7858 616.7860 611.1999 1154.8356
M9 M10 M11 t
916.2465 598.9387 725.2722 4.7038
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-495.050 -179.186 9.123 162.336 564.486
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5632.4265 1521.6958 3.701 0.000485 ***
x 233.3273 122.0614 1.912 0.060965 .
`y-1` 0.1849 0.1302 1.421 0.160892
`y-2` 0.1419 0.1292 1.098 0.276861
M1 485.5181 176.8273 2.746 0.008064 **
M2 884.0771 177.4633 4.982 6.19e-06 ***
M3 -117.3717 157.7266 -0.744 0.459844
M4 921.6543 196.2524 4.696 1.71e-05 ***
M5 567.7858 191.7118 2.962 0.004455 **
M6 616.7860 162.7950 3.789 0.000367 ***
M7 611.1999 159.6183 3.829 0.000322 ***
M8 1154.8356 157.8431 7.316 9.45e-10 ***
M9 916.2465 162.9943 5.621 5.95e-07 ***
M10 598.9387 161.3986 3.711 0.000471 ***
M11 725.2722 158.0728 4.588 2.50e-05 ***
t 4.7038 2.4533 1.917 0.060212 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 267.7 on 57 degrees of freedom
Multiple R-squared: 0.7786, Adjusted R-squared: 0.7203
F-statistic: 13.36 on 15 and 57 DF, p-value: 1.549e-13
> 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.7224170 0.5551661 0.2775830
[2,] 0.6947573 0.6104854 0.3052427
[3,] 0.5973929 0.8052142 0.4026071
[4,] 0.5839222 0.8321555 0.4160778
[5,] 0.4765265 0.9530529 0.5234735
[6,] 0.3836829 0.7673659 0.6163171
[7,] 0.5099189 0.9801621 0.4900811
[8,] 0.4313908 0.8627816 0.5686092
[9,] 0.3360678 0.6721355 0.6639322
[10,] 0.3445186 0.6890371 0.6554814
[11,] 0.4978218 0.9956436 0.5021782
[12,] 0.4629612 0.9259223 0.5370388
[13,] 0.4071675 0.8143350 0.5928325
[14,] 0.5325597 0.9348807 0.4674403
[15,] 0.5316125 0.9367751 0.4683875
[16,] 0.5640350 0.8719299 0.4359650
[17,] 0.5321932 0.9356136 0.4678068
[18,] 0.4652824 0.9305648 0.5347176
[19,] 0.5746059 0.8507883 0.4253941
[20,] 0.5132674 0.9734651 0.4867326
[21,] 0.5520885 0.8958229 0.4479115
[22,] 0.4684357 0.9368715 0.5315643
[23,] 0.4069842 0.8139684 0.5930158
[24,] 0.6944650 0.6110699 0.3055350
[25,] 0.7444970 0.5110061 0.2555030
[26,] 0.7285284 0.5429431 0.2714716
[27,] 0.6612249 0.6775503 0.3387751
[28,] 0.5919436 0.8161128 0.4080564
[29,] 0.6349659 0.7300682 0.3650341
[30,] 0.5508373 0.8983253 0.4491627
[31,] 0.7859051 0.4281898 0.2140949
[32,] 0.6951630 0.6096739 0.3048370
[33,] 0.6051517 0.7896966 0.3948483
[34,] 0.5150269 0.9699462 0.4849731
[35,] 0.5820997 0.8358007 0.4179003
[36,] 0.4083645 0.8167290 0.5916355
> postscript(file="/var/www/html/freestat/rcomp/tmp/1vqs71291455188.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/freestat/rcomp/tmp/2vqs71291455188.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/freestat/rcomp/tmp/3vqs71291455188.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/freestat/rcomp/tmp/46zrr1291455188.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/freestat/rcomp/tmp/56zrr1291455188.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 = 73
Frequency = 1
1 2 3 4 5 6
-120.6193141 249.1613305 -4.1314204 138.2578860 370.9136307 564.4864259
7 8 9 10 11 12
-26.2710182 25.5278197 157.9010551 -289.6126178 181.0169953 253.9458321
13 14 15 16 17 18
-128.3823771 -61.7291945 79.6761096 9.1229876 -307.6761521 -100.3443562
19 20 21 22 23 24
-475.3761364 96.5517818 -63.2014374 -222.0923983 36.8367888 -190.9804837
25 26 27 28 29 30
-23.3853542 23.3799595 -156.0975040 -482.4133176 249.8715288 -309.2439233
31 32 33 34 35 36
-156.8415274 229.0004964 -428.3130688 44.7430089 -207.5391829 -185.7800068
37 38 39 40 41 42
276.6716731 20.6789081 179.9287542 -141.0064747 -226.9433031 -472.2989950
43 44 45 46 47 48
353.1928848 252.1938902 71.7851467 74.8521987 -190.5764792 243.1080276
49 50 51 52 53 54
493.3755879 -180.4770989 78.6358806 313.7027182 137.5814879 0.5354674
55 56 57 58 59 60
210.5522885 -265.8554090 135.7045980 7.3564174 -114.7952743 -179.1860444
61 62 63 64 65 66
-2.6101771 -51.0139048 -178.0118200 162.3362006 -223.7471922 316.8653812
67 68 69 70 71 72
94.7435086 -337.4185790 126.1237063 384.7533910 295.0571522 58.8926751
73
-495.0500385
> postscript(file="/var/www/html/freestat/rcomp/tmp/66zrr1291455188.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 -120.6193141 NA
1 249.1613305 -120.6193141
2 -4.1314204 249.1613305
3 138.2578860 -4.1314204
4 370.9136307 138.2578860
5 564.4864259 370.9136307
6 -26.2710182 564.4864259
7 25.5278197 -26.2710182
8 157.9010551 25.5278197
9 -289.6126178 157.9010551
10 181.0169953 -289.6126178
11 253.9458321 181.0169953
12 -128.3823771 253.9458321
13 -61.7291945 -128.3823771
14 79.6761096 -61.7291945
15 9.1229876 79.6761096
16 -307.6761521 9.1229876
17 -100.3443562 -307.6761521
18 -475.3761364 -100.3443562
19 96.5517818 -475.3761364
20 -63.2014374 96.5517818
21 -222.0923983 -63.2014374
22 36.8367888 -222.0923983
23 -190.9804837 36.8367888
24 -23.3853542 -190.9804837
25 23.3799595 -23.3853542
26 -156.0975040 23.3799595
27 -482.4133176 -156.0975040
28 249.8715288 -482.4133176
29 -309.2439233 249.8715288
30 -156.8415274 -309.2439233
31 229.0004964 -156.8415274
32 -428.3130688 229.0004964
33 44.7430089 -428.3130688
34 -207.5391829 44.7430089
35 -185.7800068 -207.5391829
36 276.6716731 -185.7800068
37 20.6789081 276.6716731
38 179.9287542 20.6789081
39 -141.0064747 179.9287542
40 -226.9433031 -141.0064747
41 -472.2989950 -226.9433031
42 353.1928848 -472.2989950
43 252.1938902 353.1928848
44 71.7851467 252.1938902
45 74.8521987 71.7851467
46 -190.5764792 74.8521987
47 243.1080276 -190.5764792
48 493.3755879 243.1080276
49 -180.4770989 493.3755879
50 78.6358806 -180.4770989
51 313.7027182 78.6358806
52 137.5814879 313.7027182
53 0.5354674 137.5814879
54 210.5522885 0.5354674
55 -265.8554090 210.5522885
56 135.7045980 -265.8554090
57 7.3564174 135.7045980
58 -114.7952743 7.3564174
59 -179.1860444 -114.7952743
60 -2.6101771 -179.1860444
61 -51.0139048 -2.6101771
62 -178.0118200 -51.0139048
63 162.3362006 -178.0118200
64 -223.7471922 162.3362006
65 316.8653812 -223.7471922
66 94.7435086 316.8653812
67 -337.4185790 94.7435086
68 126.1237063 -337.4185790
69 384.7533910 126.1237063
70 295.0571522 384.7533910
71 58.8926751 295.0571522
72 -495.0500385 58.8926751
73 NA -495.0500385
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 249.1613305 -120.6193141
[2,] -4.1314204 249.1613305
[3,] 138.2578860 -4.1314204
[4,] 370.9136307 138.2578860
[5,] 564.4864259 370.9136307
[6,] -26.2710182 564.4864259
[7,] 25.5278197 -26.2710182
[8,] 157.9010551 25.5278197
[9,] -289.6126178 157.9010551
[10,] 181.0169953 -289.6126178
[11,] 253.9458321 181.0169953
[12,] -128.3823771 253.9458321
[13,] -61.7291945 -128.3823771
[14,] 79.6761096 -61.7291945
[15,] 9.1229876 79.6761096
[16,] -307.6761521 9.1229876
[17,] -100.3443562 -307.6761521
[18,] -475.3761364 -100.3443562
[19,] 96.5517818 -475.3761364
[20,] -63.2014374 96.5517818
[21,] -222.0923983 -63.2014374
[22,] 36.8367888 -222.0923983
[23,] -190.9804837 36.8367888
[24,] -23.3853542 -190.9804837
[25,] 23.3799595 -23.3853542
[26,] -156.0975040 23.3799595
[27,] -482.4133176 -156.0975040
[28,] 249.8715288 -482.4133176
[29,] -309.2439233 249.8715288
[30,] -156.8415274 -309.2439233
[31,] 229.0004964 -156.8415274
[32,] -428.3130688 229.0004964
[33,] 44.7430089 -428.3130688
[34,] -207.5391829 44.7430089
[35,] -185.7800068 -207.5391829
[36,] 276.6716731 -185.7800068
[37,] 20.6789081 276.6716731
[38,] 179.9287542 20.6789081
[39,] -141.0064747 179.9287542
[40,] -226.9433031 -141.0064747
[41,] -472.2989950 -226.9433031
[42,] 353.1928848 -472.2989950
[43,] 252.1938902 353.1928848
[44,] 71.7851467 252.1938902
[45,] 74.8521987 71.7851467
[46,] -190.5764792 74.8521987
[47,] 243.1080276 -190.5764792
[48,] 493.3755879 243.1080276
[49,] -180.4770989 493.3755879
[50,] 78.6358806 -180.4770989
[51,] 313.7027182 78.6358806
[52,] 137.5814879 313.7027182
[53,] 0.5354674 137.5814879
[54,] 210.5522885 0.5354674
[55,] -265.8554090 210.5522885
[56,] 135.7045980 -265.8554090
[57,] 7.3564174 135.7045980
[58,] -114.7952743 7.3564174
[59,] -179.1860444 -114.7952743
[60,] -2.6101771 -179.1860444
[61,] -51.0139048 -2.6101771
[62,] -178.0118200 -51.0139048
[63,] 162.3362006 -178.0118200
[64,] -223.7471922 162.3362006
[65,] 316.8653812 -223.7471922
[66,] 94.7435086 316.8653812
[67,] -337.4185790 94.7435086
[68,] 126.1237063 -337.4185790
[69,] 384.7533910 126.1237063
[70,] 295.0571522 384.7533910
[71,] 58.8926751 295.0571522
[72,] -495.0500385 58.8926751
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 249.1613305 -120.6193141
2 -4.1314204 249.1613305
3 138.2578860 -4.1314204
4 370.9136307 138.2578860
5 564.4864259 370.9136307
6 -26.2710182 564.4864259
7 25.5278197 -26.2710182
8 157.9010551 25.5278197
9 -289.6126178 157.9010551
10 181.0169953 -289.6126178
11 253.9458321 181.0169953
12 -128.3823771 253.9458321
13 -61.7291945 -128.3823771
14 79.6761096 -61.7291945
15 9.1229876 79.6761096
16 -307.6761521 9.1229876
17 -100.3443562 -307.6761521
18 -475.3761364 -100.3443562
19 96.5517818 -475.3761364
20 -63.2014374 96.5517818
21 -222.0923983 -63.2014374
22 36.8367888 -222.0923983
23 -190.9804837 36.8367888
24 -23.3853542 -190.9804837
25 23.3799595 -23.3853542
26 -156.0975040 23.3799595
27 -482.4133176 -156.0975040
28 249.8715288 -482.4133176
29 -309.2439233 249.8715288
30 -156.8415274 -309.2439233
31 229.0004964 -156.8415274
32 -428.3130688 229.0004964
33 44.7430089 -428.3130688
34 -207.5391829 44.7430089
35 -185.7800068 -207.5391829
36 276.6716731 -185.7800068
37 20.6789081 276.6716731
38 179.9287542 20.6789081
39 -141.0064747 179.9287542
40 -226.9433031 -141.0064747
41 -472.2989950 -226.9433031
42 353.1928848 -472.2989950
43 252.1938902 353.1928848
44 71.7851467 252.1938902
45 74.8521987 71.7851467
46 -190.5764792 74.8521987
47 243.1080276 -190.5764792
48 493.3755879 243.1080276
49 -180.4770989 493.3755879
50 78.6358806 -180.4770989
51 313.7027182 78.6358806
52 137.5814879 313.7027182
53 0.5354674 137.5814879
54 210.5522885 0.5354674
55 -265.8554090 210.5522885
56 135.7045980 -265.8554090
57 7.3564174 135.7045980
58 -114.7952743 7.3564174
59 -179.1860444 -114.7952743
60 -2.6101771 -179.1860444
61 -51.0139048 -2.6101771
62 -178.0118200 -51.0139048
63 162.3362006 -178.0118200
64 -223.7471922 162.3362006
65 316.8653812 -223.7471922
66 94.7435086 316.8653812
67 -337.4185790 94.7435086
68 126.1237063 -337.4185790
69 384.7533910 126.1237063
70 295.0571522 384.7533910
71 58.8926751 295.0571522
72 -495.0500385 58.8926751
> 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/freestat/rcomp/tmp/7h88c1291455188.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/freestat/rcomp/tmp/82saa1291455189.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/freestat/rcomp/tmp/92saa1291455189.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/freestat/rcomp/tmp/10djav1291455189.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11gk8j1291455189.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/freestat/rcomp/tmp/12jkpp1291455189.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/freestat/rcomp/tmp/13xcmy1291455189.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/freestat/rcomp/tmp/14jc3m1291455189.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/freestat/rcomp/tmp/15ag9y1291455189.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/freestat/rcomp/tmp/168vix1291455189.tab")
+ }
>
> try(system("convert tmp/1vqs71291455188.ps tmp/1vqs71291455188.png",intern=TRUE))
character(0)
> try(system("convert tmp/2vqs71291455188.ps tmp/2vqs71291455188.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vqs71291455188.ps tmp/3vqs71291455188.png",intern=TRUE))
character(0)
> try(system("convert tmp/46zrr1291455188.ps tmp/46zrr1291455188.png",intern=TRUE))
character(0)
> try(system("convert tmp/56zrr1291455188.ps tmp/56zrr1291455188.png",intern=TRUE))
character(0)
> try(system("convert tmp/66zrr1291455188.ps tmp/66zrr1291455188.png",intern=TRUE))
character(0)
> try(system("convert tmp/7h88c1291455188.ps tmp/7h88c1291455188.png",intern=TRUE))
character(0)
> try(system("convert tmp/82saa1291455189.ps tmp/82saa1291455189.png",intern=TRUE))
character(0)
> try(system("convert tmp/92saa1291455189.ps tmp/92saa1291455189.png",intern=TRUE))
character(0)
> try(system("convert tmp/10djav1291455189.ps tmp/10djav1291455189.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.014 2.523 4.354