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