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(31/12/1961 + ,9190 + ,0 + ,2514 + ,0 + ,2550 + ,0 + ,1512 + ,0 + ,1591 + ,0 + ,472 + ,0 + ,551 + ,0 + ,31/12/1962 + ,9251 + ,1 + ,2537 + ,2537 + ,2572 + ,2572 + ,1517 + ,1517 + ,1595 + ,1595 + ,476 + ,476 + ,554 + ,554 + ,31/12/1963 + ,9328 + ,2 + ,2564 + ,5128 + ,2597 + ,5194 + ,1525 + ,3050 + ,1602 + ,3204 + ,483 + ,966 + ,558 + ,1116 + ,31/12/1964 + ,9428 + ,0 + ,2595 + ,0 + ,2623 + ,0 + ,1540 + ,0 + ,1613 + ,0 + ,493 + ,0 + ,565 + ,0 + ,31/12/1965 + ,9499 + ,1 + ,2617 + ,2617 + ,2647 + ,2647 + ,1547 + ,1547 + ,1622 + ,1622 + ,498 + ,498 + ,568 + ,568 + ,31/12/1966 + ,9556 + ,2 + ,2638 + ,5276 + ,2670 + ,5340 + ,1547 + ,3094 + ,1627 + ,3254 + ,502 + ,1004 + ,572 + ,1144 + ,31/12/1967 + ,9606 + ,0 + ,2657 + ,0 + ,2690 + ,0 + ,1547 + ,0 + ,1632 + ,0 + ,504 + ,0 + ,575 + ,0 + ,31/12/1968 + ,9632 + ,1 + ,2668 + ,2668 + ,2705 + ,2705 + ,1547 + ,1547 + ,1634 + ,1634 + ,503 + ,503 + ,574 + ,574 + ,31/12/1969 + ,9660 + ,2 + ,2683 + ,5366 + ,2721 + ,5442 + ,1546 + ,3092 + ,1637 + ,3274 + ,501 + ,1002 + ,572 + ,1144 + ,31/12/1970 + ,9651 + ,0 + ,2687 + ,0 + ,2729 + ,0 + ,1533 + ,0 + ,1627 + ,0 + ,502 + ,0 + ,573 + ,0 + ,31/12/1971 + ,9695 + ,1 + ,2705 + ,2705 + ,2747 + ,2747 + ,1538 + ,1538 + ,1632 + ,1632 + ,502 + ,502 + ,572 + ,572 + ,31/12/1972 + ,9727 + ,2 + ,2717 + ,5434 + ,2761 + ,5522 + ,1543 + ,3086 + ,1637 + ,3274 + ,500 + ,1000 + ,569 + ,1138 + ,31/12/1973 + ,9757 + ,0 + ,2728 + ,0 + ,2773 + ,0 + ,1549 + ,0 + ,1643 + ,0 + ,498 + ,0 + ,566 + ,0 + ,31/12/1974 + ,9788 + ,1 + ,2741 + ,2741 + ,2786 + ,2786 + ,1556 + ,1556 + ,1650 + ,1650 + ,495 + ,495 + ,560 + ,560 + ,31/12/1975 + ,9813 + ,2 + ,2752 + ,5504 + ,2796 + ,5592 + ,1559 + ,3118 + ,1654 + ,3308 + ,494 + ,988 + ,557 + ,1114 + ,31/12/1976 + ,9823 + ,0 + ,2759 + ,0 + ,2807 + ,0 + ,1559 + ,0 + ,1656 + ,0 + ,490 + ,0 + ,552 + ,0 + ,31/12/1977 + ,9837 + ,1 + ,2767 + ,2767 + ,2817 + ,2817 + ,1563 + ,1563 + ,1661 + ,1661 + ,484 + ,484 + ,545 + ,545 + ,31/12/1978 + ,9842 + ,2 + ,2774 + ,5548 + ,2827 + ,5654 + ,1563 + ,3126 + ,1662 + ,3324 + ,477 + ,954 + ,539 + ,1078 + ,31/12/1979 + ,9855 + ,0 + ,2781 + ,0 + ,2838 + ,0 + ,1564 + ,0 + ,1664 + ,0 + ,474 + ,0 + ,535 + ,0 + ,31/12/1980 + ,9863 + ,1 + ,2788 + ,2788 + ,2847 + ,2847 + ,1564 + ,1564 + ,1665 + ,1665 + ,469 + ,469 + ,531 + ,531 + ,31/12/1981 + ,9855 + ,2 + ,2789 + ,5578 + ,2853 + ,5706 + ,1557 + ,3114 + ,1661 + ,3322 + ,466 + ,932 + ,528 + ,1056 + ,31/12/1982 + ,9858 + ,0 + ,2795 + ,0 + ,2860 + ,0 + ,1554 + ,0 + ,1659 + ,0 + ,464 + ,0 + ,526 + ,0 + ,31/12/1983 + ,9853 + ,1 + ,2798 + ,2798 + ,2864 + ,2864 + ,1552 + ,1552 + ,1656 + ,1656 + ,460 + ,460 + ,523 + ,523 + ,31/12/1984 + ,9858 + ,2 + ,2801 + ,5602 + ,2869 + ,5738 + ,1552 + ,3104 + ,1656 + ,3312 + ,458 + ,916 + ,521 + ,1042 + ,31/12/1985 + ,9859 + ,0 + ,2803 + ,0 + ,2873 + ,0 + ,1551 + ,0 + ,1655 + ,0 + ,457 + ,0 + ,519 + ,0 + ,31/12/1986 + ,9865 + ,1 + ,2808 + ,2808 + ,2877 + ,2877 + ,1552 + ,1552 + ,1654 + ,1654 + ,456 + ,456 + ,517 + ,517 + ,31/12/1987 + ,9876 + ,2 + ,2813 + ,5626 + ,2883 + ,5766 + ,1554 + ,3108 + ,1656 + ,3312 + ,455 + ,910 + ,515 + ,1030 + ,31/12/1988 + ,9928 + ,0 + ,2826 + ,0 + ,2896 + ,0 + ,1567 + ,0 + ,1668 + ,0 + ,456 + ,0 + ,514 + ,0 + ,31/12/1989 + ,9948 + ,1 + ,2835 + ,2835 + ,2905 + ,2905 + ,1572 + ,1572 + ,1672 + ,1672 + ,453 + ,453 + ,511 + ,511 + ,31/12/1990 + ,9987 + ,2 + ,2849 + ,5698 + ,2919 + ,5838 + ,1579 + ,3158 + ,1680 + ,3360 + ,453 + ,906 + ,508 + ,1016 + ,31/12/1991 + ,10022 + ,0 + ,2862 + ,0 + ,2933 + ,0 + ,1588 + ,0 + ,1688 + ,0 + ,449 + ,0 + ,502 + ,0 + ,31/12/1992 + ,10068 + ,1 + ,2877 + ,2877 + ,2948 + ,2948 + ,1597 + ,1597 + ,1696 + ,1696 + ,449 + ,449 + ,501 + ,501 + ,31/12/1993 + ,10101 + ,2 + ,2888 + ,5776 + ,2959 + ,5918 + ,1603 + ,3206 + ,1702 + ,3404 + ,449 + ,898 + ,500 + ,1000 + ,31/12/1994 + ,10131 + ,0 + ,2897 + ,0 + ,2969 + ,0 + ,1607 + ,0 + ,1706 + ,0 + ,452 + ,0 + ,500 + ,0 + ,31/12/1995 + ,10143 + ,1 + ,2902 + ,2902 + ,2978 + ,2978 + ,1607 + ,1607 + ,1708 + ,1708 + ,450 + ,450 + ,498 + ,498 + ,31/12/1996 + ,10170 + ,2 + ,2911 + ,5822 + ,2988 + ,5976 + ,1609 + ,3218 + ,1711 + ,3422 + ,452 + ,904 + ,499 + ,998 + ,31/12/1997 + ,10192 + ,0 + ,2917 + ,0 + ,2996 + ,0 + ,1612 + ,0 + ,1714 + ,0 + ,454 + ,0 + ,499 + ,0 + ,31/12/1998 + ,10214 + ,1 + ,2924 + ,2924 + ,3003 + ,3003 + ,1615 + ,1615 + ,1717 + ,1717 + ,455 + ,455 + ,500 + ,500 + ,31/12/1999 + ,10239 + ,2 + ,2930 + ,5860 + ,3011 + ,6022 + ,1619 + ,3238 + ,1721 + ,3442 + ,458 + ,916 + ,501 + ,1002 + ,31/12/2000 + ,10263 + ,0 + ,2935 + ,0 + ,3018 + ,0 + ,1622 + ,0 + ,1724 + ,0 + ,461 + ,0 + ,503 + ,0 + ,31/12/2001 + ,10310 + ,1 + ,2945 + ,2945 + ,3028 + ,3028 + ,1628 + ,1628 + ,1730 + ,1730 + ,469 + ,469 + ,510 + ,510 + ,31/12/2002 + ,10355 + ,2 + ,2957 + ,5914 + ,3038 + ,6076 + ,1634 + ,3268 + ,1735 + ,3470 + ,477 + ,954 + ,515 + ,1030 + ,31/12/2003 + ,10396 + ,0 + ,2967 + ,0 + ,3049 + ,0 + ,1640 + ,0 + ,1740 + ,0 + ,480 + ,0 + ,520 + ,0 + ,31/12/2004 + ,10446 + ,1 + ,2980 + ,2980 + ,3063 + ,3063 + ,1648 + ,1648 + ,1748 + ,1748 + ,484 + ,484 + ,523 + ,523 + ,31/12/2005 + ,10511 + ,2 + ,2997 + ,5994 + ,3081 + ,6162 + ,1657 + ,3314 + ,1757 + ,3514 + ,490 + ,980 + ,529 + ,1058 + ,31/12/2006 + ,10585 + ,0 + ,3017 + ,0 + ,3100 + ,0 + ,1668 + ,0 + ,1768 + ,0 + ,497 + ,0 + ,534 + ,0 + ,31/12/2007 + ,10667 + ,1 + ,3040 + ,3040 + ,3122 + ,3122 + ,1678 + ,1678 + ,1778 + ,1778 + ,506 + ,506 + ,543 + ,543 + ,31/12/2008 + ,10753 + ,2 + ,3064 + ,6128 + ,3145 + ,6290 + ,1687 + ,3374 + ,1789 + ,3578 + ,516 + ,1032 + ,553 + ,1106 + ,31/12/2009 + ,10840 + ,0 + ,3085 + ,0 + ,3167 + ,0 + ,1700 + ,0 + ,1798 + ,0 + ,527 + ,0 + ,563 + ,0 + ,31/12/2010 + ,10951 + ,1 + ,3113 + ,3113 + ,3193 + ,3193 + ,1714 + ,1714 + ,1811 + ,1811 + ,542 + ,542 + ,577 + ,577) + ,dim=c(15 + ,50) + ,dimnames=list(c('jaar' + ,'totaal' + ,'pop' + ,'vlaams_man' + ,'man_vlaams_pop' + ,'vlaams_vrouw' + ,'vrouw_vlaams_pop' + ,'waals_man' + ,'man_waals_pop' + ,'waals_vrouw' + ,'vrouw_waals_pop' + ,'brussel_man' + ,'man_brussel_pop' + ,'brussel_vrouw' + ,'vrouw_brussel_pop') + ,1:50)) > y <- array(NA,dim=c(15,50),dimnames=list(c('jaar','totaal','pop','vlaams_man','man_vlaams_pop','vlaams_vrouw','vrouw_vlaams_pop','waals_man','man_waals_pop','waals_vrouw','vrouw_waals_pop','brussel_man','man_brussel_pop','brussel_vrouw','vrouw_brussel_pop'),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 = 'Linear Trend' > par2 = 'Do not include Seasonal Dummies' > par1 = '2' > 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 totaal jaar pop vlaams_man man_vlaams_pop vlaams_vrouw 1 9190 0.001317355 0 2514 0 2550 2 9251 0.001316684 1 2537 2537 2572 3 9328 0.001316013 2 2564 5128 2597 4 9428 0.001315343 0 2595 0 2623 5 9499 0.001314673 1 2617 2617 2647 6 9556 0.001314005 2 2638 5276 2670 7 9606 0.001313337 0 2657 0 2690 8 9632 0.001312669 1 2668 2668 2705 9 9660 0.001312003 2 2683 5366 2721 10 9651 0.001311337 0 2687 0 2729 11 9695 0.001310671 1 2705 2705 2747 12 9727 0.001310007 2 2717 5434 2761 13 9757 0.001309343 0 2728 0 2773 14 9788 0.001308680 1 2741 2741 2786 15 9813 0.001308017 2 2752 5504 2796 16 9823 0.001307355 0 2759 0 2807 17 9837 0.001306694 1 2767 2767 2817 18 9842 0.001306033 2 2774 5548 2827 19 9855 0.001305373 0 2781 0 2838 20 9863 0.001304714 1 2788 2788 2847 21 9855 0.001304055 2 2789 5578 2853 22 9858 0.001303397 0 2795 0 2860 23 9853 0.001302740 1 2798 2798 2864 24 9858 0.001302083 2 2801 5602 2869 25 9859 0.001301427 0 2803 0 2873 26 9865 0.001300772 1 2808 2808 2877 27 9876 0.001300117 2 2813 5626 2883 28 9928 0.001299463 0 2826 0 2896 29 9948 0.001298810 1 2835 2835 2905 30 9987 0.001298157 2 2849 5698 2919 31 10022 0.001297505 0 2862 0 2933 32 10068 0.001296854 1 2877 2877 2948 33 10101 0.001296203 2 2888 5776 2959 34 10131 0.001295553 0 2897 0 2969 35 10143 0.001294904 1 2902 2902 2978 36 10170 0.001294255 2 2911 5822 2988 37 10192 0.001293607 0 2917 0 2996 38 10214 0.001292960 1 2924 2924 3003 39 10239 0.001292313 2 2930 5860 3011 40 10263 0.001291667 0 2935 0 3018 41 10310 0.001291021 1 2945 2945 3028 42 10355 0.001290376 2 2957 5914 3038 43 10396 0.001289732 0 2967 0 3049 44 10446 0.001289088 1 2980 2980 3063 45 10511 0.001288446 2 2997 5994 3081 46 10585 0.001287803 0 3017 0 3100 47 10667 0.001287162 1 3040 3040 3122 48 10753 0.001286521 2 3064 6128 3145 49 10840 0.001285880 0 3085 0 3167 50 10951 0.001285240 1 3113 3113 3193 vrouw_vlaams_pop waals_man man_waals_pop waals_vrouw vrouw_waals_pop 1 0 1512 0 1591 0 2 2572 1517 1517 1595 1595 3 5194 1525 3050 1602 3204 4 0 1540 0 1613 0 5 2647 1547 1547 1622 1622 6 5340 1547 3094 1627 3254 7 0 1547 0 1632 0 8 2705 1547 1547 1634 1634 9 5442 1546 3092 1637 3274 10 0 1533 0 1627 0 11 2747 1538 1538 1632 1632 12 5522 1543 3086 1637 3274 13 0 1549 0 1643 0 14 2786 1556 1556 1650 1650 15 5592 1559 3118 1654 3308 16 0 1559 0 1656 0 17 2817 1563 1563 1661 1661 18 5654 1563 3126 1662 3324 19 0 1564 0 1664 0 20 2847 1564 1564 1665 1665 21 5706 1557 3114 1661 3322 22 0 1554 0 1659 0 23 2864 1552 1552 1656 1656 24 5738 1552 3104 1656 3312 25 0 1551 0 1655 0 26 2877 1552 1552 1654 1654 27 5766 1554 3108 1656 3312 28 0 1567 0 1668 0 29 2905 1572 1572 1672 1672 30 5838 1579 3158 1680 3360 31 0 1588 0 1688 0 32 2948 1597 1597 1696 1696 33 5918 1603 3206 1702 3404 34 0 1607 0 1706 0 35 2978 1607 1607 1708 1708 36 5976 1609 3218 1711 3422 37 0 1612 0 1714 0 38 3003 1615 1615 1717 1717 39 6022 1619 3238 1721 3442 40 0 1622 0 1724 0 41 3028 1628 1628 1730 1730 42 6076 1634 3268 1735 3470 43 0 1640 0 1740 0 44 3063 1648 1648 1748 1748 45 6162 1657 3314 1757 3514 46 0 1668 0 1768 0 47 3122 1678 1678 1778 1778 48 6290 1687 3374 1789 3578 49 0 1700 0 1798 0 50 3193 1714 1714 1811 1811 brussel_man man_brussel_pop brussel_vrouw vrouw_brussel_pop t 1 472 0 551 0 1 2 476 476 554 554 2 3 483 966 558 1116 3 4 493 0 565 0 4 5 498 498 568 568 5 6 502 1004 572 1144 6 7 504 0 575 0 7 8 503 503 574 574 8 9 501 1002 572 1144 9 10 502 0 573 0 10 11 502 502 572 572 11 12 500 1000 569 1138 12 13 498 0 566 0 13 14 495 495 560 560 14 15 494 988 557 1114 15 16 490 0 552 0 16 17 484 484 545 545 17 18 477 954 539 1078 18 19 474 0 535 0 19 20 469 469 531 531 20 21 466 932 528 1056 21 22 464 0 526 0 22 23 460 460 523 523 23 24 458 916 521 1042 24 25 457 0 519 0 25 26 456 456 517 517 26 27 455 910 515 1030 27 28 456 0 514 0 28 29 453 453 511 511 29 30 453 906 508 1016 30 31 449 0 502 0 31 32 449 449 501 501 32 33 449 898 500 1000 33 34 452 0 500 0 34 35 450 450 498 498 35 36 452 904 499 998 36 37 454 0 499 0 37 38 455 455 500 500 38 39 458 916 501 1002 39 40 461 0 503 0 40 41 469 469 510 510 41 42 477 954 515 1030 42 43 480 0 520 0 43 44 484 484 523 523 44 45 490 980 529 1058 45 46 497 0 534 0 46 47 506 506 543 543 47 48 516 1032 553 1106 48 49 527 0 563 0 49 50 542 542 577 577 50 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) jaar pop vlaams_man 2.808e+04 -2.131e+07 -1.315e+01 8.365e-01 man_vlaams_pop vlaams_vrouw vrouw_vlaams_pop waals_man 4.715e-02 1.132e+00 -4.491e-02 1.112e+00 man_waals_pop waals_vrouw vrouw_waals_pop brussel_man -3.377e-02 9.287e-01 3.950e-02 9.511e-01 man_brussel_pop brussel_vrouw vrouw_brussel_pop t -2.260e-02 1.102e+00 1.389e-02 -1.382e+01 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -1.22528 -0.24954 0.03801 0.29416 1.21827 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.808e+04 1.792e+04 1.567 0.126 jaar -2.131e+07 1.361e+07 -1.567 0.126 pop -1.315e+01 2.977e+01 -0.442 0.661 vlaams_man 8.365e-01 1.061e-01 7.882 3.54e-09 *** man_vlaams_pop 4.715e-02 3.915e-02 1.204 0.237 vlaams_vrouw 1.132e+00 1.073e-01 10.552 2.89e-12 *** vrouw_vlaams_pop -4.491e-02 3.603e-02 -1.246 0.221 waals_man 1.112e+00 8.553e-02 13.003 9.46e-15 *** man_waals_pop -3.377e-02 4.728e-02 -0.714 0.480 waals_vrouw 9.287e-01 1.072e-01 8.660 4.05e-10 *** vrouw_waals_pop 3.950e-02 6.629e-02 0.596 0.555 brussel_man 9.511e-01 1.182e-01 8.044 2.24e-09 *** man_brussel_pop -2.260e-02 7.536e-02 -0.300 0.766 brussel_vrouw 1.102e+00 1.116e-01 9.870 1.63e-11 *** vrouw_brussel_pop 1.389e-02 6.772e-02 0.205 0.839 t -1.382e+01 8.781e+00 -1.574 0.125 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.6255 on 34 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 1.299e+06 on 15 and 34 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.2814685 0.562936995 0.718531502 [2,] 0.3065885 0.613176996 0.693411502 [3,] 0.1766513 0.353302670 0.823348665 [4,] 0.1764258 0.352851695 0.823574153 [5,] 0.4427707 0.885541493 0.557229253 [6,] 0.5894892 0.821021692 0.410510846 [7,] 0.8285228 0.342954325 0.171477163 [8,] 0.9509781 0.098043851 0.049021925 [9,] 0.9955757 0.008848642 0.004424321 [10,] 0.9971165 0.005767057 0.002883528 [11,] 0.9916005 0.016798998 0.008399499 [12,] 0.9834215 0.033157084 0.016578542 [13,] 0.9981110 0.003778070 0.001889035 > postscript(file="/var/wessaorg/rcomp/tmp/1d3741351885569.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/wessaorg/rcomp/tmp/2tnxr1351885569.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/wessaorg/rcomp/tmp/3iohe1351885569.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/wessaorg/rcomp/tmp/4yj3x1351885569.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/wessaorg/rcomp/tmp/547vh1351885569.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.03166215 0.27402781 -0.33394534 -0.84450862 0.10370564 0.14414127 7 8 9 10 11 12 1.21826717 0.75714246 -0.24842603 0.11180818 -0.93415407 -0.25016631 13 14 15 16 17 18 -0.18166714 -0.10043648 0.87787025 0.31729522 0.09747424 -0.24991109 19 20 21 22 23 24 -1.01210010 -1.22528282 0.65115862 -0.09679191 -0.22374878 0.64829122 25 26 27 28 29 30 0.38628299 0.66433178 -0.15155822 0.29822634 -0.51274420 -1.08478798 31 32 33 34 35 36 -0.03381595 0.05988518 0.14630376 0.28196940 0.25762758 0.51672333 37 38 39 40 41 42 0.16652147 0.43760915 -0.30159608 -0.17131536 0.08004213 -0.19611762 43 44 45 46 47 48 -0.49142259 -0.06356568 0.26960446 0.55001819 0.01613207 -0.62656388 49 50 -0.65608477 0.68991327 > postscript(file="/var/wessaorg/rcomp/tmp/6cwjp1351885569.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.03166215 NA 1 0.27402781 -0.03166215 2 -0.33394534 0.27402781 3 -0.84450862 -0.33394534 4 0.10370564 -0.84450862 5 0.14414127 0.10370564 6 1.21826717 0.14414127 7 0.75714246 1.21826717 8 -0.24842603 0.75714246 9 0.11180818 -0.24842603 10 -0.93415407 0.11180818 11 -0.25016631 -0.93415407 12 -0.18166714 -0.25016631 13 -0.10043648 -0.18166714 14 0.87787025 -0.10043648 15 0.31729522 0.87787025 16 0.09747424 0.31729522 17 -0.24991109 0.09747424 18 -1.01210010 -0.24991109 19 -1.22528282 -1.01210010 20 0.65115862 -1.22528282 21 -0.09679191 0.65115862 22 -0.22374878 -0.09679191 23 0.64829122 -0.22374878 24 0.38628299 0.64829122 25 0.66433178 0.38628299 26 -0.15155822 0.66433178 27 0.29822634 -0.15155822 28 -0.51274420 0.29822634 29 -1.08478798 -0.51274420 30 -0.03381595 -1.08478798 31 0.05988518 -0.03381595 32 0.14630376 0.05988518 33 0.28196940 0.14630376 34 0.25762758 0.28196940 35 0.51672333 0.25762758 36 0.16652147 0.51672333 37 0.43760915 0.16652147 38 -0.30159608 0.43760915 39 -0.17131536 -0.30159608 40 0.08004213 -0.17131536 41 -0.19611762 0.08004213 42 -0.49142259 -0.19611762 43 -0.06356568 -0.49142259 44 0.26960446 -0.06356568 45 0.55001819 0.26960446 46 0.01613207 0.55001819 47 -0.62656388 0.01613207 48 -0.65608477 -0.62656388 49 0.68991327 -0.65608477 50 NA 0.68991327 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] 0.27402781 -0.03166215 [2,] -0.33394534 0.27402781 [3,] -0.84450862 -0.33394534 [4,] 0.10370564 -0.84450862 [5,] 0.14414127 0.10370564 [6,] 1.21826717 0.14414127 [7,] 0.75714246 1.21826717 [8,] -0.24842603 0.75714246 [9,] 0.11180818 -0.24842603 [10,] -0.93415407 0.11180818 [11,] -0.25016631 -0.93415407 [12,] -0.18166714 -0.25016631 [13,] -0.10043648 -0.18166714 [14,] 0.87787025 -0.10043648 [15,] 0.31729522 0.87787025 [16,] 0.09747424 0.31729522 [17,] -0.24991109 0.09747424 [18,] -1.01210010 -0.24991109 [19,] -1.22528282 -1.01210010 [20,] 0.65115862 -1.22528282 [21,] -0.09679191 0.65115862 [22,] -0.22374878 -0.09679191 [23,] 0.64829122 -0.22374878 [24,] 0.38628299 0.64829122 [25,] 0.66433178 0.38628299 [26,] -0.15155822 0.66433178 [27,] 0.29822634 -0.15155822 [28,] -0.51274420 0.29822634 [29,] -1.08478798 -0.51274420 [30,] -0.03381595 -1.08478798 [31,] 0.05988518 -0.03381595 [32,] 0.14630376 0.05988518 [33,] 0.28196940 0.14630376 [34,] 0.25762758 0.28196940 [35,] 0.51672333 0.25762758 [36,] 0.16652147 0.51672333 [37,] 0.43760915 0.16652147 [38,] -0.30159608 0.43760915 [39,] -0.17131536 -0.30159608 [40,] 0.08004213 -0.17131536 [41,] -0.19611762 0.08004213 [42,] -0.49142259 -0.19611762 [43,] -0.06356568 -0.49142259 [44,] 0.26960446 -0.06356568 [45,] 0.55001819 0.26960446 [46,] 0.01613207 0.55001819 [47,] -0.62656388 0.01613207 [48,] -0.65608477 -0.62656388 [49,] 0.68991327 -0.65608477 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 0.27402781 -0.03166215 2 -0.33394534 0.27402781 3 -0.84450862 -0.33394534 4 0.10370564 -0.84450862 5 0.14414127 0.10370564 6 1.21826717 0.14414127 7 0.75714246 1.21826717 8 -0.24842603 0.75714246 9 0.11180818 -0.24842603 10 -0.93415407 0.11180818 11 -0.25016631 -0.93415407 12 -0.18166714 -0.25016631 13 -0.10043648 -0.18166714 14 0.87787025 -0.10043648 15 0.31729522 0.87787025 16 0.09747424 0.31729522 17 -0.24991109 0.09747424 18 -1.01210010 -0.24991109 19 -1.22528282 -1.01210010 20 0.65115862 -1.22528282 21 -0.09679191 0.65115862 22 -0.22374878 -0.09679191 23 0.64829122 -0.22374878 24 0.38628299 0.64829122 25 0.66433178 0.38628299 26 -0.15155822 0.66433178 27 0.29822634 -0.15155822 28 -0.51274420 0.29822634 29 -1.08478798 -0.51274420 30 -0.03381595 -1.08478798 31 0.05988518 -0.03381595 32 0.14630376 0.05988518 33 0.28196940 0.14630376 34 0.25762758 0.28196940 35 0.51672333 0.25762758 36 0.16652147 0.51672333 37 0.43760915 0.16652147 38 -0.30159608 0.43760915 39 -0.17131536 -0.30159608 40 0.08004213 -0.17131536 41 -0.19611762 0.08004213 42 -0.49142259 -0.19611762 43 -0.06356568 -0.49142259 44 0.26960446 -0.06356568 45 0.55001819 0.26960446 46 0.01613207 0.55001819 47 -0.62656388 0.01613207 48 -0.65608477 -0.62656388 49 0.68991327 -0.65608477 > 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/wessaorg/rcomp/tmp/75gph1351885569.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/wessaorg/rcomp/tmp/82o871351885569.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/wessaorg/rcomp/tmp/9ky021351885569.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/wessaorg/rcomp/tmp/103amt1351885569.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11fva11351885569.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/wessaorg/rcomp/tmp/12wuiu1351885569.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/wessaorg/rcomp/tmp/13eann1351885569.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/wessaorg/rcomp/tmp/14ak0o1351885569.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/wessaorg/rcomp/tmp/15i72a1351885569.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/wessaorg/rcomp/tmp/160w6x1351885569.tab") + } > > try(system("convert tmp/1d3741351885569.ps tmp/1d3741351885569.png",intern=TRUE)) character(0) > try(system("convert tmp/2tnxr1351885569.ps tmp/2tnxr1351885569.png",intern=TRUE)) character(0) > try(system("convert tmp/3iohe1351885569.ps tmp/3iohe1351885569.png",intern=TRUE)) character(0) > try(system("convert tmp/4yj3x1351885569.ps tmp/4yj3x1351885569.png",intern=TRUE)) character(0) > try(system("convert tmp/547vh1351885569.ps tmp/547vh1351885569.png",intern=TRUE)) character(0) > try(system("convert tmp/6cwjp1351885569.ps tmp/6cwjp1351885569.png",intern=TRUE)) character(0) > try(system("convert tmp/75gph1351885569.ps tmp/75gph1351885569.png",intern=TRUE)) character(0) > try(system("convert tmp/82o871351885569.ps tmp/82o871351885569.png",intern=TRUE)) character(0) > try(system("convert tmp/9ky021351885569.ps tmp/9ky021351885569.png",intern=TRUE)) character(0) > try(system("convert tmp/103amt1351885569.ps tmp/103amt1351885569.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 6.573 0.970 7.555