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(9700,0,9081,0,9084,0,9743,0,8587,0,9731,0,9563,0,9998,0,9437,0,10038,0,9918,0,9252,0,9737,0,9035,0,9133,0,9487,0,8700,0,9627,0,8947,0,9283,0,8829,0,9947,0,9628,0,9318,0,9605,0,8640,0,9214,0,9567,0,8547,0,9185,0,9470,0,9123,0,9278,0,10170,0,9434,0,9655,0,9429,0,8739,0,9552,0,9687,1,9019,1,9672,1,9206,1,9069,1,9788,1,10312,1,10105,1,9863,1,9656,1,9295,1,9946,1,9701,1,9049,1,10190,1,9706,1,9765,1,9893,1,9994,1,10433,1,10073,1,10112,1,9266,1,9820,1,10097,1,9115,1,10411,1,9678,1,10408,1,10153,1,10368,1,10581,1,10597,1,10680,1,9738,1,9556,1),dim=c(2,75),dimnames=list(c('births','difference'),1:75))
> y <- array(NA,dim=c(2,75),dimnames=list(c('births','difference'),1:75))
> 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 9700 0 1 0 0 0 0 0 0 0 0 0 0 1
2 9081 0 0 1 0 0 0 0 0 0 0 0 0 2
3 9084 0 0 0 1 0 0 0 0 0 0 0 0 3
4 9743 0 0 0 0 1 0 0 0 0 0 0 0 4
5 8587 0 0 0 0 0 1 0 0 0 0 0 0 5
6 9731 0 0 0 0 0 0 1 0 0 0 0 0 6
7 9563 0 0 0 0 0 0 0 1 0 0 0 0 7
8 9998 0 0 0 0 0 0 0 0 1 0 0 0 8
9 9437 0 0 0 0 0 0 0 0 0 1 0 0 9
10 10038 0 0 0 0 0 0 0 0 0 0 1 0 10
11 9918 0 0 0 0 0 0 0 0 0 0 0 1 11
12 9252 0 0 0 0 0 0 0 0 0 0 0 0 12
13 9737 0 1 0 0 0 0 0 0 0 0 0 0 13
14 9035 0 0 1 0 0 0 0 0 0 0 0 0 14
15 9133 0 0 0 1 0 0 0 0 0 0 0 0 15
16 9487 0 0 0 0 1 0 0 0 0 0 0 0 16
17 8700 0 0 0 0 0 1 0 0 0 0 0 0 17
18 9627 0 0 0 0 0 0 1 0 0 0 0 0 18
19 8947 0 0 0 0 0 0 0 1 0 0 0 0 19
20 9283 0 0 0 0 0 0 0 0 1 0 0 0 20
21 8829 0 0 0 0 0 0 0 0 0 1 0 0 21
22 9947 0 0 0 0 0 0 0 0 0 0 1 0 22
23 9628 0 0 0 0 0 0 0 0 0 0 0 1 23
24 9318 0 0 0 0 0 0 0 0 0 0 0 0 24
25 9605 0 1 0 0 0 0 0 0 0 0 0 0 25
26 8640 0 0 1 0 0 0 0 0 0 0 0 0 26
27 9214 0 0 0 1 0 0 0 0 0 0 0 0 27
28 9567 0 0 0 0 1 0 0 0 0 0 0 0 28
29 8547 0 0 0 0 0 1 0 0 0 0 0 0 29
30 9185 0 0 0 0 0 0 1 0 0 0 0 0 30
31 9470 0 0 0 0 0 0 0 1 0 0 0 0 31
32 9123 0 0 0 0 0 0 0 0 1 0 0 0 32
33 9278 0 0 0 0 0 0 0 0 0 1 0 0 33
34 10170 0 0 0 0 0 0 0 0 0 0 1 0 34
35 9434 0 0 0 0 0 0 0 0 0 0 0 1 35
36 9655 0 0 0 0 0 0 0 0 0 0 0 0 36
37 9429 0 1 0 0 0 0 0 0 0 0 0 0 37
38 8739 0 0 1 0 0 0 0 0 0 0 0 0 38
39 9552 0 0 0 1 0 0 0 0 0 0 0 0 39
40 9687 1 0 0 0 1 0 0 0 0 0 0 0 40
41 9019 1 0 0 0 0 1 0 0 0 0 0 0 41
42 9672 1 0 0 0 0 0 1 0 0 0 0 0 42
43 9206 1 0 0 0 0 0 0 1 0 0 0 0 43
44 9069 1 0 0 0 0 0 0 0 1 0 0 0 44
45 9788 1 0 0 0 0 0 0 0 0 1 0 0 45
46 10312 1 0 0 0 0 0 0 0 0 0 1 0 46
47 10105 1 0 0 0 0 0 0 0 0 0 0 1 47
48 9863 1 0 0 0 0 0 0 0 0 0 0 0 48
49 9656 1 1 0 0 0 0 0 0 0 0 0 0 49
50 9295 1 0 1 0 0 0 0 0 0 0 0 0 50
51 9946 1 0 0 1 0 0 0 0 0 0 0 0 51
52 9701 1 0 0 0 1 0 0 0 0 0 0 0 52
53 9049 1 0 0 0 0 1 0 0 0 0 0 0 53
54 10190 1 0 0 0 0 0 1 0 0 0 0 0 54
55 9706 1 0 0 0 0 0 0 1 0 0 0 0 55
56 9765 1 0 0 0 0 0 0 0 1 0 0 0 56
57 9893 1 0 0 0 0 0 0 0 0 1 0 0 57
58 9994 1 0 0 0 0 0 0 0 0 0 1 0 58
59 10433 1 0 0 0 0 0 0 0 0 0 0 1 59
60 10073 1 0 0 0 0 0 0 0 0 0 0 0 60
61 10112 1 1 0 0 0 0 0 0 0 0 0 0 61
62 9266 1 0 1 0 0 0 0 0 0 0 0 0 62
63 9820 1 0 0 1 0 0 0 0 0 0 0 0 63
64 10097 1 0 0 0 1 0 0 0 0 0 0 0 64
65 9115 1 0 0 0 0 1 0 0 0 0 0 0 65
66 10411 1 0 0 0 0 0 1 0 0 0 0 0 66
67 9678 1 0 0 0 0 0 0 1 0 0 0 0 67
68 10408 1 0 0 0 0 0 0 0 1 0 0 0 68
69 10153 1 0 0 0 0 0 0 0 0 1 0 0 69
70 10368 1 0 0 0 0 0 0 0 0 0 1 0 70
71 10581 1 0 0 0 0 0 0 0 0 0 0 1 71
72 10597 1 0 0 0 0 0 0 0 0 0 0 0 72
73 10680 1 1 0 0 0 0 0 0 0 0 0 0 73
74 9738 1 0 1 0 0 0 0 0 0 0 0 0 74
75 9556 1 0 0 1 0 0 0 0 0 0 0 0 75
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) difference M1 M2 M3 M4
9430.311 286.707 99.162 -638.203 -284.711 -37.555
M5 M6 M7 M8 M9 M10
-920.277 41.000 -338.555 -164.444 -214.333 355.611
M11 t
228.722 5.222
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-713.354 -143.871 8.146 176.584 690.355
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9430.311 140.713 67.018 < 2e-16 ***
difference 286.707 134.583 2.130 0.037185 *
M1 99.162 158.574 0.625 0.534083
M2 -638.203 158.463 -4.027 0.000159 ***
M3 -284.711 158.413 -1.797 0.077243 .
M4 -37.555 166.198 -0.226 0.821983
M5 -920.277 165.759 -5.552 6.52e-07 ***
M6 41.000 165.378 0.248 0.805030
M7 -338.555 165.055 -2.051 0.044552 *
M8 -164.444 164.790 -0.998 0.322268
M9 -214.333 164.583 -1.302 0.197717
M10 355.611 164.436 2.163 0.034501 *
M11 228.722 164.347 1.392 0.169067
t 5.222 3.116 1.676 0.098874 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 284.6 on 61 degrees of freedom
Multiple R-squared: 0.7362, Adjusted R-squared: 0.6799
F-statistic: 13.09 on 13 and 61 DF, p-value: 3.895e-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.1287041 0.2574082 0.8712959
[2,] 0.0585261 0.1170522 0.9414739
[3,] 0.3323676 0.6647352 0.6676324
[4,] 0.5508841 0.8982318 0.4491159
[5,] 0.5867567 0.8264867 0.4132433
[6,] 0.5134490 0.9731020 0.4865510
[7,] 0.4078720 0.8157440 0.5921280
[8,] 0.3705938 0.7411876 0.6294062
[9,] 0.3190366 0.6380731 0.6809634
[10,] 0.2495087 0.4990174 0.7504913
[11,] 0.2804863 0.5609726 0.7195137
[12,] 0.2423723 0.4847447 0.7576277
[13,] 0.1824973 0.3649947 0.8175027
[14,] 0.1977471 0.3954943 0.8022529
[15,] 0.2999590 0.5999181 0.7000410
[16,] 0.2877079 0.5754157 0.7122921
[17,] 0.2887686 0.5775372 0.7112314
[18,] 0.3824245 0.7648491 0.6175755
[19,] 0.3742071 0.7484142 0.6257929
[20,] 0.4429157 0.8858313 0.5570843
[21,] 0.3848646 0.7697293 0.6151354
[22,] 0.3590140 0.7180281 0.6409860
[23,] 0.4463074 0.8926148 0.5536926
[24,] 0.3751816 0.7503633 0.6248184
[25,] 0.3643615 0.7287231 0.6356385
[26,] 0.3164595 0.6329191 0.6835405
[27,] 0.2823692 0.5647384 0.7176308
[28,] 0.5455372 0.9089257 0.4544628
[29,] 0.5628523 0.8742954 0.4371477
[30,] 0.6271859 0.7456281 0.3728141
[31,] 0.5644366 0.8711269 0.4355634
[32,] 0.5031675 0.9936651 0.4968325
[33,] 0.5594251 0.8811498 0.4405749
[34,] 0.4837902 0.9675804 0.5162098
[35,] 0.7830174 0.4339651 0.2169826
[36,] 0.7053027 0.5893947 0.2946973
[37,] 0.6481885 0.7036230 0.3518115
[38,] 0.5846574 0.8306852 0.4153426
[39,] 0.5803661 0.8392678 0.4196339
[40,] 0.5584801 0.8830399 0.4415199
[41,] 0.4293485 0.8586970 0.5706515
[42,] 0.2934597 0.5869193 0.7065403
> postscript(file="/var/www/html/rcomp/tmp/1gben1291917431.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/2gben1291917431.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/3q2vq1291917431.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/4q2vq1291917431.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/5q2vq1291917431.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 = 75
Frequency = 1
1 2 3 4 5 6
165.3047619 278.4476190 -77.2666667 329.3549603 50.8549603 228.3549603
7 8 9 10 11 12
434.6882937 690.3549603 174.0216270 199.8549603 201.5216270 -240.9783730
13 14 15 16 17 18
139.6375000 169.7803571 -90.9339286 10.6876984 101.1876984 61.6876984
19 20 21 22 23 24
-243.9789683 -87.3123016 -496.6456349 46.1876984 -151.1456349 -237.6456349
25 26 27 28 29 30
-55.0297619 -287.8869048 -72.6011905 28.0204365 -114.4795635 -442.9795635
31 32 33 34 35 36
216.3537698 -309.9795635 -110.3128968 206.5204365 -407.8128968 36.6871032
37 38 39 40 41 42
-293.6970238 -251.5541667 202.7315476 -201.3537698 8.1462302 -305.3537698
43 44 45 46 47 48
-397.0204365 -713.3537698 50.3128968 -0.8537698 -86.1871032 -104.6871032
49 50 51 52 53 54
-416.0712302 -44.9283730 247.3573413 -250.0210317 -24.5210317 149.9789683
55 56 57 58 59 60
40.3123016 -80.0210317 92.6456349 -381.5210317 179.1456349 42.6456349
61 62 63 64 65 66
-22.7384921 -136.5956349 58.6900794 83.3117063 -21.1882937 308.3117063
67 68 69 70 71 72
-50.3549603 500.3117063 289.9783730 -70.1882937 264.4783730 503.9783730
73 74 75
482.5942460 272.7371032 -267.9771825
> postscript(file="/var/www/html/rcomp/tmp/61uut1291917431.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 = 75
Frequency = 1
lag(myerror, k = 1) myerror
0 165.3047619 NA
1 278.4476190 165.3047619
2 -77.2666667 278.4476190
3 329.3549603 -77.2666667
4 50.8549603 329.3549603
5 228.3549603 50.8549603
6 434.6882937 228.3549603
7 690.3549603 434.6882937
8 174.0216270 690.3549603
9 199.8549603 174.0216270
10 201.5216270 199.8549603
11 -240.9783730 201.5216270
12 139.6375000 -240.9783730
13 169.7803571 139.6375000
14 -90.9339286 169.7803571
15 10.6876984 -90.9339286
16 101.1876984 10.6876984
17 61.6876984 101.1876984
18 -243.9789683 61.6876984
19 -87.3123016 -243.9789683
20 -496.6456349 -87.3123016
21 46.1876984 -496.6456349
22 -151.1456349 46.1876984
23 -237.6456349 -151.1456349
24 -55.0297619 -237.6456349
25 -287.8869048 -55.0297619
26 -72.6011905 -287.8869048
27 28.0204365 -72.6011905
28 -114.4795635 28.0204365
29 -442.9795635 -114.4795635
30 216.3537698 -442.9795635
31 -309.9795635 216.3537698
32 -110.3128968 -309.9795635
33 206.5204365 -110.3128968
34 -407.8128968 206.5204365
35 36.6871032 -407.8128968
36 -293.6970238 36.6871032
37 -251.5541667 -293.6970238
38 202.7315476 -251.5541667
39 -201.3537698 202.7315476
40 8.1462302 -201.3537698
41 -305.3537698 8.1462302
42 -397.0204365 -305.3537698
43 -713.3537698 -397.0204365
44 50.3128968 -713.3537698
45 -0.8537698 50.3128968
46 -86.1871032 -0.8537698
47 -104.6871032 -86.1871032
48 -416.0712302 -104.6871032
49 -44.9283730 -416.0712302
50 247.3573413 -44.9283730
51 -250.0210317 247.3573413
52 -24.5210317 -250.0210317
53 149.9789683 -24.5210317
54 40.3123016 149.9789683
55 -80.0210317 40.3123016
56 92.6456349 -80.0210317
57 -381.5210317 92.6456349
58 179.1456349 -381.5210317
59 42.6456349 179.1456349
60 -22.7384921 42.6456349
61 -136.5956349 -22.7384921
62 58.6900794 -136.5956349
63 83.3117063 58.6900794
64 -21.1882937 83.3117063
65 308.3117063 -21.1882937
66 -50.3549603 308.3117063
67 500.3117063 -50.3549603
68 289.9783730 500.3117063
69 -70.1882937 289.9783730
70 264.4783730 -70.1882937
71 503.9783730 264.4783730
72 482.5942460 503.9783730
73 272.7371032 482.5942460
74 -267.9771825 272.7371032
75 NA -267.9771825
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 278.4476190 165.3047619
[2,] -77.2666667 278.4476190
[3,] 329.3549603 -77.2666667
[4,] 50.8549603 329.3549603
[5,] 228.3549603 50.8549603
[6,] 434.6882937 228.3549603
[7,] 690.3549603 434.6882937
[8,] 174.0216270 690.3549603
[9,] 199.8549603 174.0216270
[10,] 201.5216270 199.8549603
[11,] -240.9783730 201.5216270
[12,] 139.6375000 -240.9783730
[13,] 169.7803571 139.6375000
[14,] -90.9339286 169.7803571
[15,] 10.6876984 -90.9339286
[16,] 101.1876984 10.6876984
[17,] 61.6876984 101.1876984
[18,] -243.9789683 61.6876984
[19,] -87.3123016 -243.9789683
[20,] -496.6456349 -87.3123016
[21,] 46.1876984 -496.6456349
[22,] -151.1456349 46.1876984
[23,] -237.6456349 -151.1456349
[24,] -55.0297619 -237.6456349
[25,] -287.8869048 -55.0297619
[26,] -72.6011905 -287.8869048
[27,] 28.0204365 -72.6011905
[28,] -114.4795635 28.0204365
[29,] -442.9795635 -114.4795635
[30,] 216.3537698 -442.9795635
[31,] -309.9795635 216.3537698
[32,] -110.3128968 -309.9795635
[33,] 206.5204365 -110.3128968
[34,] -407.8128968 206.5204365
[35,] 36.6871032 -407.8128968
[36,] -293.6970238 36.6871032
[37,] -251.5541667 -293.6970238
[38,] 202.7315476 -251.5541667
[39,] -201.3537698 202.7315476
[40,] 8.1462302 -201.3537698
[41,] -305.3537698 8.1462302
[42,] -397.0204365 -305.3537698
[43,] -713.3537698 -397.0204365
[44,] 50.3128968 -713.3537698
[45,] -0.8537698 50.3128968
[46,] -86.1871032 -0.8537698
[47,] -104.6871032 -86.1871032
[48,] -416.0712302 -104.6871032
[49,] -44.9283730 -416.0712302
[50,] 247.3573413 -44.9283730
[51,] -250.0210317 247.3573413
[52,] -24.5210317 -250.0210317
[53,] 149.9789683 -24.5210317
[54,] 40.3123016 149.9789683
[55,] -80.0210317 40.3123016
[56,] 92.6456349 -80.0210317
[57,] -381.5210317 92.6456349
[58,] 179.1456349 -381.5210317
[59,] 42.6456349 179.1456349
[60,] -22.7384921 42.6456349
[61,] -136.5956349 -22.7384921
[62,] 58.6900794 -136.5956349
[63,] 83.3117063 58.6900794
[64,] -21.1882937 83.3117063
[65,] 308.3117063 -21.1882937
[66,] -50.3549603 308.3117063
[67,] 500.3117063 -50.3549603
[68,] 289.9783730 500.3117063
[69,] -70.1882937 289.9783730
[70,] 264.4783730 -70.1882937
[71,] 503.9783730 264.4783730
[72,] 482.5942460 503.9783730
[73,] 272.7371032 482.5942460
[74,] -267.9771825 272.7371032
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 278.4476190 165.3047619
2 -77.2666667 278.4476190
3 329.3549603 -77.2666667
4 50.8549603 329.3549603
5 228.3549603 50.8549603
6 434.6882937 228.3549603
7 690.3549603 434.6882937
8 174.0216270 690.3549603
9 199.8549603 174.0216270
10 201.5216270 199.8549603
11 -240.9783730 201.5216270
12 139.6375000 -240.9783730
13 169.7803571 139.6375000
14 -90.9339286 169.7803571
15 10.6876984 -90.9339286
16 101.1876984 10.6876984
17 61.6876984 101.1876984
18 -243.9789683 61.6876984
19 -87.3123016 -243.9789683
20 -496.6456349 -87.3123016
21 46.1876984 -496.6456349
22 -151.1456349 46.1876984
23 -237.6456349 -151.1456349
24 -55.0297619 -237.6456349
25 -287.8869048 -55.0297619
26 -72.6011905 -287.8869048
27 28.0204365 -72.6011905
28 -114.4795635 28.0204365
29 -442.9795635 -114.4795635
30 216.3537698 -442.9795635
31 -309.9795635 216.3537698
32 -110.3128968 -309.9795635
33 206.5204365 -110.3128968
34 -407.8128968 206.5204365
35 36.6871032 -407.8128968
36 -293.6970238 36.6871032
37 -251.5541667 -293.6970238
38 202.7315476 -251.5541667
39 -201.3537698 202.7315476
40 8.1462302 -201.3537698
41 -305.3537698 8.1462302
42 -397.0204365 -305.3537698
43 -713.3537698 -397.0204365
44 50.3128968 -713.3537698
45 -0.8537698 50.3128968
46 -86.1871032 -0.8537698
47 -104.6871032 -86.1871032
48 -416.0712302 -104.6871032
49 -44.9283730 -416.0712302
50 247.3573413 -44.9283730
51 -250.0210317 247.3573413
52 -24.5210317 -250.0210317
53 149.9789683 -24.5210317
54 40.3123016 149.9789683
55 -80.0210317 40.3123016
56 92.6456349 -80.0210317
57 -381.5210317 92.6456349
58 179.1456349 -381.5210317
59 42.6456349 179.1456349
60 -22.7384921 42.6456349
61 -136.5956349 -22.7384921
62 58.6900794 -136.5956349
63 83.3117063 58.6900794
64 -21.1882937 83.3117063
65 308.3117063 -21.1882937
66 -50.3549603 308.3117063
67 500.3117063 -50.3549603
68 289.9783730 500.3117063
69 -70.1882937 289.9783730
70 264.4783730 -70.1882937
71 503.9783730 264.4783730
72 482.5942460 503.9783730
73 272.7371032 482.5942460
74 -267.9771825 272.7371032
> 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/7clue1291917431.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/8clue1291917431.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/9clue1291917431.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/10ncbh1291917431.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/11i5c01291917432.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/124na61291917432.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/13ix8w1291917432.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/14lfpk1291917432.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/156g581291917432.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/16sy4w1291917432.tab")
+ }
>
> try(system("convert tmp/1gben1291917431.ps tmp/1gben1291917431.png",intern=TRUE))
character(0)
> try(system("convert tmp/2gben1291917431.ps tmp/2gben1291917431.png",intern=TRUE))
character(0)
> try(system("convert tmp/3q2vq1291917431.ps tmp/3q2vq1291917431.png",intern=TRUE))
character(0)
> try(system("convert tmp/4q2vq1291917431.ps tmp/4q2vq1291917431.png",intern=TRUE))
character(0)
> try(system("convert tmp/5q2vq1291917431.ps tmp/5q2vq1291917431.png",intern=TRUE))
character(0)
> try(system("convert tmp/61uut1291917431.ps tmp/61uut1291917431.png",intern=TRUE))
character(0)
> try(system("convert tmp/7clue1291917431.ps tmp/7clue1291917431.png",intern=TRUE))
character(0)
> try(system("convert tmp/8clue1291917431.ps tmp/8clue1291917431.png",intern=TRUE))
character(0)
> try(system("convert tmp/9clue1291917431.ps tmp/9clue1291917431.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ncbh1291917431.ps tmp/10ncbh1291917431.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.605 1.748 6.792