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