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