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(9939 + ,2462 + ,9321 + ,9769 + ,9336 + ,3695 + ,9939 + ,9321 + ,10195 + ,4831 + ,9336 + ,9939 + ,9464 + ,5134 + ,10195 + ,9336 + ,10010 + ,6250 + ,9464 + ,10195 + ,10213 + ,5760 + ,10010 + ,9464 + ,9563 + ,6249 + ,10213 + ,10010 + ,9890 + ,2917 + ,9563 + ,10213 + ,9305 + ,1741 + ,9890 + ,9563 + ,9391 + ,2359 + ,9305 + ,9890 + ,9928 + ,1511 + ,9391 + ,9305 + ,8686 + ,2059 + ,9928 + ,9391 + ,9843 + ,2635 + ,8686 + ,9928 + ,9627 + ,2867 + ,9843 + ,8686 + ,10074 + ,4403 + ,9627 + ,9843 + ,9503 + ,5720 + ,10074 + ,9627 + ,10119 + ,4502 + ,9503 + ,10074 + ,10000 + ,5749 + ,10119 + ,9503 + ,9313 + ,5627 + ,10000 + ,10119 + ,9866 + ,2846 + ,9313 + ,10000 + ,9172 + ,1762 + ,9866 + ,9313 + ,9241 + ,2429 + ,9172 + ,9866 + ,9659 + ,1169 + ,9241 + ,9172 + ,8904 + ,2154 + ,9659 + ,9241 + ,9755 + ,2249 + ,8904 + ,9659 + ,9080 + ,2687 + ,9755 + ,8904 + ,9435 + ,4359 + ,9080 + ,9755 + ,8971 + ,5382 + ,9435 + ,9080 + ,10063 + ,4459 + ,8971 + ,9435 + ,9793 + ,6398 + ,10063 + ,8971 + ,9454 + ,4596 + ,9793 + ,10063 + ,9759 + ,3024 + ,9454 + ,9793 + ,8820 + ,1887 + ,9759 + ,9454 + ,9403 + ,2070 + ,8820 + ,9759 + ,9676 + ,1351 + ,9403 + ,8820 + ,8642 + ,2218 + ,9676 + ,9403 + ,9402 + ,2461 + ,8642 + ,9676 + ,9610 + ,3028 + ,9402 + ,8642 + ,9294 + ,4784 + ,9610 + ,9402 + ,9448 + ,4975 + ,9294 + ,9610 + ,10319 + ,4607 + ,9448 + ,9294 + ,9548 + ,6249 + ,10319 + ,9448 + ,9801 + ,4809 + ,9548 + ,10319 + ,9596 + ,3157 + ,9801 + ,9548 + ,8923 + ,1910 + ,9596 + ,9801 + ,9746 + ,2228 + ,8923 + ,9596 + ,9829 + ,1594 + ,9746 + ,8923 + ,9125 + ,2467 + ,9829 + ,9746 + ,9782 + ,2222 + ,9125 + ,9829 + ,9441 + ,3607 + ,9782 + ,9125 + ,9162 + ,4685 + ,9441 + ,9782 + ,9915 + ,4962 + ,9162 + ,9441 + ,10444 + ,5770 + ,9915 + ,9162 + ,10209 + ,5480 + ,10444 + ,9915 + ,9985 + ,5000 + ,10209 + ,10444 + ,9842 + ,3228 + ,9985 + ,10209 + ,9429 + ,1993 + ,9842 + ,9985 + ,10132 + ,2288 + ,9429 + ,9842 + ,9849 + ,1580 + ,10132 + ,9429 + ,9172 + ,2111 + ,9849 + ,10132 + ,10313 + ,2192 + ,9172 + ,9849 + ,9819 + ,3601 + ,10313 + ,9172 + ,9955 + ,4665 + ,9819 + ,10313 + ,10048 + ,4876 + ,9955 + ,9819 + ,10082 + ,5813 + ,10048 + ,9955 + ,10541 + ,5589 + ,10082 + ,10048 + ,10208 + ,5331 + ,10541 + ,10082 + ,10233 + ,3075 + ,10208 + ,10541 + ,9439 + ,2002 + ,10233 + ,10208 + ,9963 + ,2306 + ,9439 + ,10233 + ,10158 + ,1507 + ,9963 + ,9439 + ,9225 + ,1992 + ,10158 + ,9963 + ,10474 + ,2487 + ,9225 + ,10158 + ,9757 + ,3490 + ,10474 + ,9225 + ,10490 + ,4647 + ,9757 + ,10474 + ,10281 + ,5594 + ,10490 + ,9757 + ,10444 + ,5611 + ,10281 + ,10490 + ,10640 + ,5788 + ,10444 + ,10281 + ,10695 + ,6204 + ,10640 + ,10444 + ,10786 + ,3013 + ,10695 + ,10640 + ,9832 + ,1931 + ,10786 + ,10695 + ,9747 + ,2549 + ,9832 + ,10786 + ,10411 + ,1504 + ,9747 + ,9832 + ,9511 + ,2090 + ,10411 + ,9747 + ,10402 + ,2702 + ,9511 + ,10411 + ,9701 + ,2939 + ,10402 + ,9511 + ,10540 + ,4500 + ,9701 + ,10402 + ,10112 + ,6208 + ,10540 + ,9701 + ,10915 + ,6415 + ,10112 + ,10540 + ,11183 + ,5657 + ,10915 + ,10112 + ,10384 + ,5964 + ,11183 + ,10915 + ,10834 + ,3163 + ,10384 + ,11183 + ,9886 + ,1997 + ,10834 + ,10384 + ,10216 + ,2422 + ,9886 + ,10834) + ,dim=c(4 + ,94) + ,dimnames=list(c('geboortes' + ,'huwelijken' + ,'geboortes-1' + ,'geboortes-2') + ,1:94)) > y <- array(NA,dim=c(4,94),dimnames=list(c('geboortes','huwelijken','geboortes-1','geboortes-2'),1:94)) > 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 geboortes huwelijken geboortes-1 geboortes-2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 1 9939 2462 9321 9769 1 0 0 0 0 0 0 0 0 0 2 9336 3695 9939 9321 0 1 0 0 0 0 0 0 0 0 3 10195 4831 9336 9939 0 0 1 0 0 0 0 0 0 0 4 9464 5134 10195 9336 0 0 0 1 0 0 0 0 0 0 5 10010 6250 9464 10195 0 0 0 0 1 0 0 0 0 0 6 10213 5760 10010 9464 0 0 0 0 0 1 0 0 0 0 7 9563 6249 10213 10010 0 0 0 0 0 0 1 0 0 0 8 9890 2917 9563 10213 0 0 0 0 0 0 0 1 0 0 9 9305 1741 9890 9563 0 0 0 0 0 0 0 0 1 0 10 9391 2359 9305 9890 0 0 0 0 0 0 0 0 0 1 11 9928 1511 9391 9305 0 0 0 0 0 0 0 0 0 0 12 8686 2059 9928 9391 0 0 0 0 0 0 0 0 0 0 13 9843 2635 8686 9928 1 0 0 0 0 0 0 0 0 0 14 9627 2867 9843 8686 0 1 0 0 0 0 0 0 0 0 15 10074 4403 9627 9843 0 0 1 0 0 0 0 0 0 0 16 9503 5720 10074 9627 0 0 0 1 0 0 0 0 0 0 17 10119 4502 9503 10074 0 0 0 0 1 0 0 0 0 0 18 10000 5749 10119 9503 0 0 0 0 0 1 0 0 0 0 19 9313 5627 10000 10119 0 0 0 0 0 0 1 0 0 0 20 9866 2846 9313 10000 0 0 0 0 0 0 0 1 0 0 21 9172 1762 9866 9313 0 0 0 0 0 0 0 0 1 0 22 9241 2429 9172 9866 0 0 0 0 0 0 0 0 0 1 23 9659 1169 9241 9172 0 0 0 0 0 0 0 0 0 0 24 8904 2154 9659 9241 0 0 0 0 0 0 0 0 0 0 25 9755 2249 8904 9659 1 0 0 0 0 0 0 0 0 0 26 9080 2687 9755 8904 0 1 0 0 0 0 0 0 0 0 27 9435 4359 9080 9755 0 0 1 0 0 0 0 0 0 0 28 8971 5382 9435 9080 0 0 0 1 0 0 0 0 0 0 29 10063 4459 8971 9435 0 0 0 0 1 0 0 0 0 0 30 9793 6398 10063 8971 0 0 0 0 0 1 0 0 0 0 31 9454 4596 9793 10063 0 0 0 0 0 0 1 0 0 0 32 9759 3024 9454 9793 0 0 0 0 0 0 0 1 0 0 33 8820 1887 9759 9454 0 0 0 0 0 0 0 0 1 0 34 9403 2070 8820 9759 0 0 0 0 0 0 0 0 0 1 35 9676 1351 9403 8820 0 0 0 0 0 0 0 0 0 0 36 8642 2218 9676 9403 0 0 0 0 0 0 0 0 0 0 37 9402 2461 8642 9676 1 0 0 0 0 0 0 0 0 0 38 9610 3028 9402 8642 0 1 0 0 0 0 0 0 0 0 39 9294 4784 9610 9402 0 0 1 0 0 0 0 0 0 0 40 9448 4975 9294 9610 0 0 0 1 0 0 0 0 0 0 41 10319 4607 9448 9294 0 0 0 0 1 0 0 0 0 0 42 9548 6249 10319 9448 0 0 0 0 0 1 0 0 0 0 43 9801 4809 9548 10319 0 0 0 0 0 0 1 0 0 0 44 9596 3157 9801 9548 0 0 0 0 0 0 0 1 0 0 45 8923 1910 9596 9801 0 0 0 0 0 0 0 0 1 0 46 9746 2228 8923 9596 0 0 0 0 0 0 0 0 0 1 47 9829 1594 9746 8923 0 0 0 0 0 0 0 0 0 0 48 9125 2467 9829 9746 0 0 0 0 0 0 0 0 0 0 49 9782 2222 9125 9829 1 0 0 0 0 0 0 0 0 0 50 9441 3607 9782 9125 0 1 0 0 0 0 0 0 0 0 51 9162 4685 9441 9782 0 0 1 0 0 0 0 0 0 0 52 9915 4962 9162 9441 0 0 0 1 0 0 0 0 0 0 53 10444 5770 9915 9162 0 0 0 0 1 0 0 0 0 0 54 10209 5480 10444 9915 0 0 0 0 0 1 0 0 0 0 55 9985 5000 10209 10444 0 0 0 0 0 0 1 0 0 0 56 9842 3228 9985 10209 0 0 0 0 0 0 0 1 0 0 57 9429 1993 9842 9985 0 0 0 0 0 0 0 0 1 0 58 10132 2288 9429 9842 0 0 0 0 0 0 0 0 0 1 59 9849 1580 10132 9429 0 0 0 0 0 0 0 0 0 0 60 9172 2111 9849 10132 0 0 0 0 0 0 0 0 0 0 61 10313 2192 9172 9849 1 0 0 0 0 0 0 0 0 0 62 9819 3601 10313 9172 0 1 0 0 0 0 0 0 0 0 63 9955 4665 9819 10313 0 0 1 0 0 0 0 0 0 0 64 10048 4876 9955 9819 0 0 0 1 0 0 0 0 0 0 65 10082 5813 10048 9955 0 0 0 0 1 0 0 0 0 0 66 10541 5589 10082 10048 0 0 0 0 0 1 0 0 0 0 67 10208 5331 10541 10082 0 0 0 0 0 0 1 0 0 0 68 10233 3075 10208 10541 0 0 0 0 0 0 0 1 0 0 69 9439 2002 10233 10208 0 0 0 0 0 0 0 0 1 0 70 9963 2306 9439 10233 0 0 0 0 0 0 0 0 0 1 71 10158 1507 9963 9439 0 0 0 0 0 0 0 0 0 0 72 9225 1992 10158 9963 0 0 0 0 0 0 0 0 0 0 73 10474 2487 9225 10158 1 0 0 0 0 0 0 0 0 0 74 9757 3490 10474 9225 0 1 0 0 0 0 0 0 0 0 75 10490 4647 9757 10474 0 0 1 0 0 0 0 0 0 0 76 10281 5594 10490 9757 0 0 0 1 0 0 0 0 0 0 77 10444 5611 10281 10490 0 0 0 0 1 0 0 0 0 0 78 10640 5788 10444 10281 0 0 0 0 0 1 0 0 0 0 79 10695 6204 10640 10444 0 0 0 0 0 0 1 0 0 0 80 10786 3013 10695 10640 0 0 0 0 0 0 0 1 0 0 81 9832 1931 10786 10695 0 0 0 0 0 0 0 0 1 0 82 9747 2549 9832 10786 0 0 0 0 0 0 0 0 0 1 83 10411 1504 9747 9832 0 0 0 0 0 0 0 0 0 0 84 9511 2090 10411 9747 0 0 0 0 0 0 0 0 0 0 85 10402 2702 9511 10411 1 0 0 0 0 0 0 0 0 0 86 9701 2939 10402 9511 0 1 0 0 0 0 0 0 0 0 87 10540 4500 9701 10402 0 0 1 0 0 0 0 0 0 0 88 10112 6208 10540 9701 0 0 0 1 0 0 0 0 0 0 89 10915 6415 10112 10540 0 0 0 0 1 0 0 0 0 0 90 11183 5657 10915 10112 0 0 0 0 0 1 0 0 0 0 91 10384 5964 11183 10915 0 0 0 0 0 0 1 0 0 0 92 10834 3163 10384 11183 0 0 0 0 0 0 0 1 0 0 93 9886 1997 10834 10384 0 0 0 0 0 0 0 0 1 0 94 10216 2422 9886 10834 0 0 0 0 0 0 0 0 0 1 M11 t 1 0 1 2 0 2 3 0 3 4 0 4 5 0 5 6 0 6 7 0 7 8 0 8 9 0 9 10 0 10 11 1 11 12 0 12 13 0 13 14 0 14 15 0 15 16 0 16 17 0 17 18 0 18 19 0 19 20 0 20 21 0 21 22 0 22 23 1 23 24 0 24 25 0 25 26 0 26 27 0 27 28 0 28 29 0 29 30 0 30 31 0 31 32 0 32 33 0 33 34 0 34 35 1 35 36 0 36 37 0 37 38 0 38 39 0 39 40 0 40 41 0 41 42 0 42 43 0 43 44 0 44 45 0 45 46 0 46 47 1 47 48 0 48 49 0 49 50 0 50 51 0 51 52 0 52 53 0 53 54 0 54 55 0 55 56 0 56 57 0 57 58 0 58 59 1 59 60 0 60 61 0 61 62 0 62 63 0 63 64 0 64 65 0 65 66 0 66 67 0 67 68 0 68 69 0 69 70 0 70 71 1 71 72 0 72 73 0 73 74 0 74 75 0 75 76 0 76 77 0 77 78 0 78 79 0 79 80 0 80 81 0 81 82 0 82 83 1 83 84 0 84 85 0 85 86 0 86 87 0 87 88 0 88 89 0 89 90 0 90 91 0 91 92 0 92 93 0 93 94 0 94 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) huwelijken `geboortes-1` `geboortes-2` M1 3636.4907 -0.1129 0.2637 0.2880 1160.8781 M2 M3 M4 M5 M6 804.7805 1154.0644 1093.9421 1625.1079 1529.2834 M7 M8 M9 M10 M11 984.2230 980.8298 147.7561 717.4243 1000.9052 t 5.0797 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -665.715 -139.046 -6.872 160.798 610.058 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3636.49072 1159.93565 3.135 0.002422 ** huwelijken -0.11289 0.08624 -1.309 0.194366 `geboortes-1` 0.26368 0.11350 2.323 0.022777 * `geboortes-2` 0.28804 0.10331 2.788 0.006656 ** M1 1160.87809 179.01480 6.485 7.39e-09 *** M2 804.78050 173.40809 4.641 1.38e-05 *** M3 1154.06440 273.38315 4.221 6.51e-05 *** M4 1093.94206 308.80929 3.542 0.000673 *** M5 1625.10788 324.63993 5.006 3.37e-06 *** M6 1529.28340 331.64297 4.611 1.54e-05 *** M7 984.22302 310.68344 3.168 0.002192 ** M8 980.82980 169.08708 5.801 1.34e-07 *** M9 147.75608 141.15335 1.047 0.298435 M10 717.42434 167.25797 4.289 5.09e-05 *** M11 1000.90518 153.90816 6.503 6.83e-09 *** t 5.07968 1.51901 3.344 0.001271 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 262.9 on 78 degrees of freedom Multiple R-squared: 0.7821, Adjusted R-squared: 0.7401 F-statistic: 18.66 on 15 and 78 DF, p-value: < 2.2e-16 > 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.20118263 0.40236527 0.79881737 [2,] 0.10477399 0.20954799 0.89522601 [3,] 0.07067286 0.14134573 0.92932714 [4,] 0.03407870 0.06815740 0.96592130 [5,] 0.03695076 0.07390152 0.96304924 [6,] 0.04639208 0.09278416 0.95360792 [7,] 0.02967776 0.05935552 0.97032224 [8,] 0.04523780 0.09047561 0.95476220 [9,] 0.25113397 0.50226794 0.74886603 [10,] 0.25320227 0.50640453 0.74679773 [11,] 0.20582747 0.41165494 0.79417253 [12,] 0.16155453 0.32310906 0.83844547 [13,] 0.13829494 0.27658987 0.86170506 [14,] 0.10326471 0.20652942 0.89673529 [15,] 0.07557841 0.15115683 0.92442159 [16,] 0.09413043 0.18826086 0.90586957 [17,] 0.06699254 0.13398509 0.93300746 [18,] 0.04819686 0.09639372 0.95180314 [19,] 0.03791808 0.07583616 0.96208192 [20,] 0.14332107 0.28664213 0.85667893 [21,] 0.15404148 0.30808296 0.84595852 [22,] 0.16888860 0.33777720 0.83111140 [23,] 0.24409699 0.48819397 0.75590301 [24,] 0.26481029 0.52962058 0.73518971 [25,] 0.32993928 0.65987856 0.67006072 [26,] 0.33359351 0.66718701 0.66640649 [27,] 0.30611115 0.61222230 0.69388885 [28,] 0.43435304 0.86870607 0.56564696 [29,] 0.41836568 0.83673135 0.58163432 [30,] 0.51420135 0.97159731 0.48579865 [31,] 0.46046593 0.92093186 0.53953407 [32,] 0.40783131 0.81566262 0.59216869 [33,] 0.79173253 0.41653494 0.20826747 [34,] 0.84327030 0.31345940 0.15672970 [35,] 0.87686327 0.24627345 0.12313673 [36,] 0.85261189 0.29477621 0.14738811 [37,] 0.82954346 0.34091308 0.17045654 [38,] 0.88657633 0.22684735 0.11342367 [39,] 0.86033235 0.27933531 0.13966765 [40,] 0.93399823 0.13200354 0.06600177 [41,] 0.90925220 0.18149560 0.09074780 [42,] 0.87633145 0.24733711 0.12366855 [43,] 0.86423681 0.27152639 0.13576319 [44,] 0.89839135 0.20321731 0.10160865 [45,] 0.87032022 0.25935957 0.12967978 [46,] 0.85224130 0.29551740 0.14775870 [47,] 0.87946570 0.24106861 0.12053430 [48,] 0.83610760 0.32778480 0.16389240 [49,] 0.79984813 0.40030374 0.20015187 [50,] 0.83872709 0.32254583 0.16127291 [51,] 0.84815085 0.30369831 0.15184915 [52,] 0.77380667 0.45238665 0.22619333 [53,] 0.71243413 0.57513173 0.28756587 [54,] 0.59773378 0.80453245 0.40226622 [55,] 0.48749745 0.97499490 0.51250255 [56,] 0.35628780 0.71257561 0.64371220 [57,] 0.27725919 0.55451838 0.72274081 > postscript(file="/var/www/html/rcomp/tmp/1ri7j1291545234.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/2ri7j1291545234.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/3ri7j1291545234.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/4jr6m1291545234.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/5jr6m1291545234.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 = 94 Frequency = 1 1 2 3 4 5 6 142.881936 -3.820038 610.057756 -84.511899 -3.438908 301.571139 7 8 9 10 11 12 35.961299 98.039513 309.270237 -49.642739 248.887644 -101.790744 13 14 15 16 17 18 127.097591 340.964712 330.703005 -92.225554 -128.162874 -13.601740 19 20 21 22 23 24 -320.445398 132.340361 196.022049 -210.713676 -41.816605 180.114047 25 26 27 28 29 30 -45.436655 -326.899711 -204.638466 -397.290360 74.360779 -40.289480 31 32 33 34 35 36 -286.081469 6.923025 -215.221397 -26.561953 -6.554596 -186.761349 37 38 39 40 41 42 -371.271110 349.186054 -396.691028 -142.673230 200.948673 -567.962513 43 44 45 46 47 48 14.873898 -222.947939 -227.548905 293.109275 -7.188936 124.253045 49 50 51 52 53 54 -250.637181 -54.726072 -665.714509 345.387184 311.167156 -222.206117 55 56 57 58 59 60 -50.818993 -268.797876 109.000361 420.646106 -297.253338 -46.348249 61 62 63 64 65 66 197.866027 108.086792 -188.547683 89.743851 -370.416585 118.287730 67 68 69 70 71 72 165.318184 -110.455797 -108.271919 77.463333 -15.768525 -100.538883 73 74 75 76 77 78 228.235036 -85.119393 253.438753 219.632137 -307.714261 16.231528 79 80 81 82 83 84 559.543849 217.659312 -70.333676 -434.971095 119.694356 131.072133 85 86 87 88 89 90 -28.735644 -327.672342 261.392172 61.937870 223.256020 407.969454 91 92 93 94 -118.351371 147.239402 7.083251 -69.329252 > postscript(file="/var/www/html/rcomp/tmp/6jr6m1291545234.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 = 94 Frequency = 1 lag(myerror, k = 1) myerror 0 142.881936 NA 1 -3.820038 142.881936 2 610.057756 -3.820038 3 -84.511899 610.057756 4 -3.438908 -84.511899 5 301.571139 -3.438908 6 35.961299 301.571139 7 98.039513 35.961299 8 309.270237 98.039513 9 -49.642739 309.270237 10 248.887644 -49.642739 11 -101.790744 248.887644 12 127.097591 -101.790744 13 340.964712 127.097591 14 330.703005 340.964712 15 -92.225554 330.703005 16 -128.162874 -92.225554 17 -13.601740 -128.162874 18 -320.445398 -13.601740 19 132.340361 -320.445398 20 196.022049 132.340361 21 -210.713676 196.022049 22 -41.816605 -210.713676 23 180.114047 -41.816605 24 -45.436655 180.114047 25 -326.899711 -45.436655 26 -204.638466 -326.899711 27 -397.290360 -204.638466 28 74.360779 -397.290360 29 -40.289480 74.360779 30 -286.081469 -40.289480 31 6.923025 -286.081469 32 -215.221397 6.923025 33 -26.561953 -215.221397 34 -6.554596 -26.561953 35 -186.761349 -6.554596 36 -371.271110 -186.761349 37 349.186054 -371.271110 38 -396.691028 349.186054 39 -142.673230 -396.691028 40 200.948673 -142.673230 41 -567.962513 200.948673 42 14.873898 -567.962513 43 -222.947939 14.873898 44 -227.548905 -222.947939 45 293.109275 -227.548905 46 -7.188936 293.109275 47 124.253045 -7.188936 48 -250.637181 124.253045 49 -54.726072 -250.637181 50 -665.714509 -54.726072 51 345.387184 -665.714509 52 311.167156 345.387184 53 -222.206117 311.167156 54 -50.818993 -222.206117 55 -268.797876 -50.818993 56 109.000361 -268.797876 57 420.646106 109.000361 58 -297.253338 420.646106 59 -46.348249 -297.253338 60 197.866027 -46.348249 61 108.086792 197.866027 62 -188.547683 108.086792 63 89.743851 -188.547683 64 -370.416585 89.743851 65 118.287730 -370.416585 66 165.318184 118.287730 67 -110.455797 165.318184 68 -108.271919 -110.455797 69 77.463333 -108.271919 70 -15.768525 77.463333 71 -100.538883 -15.768525 72 228.235036 -100.538883 73 -85.119393 228.235036 74 253.438753 -85.119393 75 219.632137 253.438753 76 -307.714261 219.632137 77 16.231528 -307.714261 78 559.543849 16.231528 79 217.659312 559.543849 80 -70.333676 217.659312 81 -434.971095 -70.333676 82 119.694356 -434.971095 83 131.072133 119.694356 84 -28.735644 131.072133 85 -327.672342 -28.735644 86 261.392172 -327.672342 87 61.937870 261.392172 88 223.256020 61.937870 89 407.969454 223.256020 90 -118.351371 407.969454 91 147.239402 -118.351371 92 7.083251 147.239402 93 -69.329252 7.083251 94 NA -69.329252 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -3.820038 142.881936 [2,] 610.057756 -3.820038 [3,] -84.511899 610.057756 [4,] -3.438908 -84.511899 [5,] 301.571139 -3.438908 [6,] 35.961299 301.571139 [7,] 98.039513 35.961299 [8,] 309.270237 98.039513 [9,] -49.642739 309.270237 [10,] 248.887644 -49.642739 [11,] -101.790744 248.887644 [12,] 127.097591 -101.790744 [13,] 340.964712 127.097591 [14,] 330.703005 340.964712 [15,] -92.225554 330.703005 [16,] -128.162874 -92.225554 [17,] -13.601740 -128.162874 [18,] -320.445398 -13.601740 [19,] 132.340361 -320.445398 [20,] 196.022049 132.340361 [21,] -210.713676 196.022049 [22,] -41.816605 -210.713676 [23,] 180.114047 -41.816605 [24,] -45.436655 180.114047 [25,] -326.899711 -45.436655 [26,] -204.638466 -326.899711 [27,] -397.290360 -204.638466 [28,] 74.360779 -397.290360 [29,] -40.289480 74.360779 [30,] -286.081469 -40.289480 [31,] 6.923025 -286.081469 [32,] -215.221397 6.923025 [33,] -26.561953 -215.221397 [34,] -6.554596 -26.561953 [35,] -186.761349 -6.554596 [36,] -371.271110 -186.761349 [37,] 349.186054 -371.271110 [38,] -396.691028 349.186054 [39,] -142.673230 -396.691028 [40,] 200.948673 -142.673230 [41,] -567.962513 200.948673 [42,] 14.873898 -567.962513 [43,] -222.947939 14.873898 [44,] -227.548905 -222.947939 [45,] 293.109275 -227.548905 [46,] -7.188936 293.109275 [47,] 124.253045 -7.188936 [48,] -250.637181 124.253045 [49,] -54.726072 -250.637181 [50,] -665.714509 -54.726072 [51,] 345.387184 -665.714509 [52,] 311.167156 345.387184 [53,] -222.206117 311.167156 [54,] -50.818993 -222.206117 [55,] -268.797876 -50.818993 [56,] 109.000361 -268.797876 [57,] 420.646106 109.000361 [58,] -297.253338 420.646106 [59,] -46.348249 -297.253338 [60,] 197.866027 -46.348249 [61,] 108.086792 197.866027 [62,] -188.547683 108.086792 [63,] 89.743851 -188.547683 [64,] -370.416585 89.743851 [65,] 118.287730 -370.416585 [66,] 165.318184 118.287730 [67,] -110.455797 165.318184 [68,] -108.271919 -110.455797 [69,] 77.463333 -108.271919 [70,] -15.768525 77.463333 [71,] -100.538883 -15.768525 [72,] 228.235036 -100.538883 [73,] -85.119393 228.235036 [74,] 253.438753 -85.119393 [75,] 219.632137 253.438753 [76,] -307.714261 219.632137 [77,] 16.231528 -307.714261 [78,] 559.543849 16.231528 [79,] 217.659312 559.543849 [80,] -70.333676 217.659312 [81,] -434.971095 -70.333676 [82,] 119.694356 -434.971095 [83,] 131.072133 119.694356 [84,] -28.735644 131.072133 [85,] -327.672342 -28.735644 [86,] 261.392172 -327.672342 [87,] 61.937870 261.392172 [88,] 223.256020 61.937870 [89,] 407.969454 223.256020 [90,] -118.351371 407.969454 [91,] 147.239402 -118.351371 [92,] 7.083251 147.239402 [93,] -69.329252 7.083251 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -3.820038 142.881936 2 610.057756 -3.820038 3 -84.511899 610.057756 4 -3.438908 -84.511899 5 301.571139 -3.438908 6 35.961299 301.571139 7 98.039513 35.961299 8 309.270237 98.039513 9 -49.642739 309.270237 10 248.887644 -49.642739 11 -101.790744 248.887644 12 127.097591 -101.790744 13 340.964712 127.097591 14 330.703005 340.964712 15 -92.225554 330.703005 16 -128.162874 -92.225554 17 -13.601740 -128.162874 18 -320.445398 -13.601740 19 132.340361 -320.445398 20 196.022049 132.340361 21 -210.713676 196.022049 22 -41.816605 -210.713676 23 180.114047 -41.816605 24 -45.436655 180.114047 25 -326.899711 -45.436655 26 -204.638466 -326.899711 27 -397.290360 -204.638466 28 74.360779 -397.290360 29 -40.289480 74.360779 30 -286.081469 -40.289480 31 6.923025 -286.081469 32 -215.221397 6.923025 33 -26.561953 -215.221397 34 -6.554596 -26.561953 35 -186.761349 -6.554596 36 -371.271110 -186.761349 37 349.186054 -371.271110 38 -396.691028 349.186054 39 -142.673230 -396.691028 40 200.948673 -142.673230 41 -567.962513 200.948673 42 14.873898 -567.962513 43 -222.947939 14.873898 44 -227.548905 -222.947939 45 293.109275 -227.548905 46 -7.188936 293.109275 47 124.253045 -7.188936 48 -250.637181 124.253045 49 -54.726072 -250.637181 50 -665.714509 -54.726072 51 345.387184 -665.714509 52 311.167156 345.387184 53 -222.206117 311.167156 54 -50.818993 -222.206117 55 -268.797876 -50.818993 56 109.000361 -268.797876 57 420.646106 109.000361 58 -297.253338 420.646106 59 -46.348249 -297.253338 60 197.866027 -46.348249 61 108.086792 197.866027 62 -188.547683 108.086792 63 89.743851 -188.547683 64 -370.416585 89.743851 65 118.287730 -370.416585 66 165.318184 118.287730 67 -110.455797 165.318184 68 -108.271919 -110.455797 69 77.463333 -108.271919 70 -15.768525 77.463333 71 -100.538883 -15.768525 72 228.235036 -100.538883 73 -85.119393 228.235036 74 253.438753 -85.119393 75 219.632137 253.438753 76 -307.714261 219.632137 77 16.231528 -307.714261 78 559.543849 16.231528 79 217.659312 559.543849 80 -70.333676 217.659312 81 -434.971095 -70.333676 82 119.694356 -434.971095 83 131.072133 119.694356 84 -28.735644 131.072133 85 -327.672342 -28.735644 86 261.392172 -327.672342 87 61.937870 261.392172 88 223.256020 61.937870 89 407.969454 223.256020 90 -118.351371 407.969454 91 147.239402 -118.351371 92 7.083251 147.239402 93 -69.329252 7.083251 > 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/7c05p1291545234.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/8c05p1291545234.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/9mr591291545234.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/10mr591291545234.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/11qs4g1291545234.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/12tb231291545234.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/13iuzx1291545234.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/14blyi1291545234.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/15e3fo1291545234.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/16sdcx1291545234.tab") + } > > try(system("convert tmp/1ri7j1291545234.ps tmp/1ri7j1291545234.png",intern=TRUE)) character(0) > try(system("convert tmp/2ri7j1291545234.ps tmp/2ri7j1291545234.png",intern=TRUE)) character(0) > try(system("convert tmp/3ri7j1291545234.ps tmp/3ri7j1291545234.png",intern=TRUE)) character(0) > try(system("convert tmp/4jr6m1291545234.ps tmp/4jr6m1291545234.png",intern=TRUE)) character(0) > try(system("convert tmp/5jr6m1291545234.ps tmp/5jr6m1291545234.png",intern=TRUE)) character(0) > try(system("convert tmp/6jr6m1291545234.ps tmp/6jr6m1291545234.png",intern=TRUE)) character(0) > try(system("convert tmp/7c05p1291545234.ps tmp/7c05p1291545234.png",intern=TRUE)) character(0) > try(system("convert tmp/8c05p1291545234.ps tmp/8c05p1291545234.png",intern=TRUE)) character(0) > try(system("convert tmp/9mr591291545234.ps tmp/9mr591291545234.png",intern=TRUE)) character(0) > try(system("convert tmp/10mr591291545234.ps tmp/10mr591291545234.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.023 1.679 7.898