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(8587
+ ,0
+ ,9743
+ ,9084
+ ,9081
+ ,9700
+ ,9731
+ ,0
+ ,8587
+ ,9743
+ ,9084
+ ,9081
+ ,9563
+ ,0
+ ,9731
+ ,8587
+ ,9743
+ ,9084
+ ,9998
+ ,0
+ ,9563
+ ,9731
+ ,8587
+ ,9743
+ ,9437
+ ,0
+ ,9998
+ ,9563
+ ,9731
+ ,8587
+ ,10038
+ ,0
+ ,9437
+ ,9998
+ ,9563
+ ,9731
+ ,9918
+ ,0
+ ,10038
+ ,9437
+ ,9998
+ ,9563
+ ,9252
+ ,0
+ ,9918
+ ,10038
+ ,9437
+ ,9998
+ ,9737
+ ,0
+ ,9252
+ ,9918
+ ,10038
+ ,9437
+ ,9035
+ ,0
+ ,9737
+ ,9252
+ ,9918
+ ,10038
+ ,9133
+ ,0
+ ,9035
+ ,9737
+ ,9252
+ ,9918
+ ,9487
+ ,0
+ ,9133
+ ,9035
+ ,9737
+ ,9252
+ ,8700
+ ,0
+ ,9487
+ ,9133
+ ,9035
+ ,9737
+ ,9627
+ ,0
+ ,8700
+ ,9487
+ ,9133
+ ,9035
+ ,8947
+ ,0
+ ,9627
+ ,8700
+ ,9487
+ ,9133
+ ,9283
+ ,0
+ ,8947
+ ,9627
+ ,8700
+ ,9487
+ ,8829
+ ,0
+ ,9283
+ ,8947
+ ,9627
+ ,8700
+ ,9947
+ ,0
+ ,8829
+ ,9283
+ ,8947
+ ,9627
+ ,9628
+ ,0
+ ,9947
+ ,8829
+ ,9283
+ ,8947
+ ,9318
+ ,0
+ ,9628
+ ,9947
+ ,8829
+ ,9283
+ ,9605
+ ,0
+ ,9318
+ ,9628
+ ,9947
+ ,8829
+ ,8640
+ ,0
+ ,9605
+ ,9318
+ ,9628
+ ,9947
+ ,9214
+ ,0
+ ,8640
+ ,9605
+ ,9318
+ ,9628
+ ,9567
+ ,0
+ ,9214
+ ,8640
+ ,9605
+ ,9318
+ ,8547
+ ,0
+ ,9567
+ ,9214
+ ,8640
+ ,9605
+ ,9185
+ ,0
+ ,8547
+ ,9567
+ ,9214
+ ,8640
+ ,9470
+ ,0
+ ,9185
+ ,8547
+ ,9567
+ ,9214
+ ,9123
+ ,0
+ ,9470
+ ,9185
+ ,8547
+ ,9567
+ ,9278
+ ,0
+ ,9123
+ ,9470
+ ,9185
+ ,8547
+ ,10170
+ ,0
+ ,9278
+ ,9123
+ ,9470
+ ,9185
+ ,9434
+ ,0
+ ,10170
+ ,9278
+ ,9123
+ ,9470
+ ,9655
+ ,0
+ ,9434
+ ,10170
+ ,9278
+ ,9123
+ ,9429
+ ,0
+ ,9655
+ ,9434
+ ,10170
+ ,9278
+ ,8739
+ ,0
+ ,9429
+ ,9655
+ ,9434
+ ,10170
+ ,9552
+ ,0
+ ,8739
+ ,9429
+ ,9655
+ ,9434
+ ,9687
+ ,1
+ ,9552
+ ,8739
+ ,9429
+ ,9655
+ ,9019
+ ,1
+ ,9687
+ ,9552
+ ,8739
+ ,9429
+ ,9672
+ ,1
+ ,9019
+ ,9687
+ ,9552
+ ,8739
+ ,9206
+ ,1
+ ,9672
+ ,9019
+ ,9687
+ ,9552
+ ,9069
+ ,1
+ ,9206
+ ,9672
+ ,9019
+ ,9687
+ ,9788
+ ,1
+ ,9069
+ ,9206
+ ,9672
+ ,9019
+ ,10312
+ ,1
+ ,9788
+ ,9069
+ ,9206
+ ,9672
+ ,10105
+ ,1
+ ,10312
+ ,9788
+ ,9069
+ ,9206
+ ,9863
+ ,1
+ ,10105
+ ,10312
+ ,9788
+ ,9069
+ ,9656
+ ,1
+ ,9863
+ ,10105
+ ,10312
+ ,9788
+ ,9295
+ ,1
+ ,9656
+ ,9863
+ ,10105
+ ,10312
+ ,9946
+ ,1
+ ,9295
+ ,9656
+ ,9863
+ ,10105
+ ,9701
+ ,1
+ ,9946
+ ,9295
+ ,9656
+ ,9863
+ ,9049
+ ,1
+ ,9701
+ ,9946
+ ,9295
+ ,9656
+ ,10190
+ ,1
+ ,9049
+ ,9701
+ ,9946
+ ,9295
+ ,9706
+ ,1
+ ,10190
+ ,9049
+ ,9701
+ ,9946
+ ,9765
+ ,1
+ ,9706
+ ,10190
+ ,9049
+ ,9701
+ ,9893
+ ,1
+ ,9765
+ ,9706
+ ,10190
+ ,9049
+ ,9994
+ ,1
+ ,9893
+ ,9765
+ ,9706
+ ,10190
+ ,10433
+ ,1
+ ,9994
+ ,9893
+ ,9765
+ ,9706
+ ,10073
+ ,1
+ ,10433
+ ,9994
+ ,9893
+ ,9765
+ ,10112
+ ,1
+ ,10073
+ ,10433
+ ,9994
+ ,9893
+ ,9266
+ ,1
+ ,10112
+ ,10073
+ ,10433
+ ,9994
+ ,9820
+ ,1
+ ,9266
+ ,10112
+ ,10073
+ ,10433
+ ,10097
+ ,1
+ ,9820
+ ,9266
+ ,10112
+ ,10073
+ ,9115
+ ,1
+ ,10097
+ ,9820
+ ,9266
+ ,10112
+ ,10411
+ ,1
+ ,9115
+ ,10097
+ ,9820
+ ,9266
+ ,9678
+ ,1
+ ,10411
+ ,9115
+ ,10097
+ ,9820
+ ,10408
+ ,1
+ ,9678
+ ,10411
+ ,9115
+ ,10097
+ ,10153
+ ,1
+ ,10408
+ ,9678
+ ,10411
+ ,9115
+ ,10368
+ ,1
+ ,10153
+ ,10408
+ ,9678
+ ,10411
+ ,10581
+ ,1
+ ,10368
+ ,10153
+ ,10408
+ ,9678
+ ,10597
+ ,1
+ ,10581
+ ,10368
+ ,10153
+ ,10408
+ ,10680
+ ,1
+ ,10597
+ ,10581
+ ,10368
+ ,10153
+ ,9738
+ ,1
+ ,10680
+ ,10597
+ ,10581
+ ,10368
+ ,9556
+ ,1
+ ,9738
+ ,10680
+ ,10597
+ ,10581)
+ ,dim=c(6
+ ,71)
+ ,dimnames=list(c('births'
+ ,'difference'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3'
+ ,'Y4')
+ ,1:71))
> y <- array(NA,dim=c(6,71),dimnames=list(c('births','difference','Y1','Y2','Y3','Y4'),1:71))
> 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 difference Y1 Y2 Y3 Y4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 8587 0 9743 9084 9081 9700 1 0 0 0 0 0 0 0 0 0 0
2 9731 0 8587 9743 9084 9081 0 1 0 0 0 0 0 0 0 0 0
3 9563 0 9731 8587 9743 9084 0 0 1 0 0 0 0 0 0 0 0
4 9998 0 9563 9731 8587 9743 0 0 0 1 0 0 0 0 0 0 0
5 9437 0 9998 9563 9731 8587 0 0 0 0 1 0 0 0 0 0 0
6 10038 0 9437 9998 9563 9731 0 0 0 0 0 1 0 0 0 0 0
7 9918 0 10038 9437 9998 9563 0 0 0 0 0 0 1 0 0 0 0
8 9252 0 9918 10038 9437 9998 0 0 0 0 0 0 0 1 0 0 0
9 9737 0 9252 9918 10038 9437 0 0 0 0 0 0 0 0 1 0 0
10 9035 0 9737 9252 9918 10038 0 0 0 0 0 0 0 0 0 1 0
11 9133 0 9035 9737 9252 9918 0 0 0 0 0 0 0 0 0 0 1
12 9487 0 9133 9035 9737 9252 0 0 0 0 0 0 0 0 0 0 0
13 8700 0 9487 9133 9035 9737 1 0 0 0 0 0 0 0 0 0 0
14 9627 0 8700 9487 9133 9035 0 1 0 0 0 0 0 0 0 0 0
15 8947 0 9627 8700 9487 9133 0 0 1 0 0 0 0 0 0 0 0
16 9283 0 8947 9627 8700 9487 0 0 0 1 0 0 0 0 0 0 0
17 8829 0 9283 8947 9627 8700 0 0 0 0 1 0 0 0 0 0 0
18 9947 0 8829 9283 8947 9627 0 0 0 0 0 1 0 0 0 0 0
19 9628 0 9947 8829 9283 8947 0 0 0 0 0 0 1 0 0 0 0
20 9318 0 9628 9947 8829 9283 0 0 0 0 0 0 0 1 0 0 0
21 9605 0 9318 9628 9947 8829 0 0 0 0 0 0 0 0 1 0 0
22 8640 0 9605 9318 9628 9947 0 0 0 0 0 0 0 0 0 1 0
23 9214 0 8640 9605 9318 9628 0 0 0 0 0 0 0 0 0 0 1
24 9567 0 9214 8640 9605 9318 0 0 0 0 0 0 0 0 0 0 0
25 8547 0 9567 9214 8640 9605 1 0 0 0 0 0 0 0 0 0 0
26 9185 0 8547 9567 9214 8640 0 1 0 0 0 0 0 0 0 0 0
27 9470 0 9185 8547 9567 9214 0 0 1 0 0 0 0 0 0 0 0
28 9123 0 9470 9185 8547 9567 0 0 0 1 0 0 0 0 0 0 0
29 9278 0 9123 9470 9185 8547 0 0 0 0 1 0 0 0 0 0 0
30 10170 0 9278 9123 9470 9185 0 0 0 0 0 1 0 0 0 0 0
31 9434 0 10170 9278 9123 9470 0 0 0 0 0 0 1 0 0 0 0
32 9655 0 9434 10170 9278 9123 0 0 0 0 0 0 0 1 0 0 0
33 9429 0 9655 9434 10170 9278 0 0 0 0 0 0 0 0 1 0 0
34 8739 0 9429 9655 9434 10170 0 0 0 0 0 0 0 0 0 1 0
35 9552 0 8739 9429 9655 9434 0 0 0 0 0 0 0 0 0 0 1
36 9687 1 9552 8739 9429 9655 0 0 0 0 0 0 0 0 0 0 0
37 9019 1 9687 9552 8739 9429 1 0 0 0 0 0 0 0 0 0 0
38 9672 1 9019 9687 9552 8739 0 1 0 0 0 0 0 0 0 0 0
39 9206 1 9672 9019 9687 9552 0 0 1 0 0 0 0 0 0 0 0
40 9069 1 9206 9672 9019 9687 0 0 0 1 0 0 0 0 0 0 0
41 9788 1 9069 9206 9672 9019 0 0 0 0 1 0 0 0 0 0 0
42 10312 1 9788 9069 9206 9672 0 0 0 0 0 1 0 0 0 0 0
43 10105 1 10312 9788 9069 9206 0 0 0 0 0 0 1 0 0 0 0
44 9863 1 10105 10312 9788 9069 0 0 0 0 0 0 0 1 0 0 0
45 9656 1 9863 10105 10312 9788 0 0 0 0 0 0 0 0 1 0 0
46 9295 1 9656 9863 10105 10312 0 0 0 0 0 0 0 0 0 1 0
47 9946 1 9295 9656 9863 10105 0 0 0 0 0 0 0 0 0 0 1
48 9701 1 9946 9295 9656 9863 0 0 0 0 0 0 0 0 0 0 0
49 9049 1 9701 9946 9295 9656 1 0 0 0 0 0 0 0 0 0 0
50 10190 1 9049 9701 9946 9295 0 1 0 0 0 0 0 0 0 0 0
51 9706 1 10190 9049 9701 9946 0 0 1 0 0 0 0 0 0 0 0
52 9765 1 9706 10190 9049 9701 0 0 0 1 0 0 0 0 0 0 0
53 9893 1 9765 9706 10190 9049 0 0 0 0 1 0 0 0 0 0 0
54 9994 1 9893 9765 9706 10190 0 0 0 0 0 1 0 0 0 0 0
55 10433 1 9994 9893 9765 9706 0 0 0 0 0 0 1 0 0 0 0
56 10073 1 10433 9994 9893 9765 0 0 0 0 0 0 0 1 0 0 0
57 10112 1 10073 10433 9994 9893 0 0 0 0 0 0 0 0 1 0 0
58 9266 1 10112 10073 10433 9994 0 0 0 0 0 0 0 0 0 1 0
59 9820 1 9266 10112 10073 10433 0 0 0 0 0 0 0 0 0 0 1
60 10097 1 9820 9266 10112 10073 0 0 0 0 0 0 0 0 0 0 0
61 9115 1 10097 9820 9266 10112 1 0 0 0 0 0 0 0 0 0 0
62 10411 1 9115 10097 9820 9266 0 1 0 0 0 0 0 0 0 0 0
63 9678 1 10411 9115 10097 9820 0 0 1 0 0 0 0 0 0 0 0
64 10408 1 9678 10411 9115 10097 0 0 0 1 0 0 0 0 0 0 0
65 10153 1 10408 9678 10411 9115 0 0 0 0 1 0 0 0 0 0 0
66 10368 1 10153 10408 9678 10411 0 0 0 0 0 1 0 0 0 0 0
67 10581 1 10368 10153 10408 9678 0 0 0 0 0 0 1 0 0 0 0
68 10597 1 10581 10368 10153 10408 0 0 0 0 0 0 0 1 0 0 0
69 10680 1 10597 10581 10368 10153 0 0 0 0 0 0 0 0 1 0 0
70 9738 1 10680 10597 10581 10368 0 0 0 0 0 0 0 0 0 1 0
71 9556 1 9738 10680 10597 10581 0 0 0 0 0 0 0 0 0 0 1
t
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26
27 27
28 28
29 29
30 30
31 31
32 32
33 33
34 34
35 35
36 36
37 37
38 38
39 39
40 40
41 41
42 42
43 43
44 44
45 45
46 46
47 47
48 48
49 49
50 50
51 51
52 52
53 53
54 54
55 55
56 56
57 57
58 58
59 59
60 60
61 61
62 62
63 63
64 64
65 65
66 66
67 67
68 68
69 69
70 70
71 71
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) difference Y1 Y2 Y3 Y4
4290.89132 131.65152 0.10820 0.15434 0.23796 0.05107
M1 M2 M3 M4 M5 M6
-770.95853 176.62548 -253.79240 9.40421 -185.28757 403.52060
M7 M8 M9 M10 M11 t
199.73483 -101.13082 -119.14871 -847.61638 -304.04768 3.22630
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-621.77 -118.95 26.62 138.52 607.23
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4290.89132 1655.21505 2.592 0.012292 *
difference 131.65152 139.15939 0.946 0.348417
Y1 0.10820 0.14475 0.747 0.458083
Y2 0.15434 0.13694 1.127 0.264794
Y3 0.23796 0.13725 1.734 0.088762 .
Y4 0.05107 0.14601 0.350 0.727877
M1 -770.95853 209.51457 -3.680 0.000547 ***
M2 176.62548 236.03024 0.748 0.457577
M3 -253.79240 174.21628 -1.457 0.151080
M4 9.40421 240.16382 0.039 0.968912
M5 -185.28757 219.43780 -0.844 0.402256
M6 403.52060 191.71875 2.105 0.040071 *
M7 199.73483 210.27704 0.950 0.346492
M8 -101.13082 238.13587 -0.425 0.672791
M9 -119.14871 218.23283 -0.546 0.587377
M10 -847.61638 196.39695 -4.316 6.98e-05 ***
M11 -304.04768 222.05351 -1.369 0.176694
t 3.22630 3.55108 0.909 0.367706
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 271.5 on 53 degrees of freedom
Multiple R-squared: 0.7847, Adjusted R-squared: 0.7157
F-statistic: 11.36 on 17 and 53 DF, p-value: 3.732e-12
> 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.8221225 0.3557551 0.1778775
[2,] 0.7262391 0.5475218 0.2737609
[3,] 0.7340835 0.5318330 0.2659165
[4,] 0.6974267 0.6051467 0.3025733
[5,] 0.6212075 0.7575850 0.3787925
[6,] 0.5674769 0.8650462 0.4325231
[7,] 0.7776710 0.4446579 0.2223290
[8,] 0.7377039 0.5245923 0.2622961
[9,] 0.6740844 0.6518312 0.3259156
[10,] 0.8519401 0.2961198 0.1480599
[11,] 0.8132355 0.3735290 0.1867645
[12,] 0.8037901 0.3924198 0.1962099
[13,] 0.7292036 0.5415928 0.2707964
[14,] 0.7155454 0.5689092 0.2844546
[15,] 0.6986932 0.6026135 0.3013068
[16,] 0.6115776 0.7768448 0.3884224
[17,] 0.5582456 0.8835087 0.4417544
[18,] 0.4956086 0.9912171 0.5043914
[19,] 0.4376743 0.8753487 0.5623257
[20,] 0.6932033 0.6135934 0.3067967
[21,] 0.7398541 0.5202918 0.2601459
[22,] 0.7078132 0.5843736 0.2921868
[23,] 0.6741270 0.6517460 0.3258730
[24,] 0.6138300 0.7723400 0.3861700
[25,] 0.5357936 0.9284128 0.4642064
[26,] 0.4409962 0.8819924 0.5590038
[27,] 0.7147644 0.5704711 0.2852356
[28,] 0.5866265 0.8267471 0.4133735
[29,] 0.8226081 0.3547837 0.1773919
[30,] 0.6964647 0.6070706 0.3035353
> postscript(file="/var/www/html/freestat/rcomp/tmp/1rqbs1291922150.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/2rqbs1291922150.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/3kztd1291922150.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/4kztd1291922150.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/5kztd1291922150.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 = 71
Frequency = 1
1 2 3 4 5 6
-48.697573 198.755982 355.622363 607.232254 3.376573 -12.549489
7 8 9 10 11 12
-5.361978 -342.220656 133.789414 205.207941 -77.879277 -14.803346
13 14 15 16 17 18
54.778859 74.015372 -246.866371 -77.597124 -451.926468 185.768766
19 20 21 22 23 24
-28.788524 -88.316867 53.399343 -150.756249 26.622088 106.722542
25 26 27 28 29 30
-57.358384 -401.594289 285.681776 -232.357925 7.940514 244.289216
31 32 33 34 35 36
-343.569790 97.866998 -243.832946 -88.667951 272.075207 29.174218
37 38 39 40 41 42
164.553867 -240.038080 -320.045018 -621.772014 254.169534 207.025008
43 44 45 46 47 48
89.317443 -77.618352 -373.106874 73.377468 316.748984 -188.625742
49 50 51 52 53 54
-50.386958 111.686668 57.108888 -106.389455 43.182125 -413.909932
55 56 57 58 59 60
205.645654 46.725806 41.140268 -137.897324 18.067048 67.532327
61 62 63 64 65 66
-62.889810 257.174347 -131.501638 430.884264 143.257723 -210.623569
67 68 69 70 71
82.757195 363.563072 388.610795 98.736114 -555.634049
> postscript(file="/var/www/html/freestat/rcomp/tmp/6d8af1291922150.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 = 71
Frequency = 1
lag(myerror, k = 1) myerror
0 -48.697573 NA
1 198.755982 -48.697573
2 355.622363 198.755982
3 607.232254 355.622363
4 3.376573 607.232254
5 -12.549489 3.376573
6 -5.361978 -12.549489
7 -342.220656 -5.361978
8 133.789414 -342.220656
9 205.207941 133.789414
10 -77.879277 205.207941
11 -14.803346 -77.879277
12 54.778859 -14.803346
13 74.015372 54.778859
14 -246.866371 74.015372
15 -77.597124 -246.866371
16 -451.926468 -77.597124
17 185.768766 -451.926468
18 -28.788524 185.768766
19 -88.316867 -28.788524
20 53.399343 -88.316867
21 -150.756249 53.399343
22 26.622088 -150.756249
23 106.722542 26.622088
24 -57.358384 106.722542
25 -401.594289 -57.358384
26 285.681776 -401.594289
27 -232.357925 285.681776
28 7.940514 -232.357925
29 244.289216 7.940514
30 -343.569790 244.289216
31 97.866998 -343.569790
32 -243.832946 97.866998
33 -88.667951 -243.832946
34 272.075207 -88.667951
35 29.174218 272.075207
36 164.553867 29.174218
37 -240.038080 164.553867
38 -320.045018 -240.038080
39 -621.772014 -320.045018
40 254.169534 -621.772014
41 207.025008 254.169534
42 89.317443 207.025008
43 -77.618352 89.317443
44 -373.106874 -77.618352
45 73.377468 -373.106874
46 316.748984 73.377468
47 -188.625742 316.748984
48 -50.386958 -188.625742
49 111.686668 -50.386958
50 57.108888 111.686668
51 -106.389455 57.108888
52 43.182125 -106.389455
53 -413.909932 43.182125
54 205.645654 -413.909932
55 46.725806 205.645654
56 41.140268 46.725806
57 -137.897324 41.140268
58 18.067048 -137.897324
59 67.532327 18.067048
60 -62.889810 67.532327
61 257.174347 -62.889810
62 -131.501638 257.174347
63 430.884264 -131.501638
64 143.257723 430.884264
65 -210.623569 143.257723
66 82.757195 -210.623569
67 363.563072 82.757195
68 388.610795 363.563072
69 98.736114 388.610795
70 -555.634049 98.736114
71 NA -555.634049
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 198.755982 -48.697573
[2,] 355.622363 198.755982
[3,] 607.232254 355.622363
[4,] 3.376573 607.232254
[5,] -12.549489 3.376573
[6,] -5.361978 -12.549489
[7,] -342.220656 -5.361978
[8,] 133.789414 -342.220656
[9,] 205.207941 133.789414
[10,] -77.879277 205.207941
[11,] -14.803346 -77.879277
[12,] 54.778859 -14.803346
[13,] 74.015372 54.778859
[14,] -246.866371 74.015372
[15,] -77.597124 -246.866371
[16,] -451.926468 -77.597124
[17,] 185.768766 -451.926468
[18,] -28.788524 185.768766
[19,] -88.316867 -28.788524
[20,] 53.399343 -88.316867
[21,] -150.756249 53.399343
[22,] 26.622088 -150.756249
[23,] 106.722542 26.622088
[24,] -57.358384 106.722542
[25,] -401.594289 -57.358384
[26,] 285.681776 -401.594289
[27,] -232.357925 285.681776
[28,] 7.940514 -232.357925
[29,] 244.289216 7.940514
[30,] -343.569790 244.289216
[31,] 97.866998 -343.569790
[32,] -243.832946 97.866998
[33,] -88.667951 -243.832946
[34,] 272.075207 -88.667951
[35,] 29.174218 272.075207
[36,] 164.553867 29.174218
[37,] -240.038080 164.553867
[38,] -320.045018 -240.038080
[39,] -621.772014 -320.045018
[40,] 254.169534 -621.772014
[41,] 207.025008 254.169534
[42,] 89.317443 207.025008
[43,] -77.618352 89.317443
[44,] -373.106874 -77.618352
[45,] 73.377468 -373.106874
[46,] 316.748984 73.377468
[47,] -188.625742 316.748984
[48,] -50.386958 -188.625742
[49,] 111.686668 -50.386958
[50,] 57.108888 111.686668
[51,] -106.389455 57.108888
[52,] 43.182125 -106.389455
[53,] -413.909932 43.182125
[54,] 205.645654 -413.909932
[55,] 46.725806 205.645654
[56,] 41.140268 46.725806
[57,] -137.897324 41.140268
[58,] 18.067048 -137.897324
[59,] 67.532327 18.067048
[60,] -62.889810 67.532327
[61,] 257.174347 -62.889810
[62,] -131.501638 257.174347
[63,] 430.884264 -131.501638
[64,] 143.257723 430.884264
[65,] -210.623569 143.257723
[66,] 82.757195 -210.623569
[67,] 363.563072 82.757195
[68,] 388.610795 363.563072
[69,] 98.736114 388.610795
[70,] -555.634049 98.736114
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 198.755982 -48.697573
2 355.622363 198.755982
3 607.232254 355.622363
4 3.376573 607.232254
5 -12.549489 3.376573
6 -5.361978 -12.549489
7 -342.220656 -5.361978
8 133.789414 -342.220656
9 205.207941 133.789414
10 -77.879277 205.207941
11 -14.803346 -77.879277
12 54.778859 -14.803346
13 74.015372 54.778859
14 -246.866371 74.015372
15 -77.597124 -246.866371
16 -451.926468 -77.597124
17 185.768766 -451.926468
18 -28.788524 185.768766
19 -88.316867 -28.788524
20 53.399343 -88.316867
21 -150.756249 53.399343
22 26.622088 -150.756249
23 106.722542 26.622088
24 -57.358384 106.722542
25 -401.594289 -57.358384
26 285.681776 -401.594289
27 -232.357925 285.681776
28 7.940514 -232.357925
29 244.289216 7.940514
30 -343.569790 244.289216
31 97.866998 -343.569790
32 -243.832946 97.866998
33 -88.667951 -243.832946
34 272.075207 -88.667951
35 29.174218 272.075207
36 164.553867 29.174218
37 -240.038080 164.553867
38 -320.045018 -240.038080
39 -621.772014 -320.045018
40 254.169534 -621.772014
41 207.025008 254.169534
42 89.317443 207.025008
43 -77.618352 89.317443
44 -373.106874 -77.618352
45 73.377468 -373.106874
46 316.748984 73.377468
47 -188.625742 316.748984
48 -50.386958 -188.625742
49 111.686668 -50.386958
50 57.108888 111.686668
51 -106.389455 57.108888
52 43.182125 -106.389455
53 -413.909932 43.182125
54 205.645654 -413.909932
55 46.725806 205.645654
56 41.140268 46.725806
57 -137.897324 41.140268
58 18.067048 -137.897324
59 67.532327 18.067048
60 -62.889810 67.532327
61 257.174347 -62.889810
62 -131.501638 257.174347
63 430.884264 -131.501638
64 143.257723 430.884264
65 -210.623569 143.257723
66 82.757195 -210.623569
67 363.563072 82.757195
68 388.610795 363.563072
69 98.736114 388.610795
70 -555.634049 98.736114
> 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/7nz901291922150.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/8nz901291922150.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/9nz901291922150.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/10y9ql1291922150.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/11297r1291922150.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/12na5f1291922150.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/131jl61291922150.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/14522c1291922150.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/15q3ii1291922150.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/16t3z51291922150.tab")
+ }
>
> try(system("convert tmp/1rqbs1291922150.ps tmp/1rqbs1291922150.png",intern=TRUE))
character(0)
> try(system("convert tmp/2rqbs1291922150.ps tmp/2rqbs1291922150.png",intern=TRUE))
character(0)
> try(system("convert tmp/3kztd1291922150.ps tmp/3kztd1291922150.png",intern=TRUE))
character(0)
> try(system("convert tmp/4kztd1291922150.ps tmp/4kztd1291922150.png",intern=TRUE))
character(0)
> try(system("convert tmp/5kztd1291922150.ps tmp/5kztd1291922150.png",intern=TRUE))
character(0)
> try(system("convert tmp/6d8af1291922150.ps tmp/6d8af1291922150.png",intern=TRUE))
character(0)
> try(system("convert tmp/7nz901291922150.ps tmp/7nz901291922150.png",intern=TRUE))
character(0)
> try(system("convert tmp/8nz901291922150.ps tmp/8nz901291922150.png",intern=TRUE))
character(0)
> try(system("convert tmp/9nz901291922150.ps tmp/9nz901291922150.png",intern=TRUE))
character(0)
> try(system("convert tmp/10y9ql1291922150.ps tmp/10y9ql1291922150.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.017 2.539 4.410