R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i686-pc-linux-gnu (32-bit) 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(10951 + ,3113 + ,3193 + ,1714 + ,1811 + ,542 + ,577 + ,10840 + ,3085 + ,3167 + ,1700 + ,1798 + ,527 + ,563 + ,10753 + ,3064 + ,3145 + ,1687 + ,1789 + ,516 + ,553 + ,10667 + ,3040 + ,3122 + ,1678 + ,1778 + ,506 + ,543 + ,10585 + ,3017 + ,3100 + ,1668 + ,1768 + ,497 + ,534 + ,10511 + ,2997 + ,3081 + ,1657 + ,1757 + ,490 + ,529 + ,10446 + ,2980 + ,3063 + ,1648 + ,1748 + ,484 + ,523 + ,10396 + ,2967 + ,3049 + ,1640 + ,1740 + ,480 + ,520 + ,10355 + ,2957 + ,3038 + ,1634 + ,1735 + ,477 + ,515 + ,10310 + ,2945 + ,3028 + ,1628 + ,1730 + ,469 + ,510 + ,10263 + ,2935 + ,3018 + ,1622 + ,1724 + ,461 + ,503 + ,10239 + ,2930 + ,3011 + ,1619 + ,1721 + ,458 + ,501 + ,10214 + ,2924 + ,3003 + ,1615 + ,1717 + ,455 + ,500 + ,10192 + ,2917 + ,2996 + ,1612 + ,1714 + ,454 + ,499 + ,10170 + ,2911 + ,2988 + ,1609 + ,1711 + ,452 + ,499 + ,10143 + ,2902 + ,2978 + ,1607 + ,1708 + ,450 + ,498 + ,10131 + ,2897 + ,2969 + ,1607 + ,1706 + ,452 + ,500 + ,10101 + ,2888 + ,2959 + ,1603 + ,1702 + ,449 + ,500 + ,10068 + ,2877 + ,2948 + ,1597 + ,1696 + ,449 + ,501 + ,10022 + ,2862 + ,2933 + ,1588 + ,1688 + ,449 + ,502 + ,9987 + ,2849 + ,2919 + ,1579 + ,1680 + ,453 + ,508 + ,9948 + ,2835 + ,2905 + ,1572 + ,1672 + ,453 + ,511 + ,9928 + ,2826 + ,2896 + ,1567 + ,1668 + ,456 + ,514 + ,9876 + ,2813 + ,2883 + ,1554 + ,1656 + ,455 + ,515 + ,9865 + ,2808 + ,2877 + ,1552 + ,1654 + ,456 + ,517 + ,9859 + ,2803 + ,2873 + ,1551 + ,1655 + ,457 + ,519 + ,9858 + ,2801 + ,2869 + ,1552 + ,1656 + ,458 + ,521 + ,9853 + ,2798 + ,2864 + ,1552 + ,1656 + ,460 + ,523 + ,9858 + ,2795 + ,2860 + ,1554 + ,1659 + ,464 + ,526 + ,9855 + ,2789 + ,2853 + ,1557 + ,1661 + ,466 + ,528 + ,9863 + ,2788 + ,2847 + ,1564 + ,1665 + ,469 + ,531 + ,9855 + ,2781 + ,2838 + ,1564 + ,1664 + ,474 + ,535 + ,9842 + ,2774 + ,2827 + ,1563 + ,1662 + ,477 + ,539 + ,9837 + ,2767 + ,2817 + ,1563 + ,1661 + ,484 + ,545 + ,9823 + ,2759 + ,2807 + ,1559 + ,1656 + ,490 + ,552 + ,9813 + ,2752 + ,2796 + ,1559 + ,1654 + ,494 + ,557 + ,9788 + ,2741 + ,2786 + ,1556 + ,1650 + ,495 + ,560 + ,9757 + ,2728 + ,2773 + ,1549 + ,1643 + ,498 + ,566 + ,9727 + ,2717 + ,2761 + ,1543 + ,1637 + ,500 + ,569 + ,9695 + ,2705 + ,2747 + ,1538 + ,1632 + ,502 + ,572 + ,9651 + ,2687 + ,2729 + ,1533 + ,1627 + ,502 + ,573 + ,9660 + ,2683 + ,2721 + ,1546 + ,1637 + ,501 + ,572 + ,9632 + ,2668 + ,2705 + ,1547 + ,1634 + ,503 + ,574 + ,9606 + ,2657 + ,2690 + ,1547 + ,1632 + ,504 + ,575 + ,9556 + ,2638 + ,2670 + ,1547 + ,1627 + ,502 + ,572 + ,9499 + ,2617 + ,2647 + ,1547 + ,1622 + ,498 + ,568 + ,9428 + ,2595 + ,2623 + ,1540 + ,1613 + ,493 + ,565 + ,9328 + ,2564 + ,2597 + ,1525 + ,1602 + ,483 + ,558 + ,9251 + ,2537 + ,2572 + ,1517 + ,1595 + ,476 + ,554 + ,9190 + ,2514 + ,2550 + ,1512 + ,1591 + ,472 + ,551) + ,dim=c(7 + ,50) + ,dimnames=list(c('Totale_bevolking' + ,'Vlaams_Gewest_mannen' + ,'Vlaams_Gewest_vrouwen' + ,'Waals_Gewest_manen' + ,'Waals_Gewest_vrouwen' + ,'Brussel_mannen' + ,'Brussel_vrouwen') + ,1:50)) > y <- array(NA,dim=c(7,50),dimnames=list(c('Totale_bevolking','Vlaams_Gewest_mannen','Vlaams_Gewest_vrouwen','Waals_Gewest_manen','Waals_Gewest_vrouwen','Brussel_mannen','Brussel_vrouwen'),1:50)) > for (i in 1:dim(x)[1]) + { + for (j in 1:dim(x)[2]) + { + y[i,j] <- as.numeric(x[i,j]) + } + } > par3 = 'No Linear Trend' > par2 = 'Do not include Seasonal Dummies' > par1 = '1' > par3 <- 'No Linear Trend' > par2 <- 'Do not include Seasonal 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, 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 Totale_bevolking Vlaams_Gewest_mannen Vlaams_Gewest_vrouwen 1 10951 3113 3193 2 10840 3085 3167 3 10753 3064 3145 4 10667 3040 3122 5 10585 3017 3100 6 10511 2997 3081 7 10446 2980 3063 8 10396 2967 3049 9 10355 2957 3038 10 10310 2945 3028 11 10263 2935 3018 12 10239 2930 3011 13 10214 2924 3003 14 10192 2917 2996 15 10170 2911 2988 16 10143 2902 2978 17 10131 2897 2969 18 10101 2888 2959 19 10068 2877 2948 20 10022 2862 2933 21 9987 2849 2919 22 9948 2835 2905 23 9928 2826 2896 24 9876 2813 2883 25 9865 2808 2877 26 9859 2803 2873 27 9858 2801 2869 28 9853 2798 2864 29 9858 2795 2860 30 9855 2789 2853 31 9863 2788 2847 32 9855 2781 2838 33 9842 2774 2827 34 9837 2767 2817 35 9823 2759 2807 36 9813 2752 2796 37 9788 2741 2786 38 9757 2728 2773 39 9727 2717 2761 40 9695 2705 2747 41 9651 2687 2729 42 9660 2683 2721 43 9632 2668 2705 44 9606 2657 2690 45 9556 2638 2670 46 9499 2617 2647 47 9428 2595 2623 48 9328 2564 2597 49 9251 2537 2572 50 9190 2514 2550 Waals_Gewest_manen Waals_Gewest_vrouwen Brussel_mannen Brussel_vrouwen 1 1714 1811 542 577 2 1700 1798 527 563 3 1687 1789 516 553 4 1678 1778 506 543 5 1668 1768 497 534 6 1657 1757 490 529 7 1648 1748 484 523 8 1640 1740 480 520 9 1634 1735 477 515 10 1628 1730 469 510 11 1622 1724 461 503 12 1619 1721 458 501 13 1615 1717 455 500 14 1612 1714 454 499 15 1609 1711 452 499 16 1607 1708 450 498 17 1607 1706 452 500 18 1603 1702 449 500 19 1597 1696 449 501 20 1588 1688 449 502 21 1579 1680 453 508 22 1572 1672 453 511 23 1567 1668 456 514 24 1554 1656 455 515 25 1552 1654 456 517 26 1551 1655 457 519 27 1552 1656 458 521 28 1552 1656 460 523 29 1554 1659 464 526 30 1557 1661 466 528 31 1564 1665 469 531 32 1564 1664 474 535 33 1563 1662 477 539 34 1563 1661 484 545 35 1559 1656 490 552 36 1559 1654 494 557 37 1556 1650 495 560 38 1549 1643 498 566 39 1543 1637 500 569 40 1538 1632 502 572 41 1533 1627 502 573 42 1546 1637 501 572 43 1547 1634 503 574 44 1547 1632 504 575 45 1547 1627 502 572 46 1547 1622 498 568 47 1540 1613 493 565 48 1525 1602 483 558 49 1517 1595 476 554 50 1512 1591 472 551 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) Vlaams_Gewest_mannen Vlaams_Gewest_vrouwen -15.6002 0.9765 1.0300 Waals_Gewest_manen Waals_Gewest_vrouwen Brussel_mannen 1.0199 0.9740 0.9301 Brussel_vrouwen 1.0769 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -1.23005 -0.27032 0.08218 0.21548 1.15485 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -15.60018 22.71234 -0.687 0.496 Vlaams_Gewest_mannen 0.97647 0.03002 32.530 <2e-16 *** Vlaams_Gewest_vrouwen 1.02997 0.02795 36.847 <2e-16 *** Waals_Gewest_manen 1.01993 0.03544 28.780 <2e-16 *** Waals_Gewest_vrouwen 0.97402 0.05014 19.426 <2e-16 *** Brussel_mannen 0.93010 0.05695 16.331 <2e-16 *** Brussel_vrouwen 1.07691 0.05123 21.020 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.6203 on 43 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 3.302e+06 on 6 and 43 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.48509082 0.97018163 0.5149092 [2,] 0.40562242 0.81124484 0.5943776 [3,] 0.53858489 0.92283023 0.4614151 [4,] 0.43066581 0.86133162 0.5693342 [5,] 0.31014821 0.62029642 0.6898518 [6,] 0.24031791 0.48063583 0.7596821 [7,] 0.23704799 0.47409598 0.7629520 [8,] 0.18626734 0.37253468 0.8137327 [9,] 0.13690244 0.27380487 0.8630976 [10,] 0.08596361 0.17192721 0.9140364 [11,] 0.05883437 0.11766875 0.9411656 [12,] 0.06083335 0.12166669 0.9391667 [13,] 0.04597138 0.09194276 0.9540286 [14,] 0.13366985 0.26733971 0.8663301 [15,] 0.09516927 0.19033853 0.9048307 [16,] 0.12391893 0.24783787 0.8760811 [17,] 0.12708815 0.25417631 0.8729118 [18,] 0.17748324 0.35496647 0.8225168 [19,] 0.19047802 0.38095605 0.8095220 [20,] 0.15403779 0.30807558 0.8459622 [21,] 0.49698637 0.99397273 0.5030136 [22,] 0.54724558 0.90550883 0.4527544 [23,] 0.60779939 0.78440122 0.3922006 [24,] 0.58011271 0.83977457 0.4198873 [25,] 0.59661175 0.80677651 0.4033883 [26,] 0.62594685 0.74810631 0.3740532 [27,] 0.81875682 0.36248635 0.1812432 [28,] 0.81401755 0.37196491 0.1859825 [29,] 0.73409102 0.53181796 0.2659090 [30,] 0.89724824 0.20550351 0.1027518 [31,] 0.81232716 0.37534569 0.1876728 > postscript(file="/var/fisher/rcomp/tmp/11o9v1351956747.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/fisher/rcomp/tmp/2qsv01351956747.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/fisher/rcomp/tmp/3o4sf1351956747.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/fisher/rcomp/tmp/4pqbt1351956747.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/fisher/rcomp/tmp/5g3fo1351956747.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 = 50 Frequency = 1 1 2 3 4 5 6 0.558206288 -0.351932953 -1.161299916 -0.073039052 1.047677082 -0.024823157 7 8 9 10 11 12 0.102207720 0.118600257 -0.622520863 0.209834186 0.217081087 -0.764824642 13 14 15 16 17 18 0.176746313 0.210680134 0.151296765 0.138254265 0.224351861 0.078369366 19 20 21 22 23 24 0.035985226 0.027182190 -1.069482031 -0.278396962 0.754253970 -0.361551315 25 26 27 28 29 30 0.604594115 0.568819760 0.563772909 -0.370986545 -0.234732021 0.812048948 31 32 33 34 35 36 -1.088257339 -0.967349084 0.067637581 0.204495157 0.146797034 1.154849226 37 38 39 40 41 42 -0.009237233 -0.219640408 -0.246075228 -1.230046658 -0.221287156 -0.067872566 43 44 45 46 47 48 0.946830002 1.078584289 0.191963126 0.285300117 -0.726174890 -0.823871719 49 50 0.085983989 0.150998776 > postscript(file="/var/fisher/rcomp/tmp/6bggf1351956747.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 = 50 Frequency = 1 lag(myerror, k = 1) myerror 0 0.558206288 NA 1 -0.351932953 0.558206288 2 -1.161299916 -0.351932953 3 -0.073039052 -1.161299916 4 1.047677082 -0.073039052 5 -0.024823157 1.047677082 6 0.102207720 -0.024823157 7 0.118600257 0.102207720 8 -0.622520863 0.118600257 9 0.209834186 -0.622520863 10 0.217081087 0.209834186 11 -0.764824642 0.217081087 12 0.176746313 -0.764824642 13 0.210680134 0.176746313 14 0.151296765 0.210680134 15 0.138254265 0.151296765 16 0.224351861 0.138254265 17 0.078369366 0.224351861 18 0.035985226 0.078369366 19 0.027182190 0.035985226 20 -1.069482031 0.027182190 21 -0.278396962 -1.069482031 22 0.754253970 -0.278396962 23 -0.361551315 0.754253970 24 0.604594115 -0.361551315 25 0.568819760 0.604594115 26 0.563772909 0.568819760 27 -0.370986545 0.563772909 28 -0.234732021 -0.370986545 29 0.812048948 -0.234732021 30 -1.088257339 0.812048948 31 -0.967349084 -1.088257339 32 0.067637581 -0.967349084 33 0.204495157 0.067637581 34 0.146797034 0.204495157 35 1.154849226 0.146797034 36 -0.009237233 1.154849226 37 -0.219640408 -0.009237233 38 -0.246075228 -0.219640408 39 -1.230046658 -0.246075228 40 -0.221287156 -1.230046658 41 -0.067872566 -0.221287156 42 0.946830002 -0.067872566 43 1.078584289 0.946830002 44 0.191963126 1.078584289 45 0.285300117 0.191963126 46 -0.726174890 0.285300117 47 -0.823871719 -0.726174890 48 0.085983989 -0.823871719 49 0.150998776 0.085983989 50 NA 0.150998776 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -0.351932953 0.558206288 [2,] -1.161299916 -0.351932953 [3,] -0.073039052 -1.161299916 [4,] 1.047677082 -0.073039052 [5,] -0.024823157 1.047677082 [6,] 0.102207720 -0.024823157 [7,] 0.118600257 0.102207720 [8,] -0.622520863 0.118600257 [9,] 0.209834186 -0.622520863 [10,] 0.217081087 0.209834186 [11,] -0.764824642 0.217081087 [12,] 0.176746313 -0.764824642 [13,] 0.210680134 0.176746313 [14,] 0.151296765 0.210680134 [15,] 0.138254265 0.151296765 [16,] 0.224351861 0.138254265 [17,] 0.078369366 0.224351861 [18,] 0.035985226 0.078369366 [19,] 0.027182190 0.035985226 [20,] -1.069482031 0.027182190 [21,] -0.278396962 -1.069482031 [22,] 0.754253970 -0.278396962 [23,] -0.361551315 0.754253970 [24,] 0.604594115 -0.361551315 [25,] 0.568819760 0.604594115 [26,] 0.563772909 0.568819760 [27,] -0.370986545 0.563772909 [28,] -0.234732021 -0.370986545 [29,] 0.812048948 -0.234732021 [30,] -1.088257339 0.812048948 [31,] -0.967349084 -1.088257339 [32,] 0.067637581 -0.967349084 [33,] 0.204495157 0.067637581 [34,] 0.146797034 0.204495157 [35,] 1.154849226 0.146797034 [36,] -0.009237233 1.154849226 [37,] -0.219640408 -0.009237233 [38,] -0.246075228 -0.219640408 [39,] -1.230046658 -0.246075228 [40,] -0.221287156 -1.230046658 [41,] -0.067872566 -0.221287156 [42,] 0.946830002 -0.067872566 [43,] 1.078584289 0.946830002 [44,] 0.191963126 1.078584289 [45,] 0.285300117 0.191963126 [46,] -0.726174890 0.285300117 [47,] -0.823871719 -0.726174890 [48,] 0.085983989 -0.823871719 [49,] 0.150998776 0.085983989 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -0.351932953 0.558206288 2 -1.161299916 -0.351932953 3 -0.073039052 -1.161299916 4 1.047677082 -0.073039052 5 -0.024823157 1.047677082 6 0.102207720 -0.024823157 7 0.118600257 0.102207720 8 -0.622520863 0.118600257 9 0.209834186 -0.622520863 10 0.217081087 0.209834186 11 -0.764824642 0.217081087 12 0.176746313 -0.764824642 13 0.210680134 0.176746313 14 0.151296765 0.210680134 15 0.138254265 0.151296765 16 0.224351861 0.138254265 17 0.078369366 0.224351861 18 0.035985226 0.078369366 19 0.027182190 0.035985226 20 -1.069482031 0.027182190 21 -0.278396962 -1.069482031 22 0.754253970 -0.278396962 23 -0.361551315 0.754253970 24 0.604594115 -0.361551315 25 0.568819760 0.604594115 26 0.563772909 0.568819760 27 -0.370986545 0.563772909 28 -0.234732021 -0.370986545 29 0.812048948 -0.234732021 30 -1.088257339 0.812048948 31 -0.967349084 -1.088257339 32 0.067637581 -0.967349084 33 0.204495157 0.067637581 34 0.146797034 0.204495157 35 1.154849226 0.146797034 36 -0.009237233 1.154849226 37 -0.219640408 -0.009237233 38 -0.246075228 -0.219640408 39 -1.230046658 -0.246075228 40 -0.221287156 -1.230046658 41 -0.067872566 -0.221287156 42 0.946830002 -0.067872566 43 1.078584289 0.946830002 44 0.191963126 1.078584289 45 0.285300117 0.191963126 46 -0.726174890 0.285300117 47 -0.823871719 -0.726174890 48 0.085983989 -0.823871719 49 0.150998776 0.085983989 > 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/fisher/rcomp/tmp/7kq5w1351956747.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/fisher/rcomp/tmp/8ou4a1351956747.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/fisher/rcomp/tmp/91mfh1351956747.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/fisher/rcomp/tmp/10pnsa1351956747.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/fisher/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/fisher/rcomp/tmp/11hddx1351956747.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/fisher/rcomp/tmp/12mc6e1351956747.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/fisher/rcomp/tmp/13w9141351956747.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/fisher/rcomp/tmp/14a8uu1351956747.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/fisher/rcomp/tmp/15hh9z1351956747.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/fisher/rcomp/tmp/16gsum1351956747.tab") + } > > try(system("convert tmp/11o9v1351956747.ps tmp/11o9v1351956747.png",intern=TRUE)) character(0) > try(system("convert tmp/2qsv01351956747.ps tmp/2qsv01351956747.png",intern=TRUE)) character(0) > try(system("convert tmp/3o4sf1351956747.ps tmp/3o4sf1351956747.png",intern=TRUE)) character(0) > try(system("convert tmp/4pqbt1351956747.ps tmp/4pqbt1351956747.png",intern=TRUE)) character(0) > try(system("convert tmp/5g3fo1351956747.ps tmp/5g3fo1351956747.png",intern=TRUE)) character(0) > try(system("convert tmp/6bggf1351956747.ps tmp/6bggf1351956747.png",intern=TRUE)) character(0) > try(system("convert tmp/7kq5w1351956747.ps tmp/7kq5w1351956747.png",intern=TRUE)) character(0) > try(system("convert tmp/8ou4a1351956747.ps tmp/8ou4a1351956747.png",intern=TRUE)) character(0) > try(system("convert tmp/91mfh1351956747.ps tmp/91mfh1351956747.png",intern=TRUE)) character(0) > try(system("convert tmp/10pnsa1351956747.ps tmp/10pnsa1351956747.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 5.716 1.149 6.860