R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
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(695
+ ,0
+ ,641
+ ,696
+ ,729
+ ,839
+ ,627
+ ,608
+ ,638
+ ,0
+ ,695
+ ,641
+ ,696
+ ,729
+ ,696
+ ,651
+ ,762
+ ,0
+ ,638
+ ,695
+ ,641
+ ,696
+ ,825
+ ,691
+ ,635
+ ,0
+ ,762
+ ,638
+ ,695
+ ,641
+ ,677
+ ,627
+ ,721
+ ,0
+ ,635
+ ,762
+ ,638
+ ,695
+ ,656
+ ,634
+ ,854
+ ,0
+ ,721
+ ,635
+ ,762
+ ,638
+ ,785
+ ,731
+ ,418
+ ,0
+ ,854
+ ,721
+ ,635
+ ,762
+ ,412
+ ,475
+ ,367
+ ,0
+ ,418
+ ,854
+ ,721
+ ,635
+ ,352
+ ,337
+ ,824
+ ,0
+ ,367
+ ,418
+ ,854
+ ,721
+ ,839
+ ,803
+ ,687
+ ,0
+ ,824
+ ,367
+ ,418
+ ,854
+ ,729
+ ,722
+ ,601
+ ,0
+ ,687
+ ,824
+ ,367
+ ,418
+ ,696
+ ,590
+ ,676
+ ,0
+ ,601
+ ,687
+ ,824
+ ,367
+ ,641
+ ,724
+ ,740
+ ,0
+ ,676
+ ,601
+ ,687
+ ,824
+ ,695
+ ,627
+ ,691
+ ,0
+ ,740
+ ,676
+ ,601
+ ,687
+ ,638
+ ,696
+ ,683
+ ,0
+ ,691
+ ,740
+ ,676
+ ,601
+ ,762
+ ,825
+ ,594
+ ,0
+ ,683
+ ,691
+ ,740
+ ,676
+ ,635
+ ,677
+ ,729
+ ,0
+ ,594
+ ,683
+ ,691
+ ,740
+ ,721
+ ,656
+ ,731
+ ,0
+ ,729
+ ,594
+ ,683
+ ,691
+ ,854
+ ,785
+ ,386
+ ,0
+ ,731
+ ,729
+ ,594
+ ,683
+ ,418
+ ,412
+ ,331
+ ,0
+ ,386
+ ,731
+ ,729
+ ,594
+ ,367
+ ,352
+ ,706
+ ,0
+ ,331
+ ,386
+ ,731
+ ,729
+ ,824
+ ,839
+ ,715
+ ,0
+ ,706
+ ,331
+ ,386
+ ,731
+ ,687
+ ,729
+ ,657
+ ,0
+ ,715
+ ,706
+ ,331
+ ,386
+ ,601
+ ,696
+ ,653
+ ,0
+ ,657
+ ,715
+ ,706
+ ,331
+ ,676
+ ,641
+ ,642
+ ,0
+ ,653
+ ,657
+ ,715
+ ,706
+ ,740
+ ,695
+ ,643
+ ,0
+ ,642
+ ,653
+ ,657
+ ,715
+ ,691
+ ,638
+ ,718
+ ,0
+ ,643
+ ,642
+ ,653
+ ,657
+ ,683
+ ,762
+ ,654
+ ,0
+ ,718
+ ,643
+ ,642
+ ,653
+ ,594
+ ,635
+ ,632
+ ,0
+ ,654
+ ,718
+ ,643
+ ,642
+ ,729
+ ,721
+ ,731
+ ,0
+ ,632
+ ,654
+ ,718
+ ,643
+ ,731
+ ,854
+ ,392
+ ,1
+ ,731
+ ,632
+ ,654
+ ,718
+ ,386
+ ,418
+ ,344
+ ,1
+ ,392
+ ,731
+ ,632
+ ,654
+ ,331
+ ,367
+ ,792
+ ,1
+ ,344
+ ,392
+ ,731
+ ,632
+ ,706
+ ,824
+ ,852
+ ,1
+ ,792
+ ,344
+ ,392
+ ,731
+ ,715
+ ,687
+ ,649
+ ,1
+ ,852
+ ,792
+ ,344
+ ,392
+ ,657
+ ,601
+ ,629
+ ,1
+ ,649
+ ,852
+ ,792
+ ,344
+ ,653
+ ,676
+ ,685
+ ,1
+ ,629
+ ,649
+ ,852
+ ,792
+ ,642
+ ,740
+ ,617
+ ,1
+ ,685
+ ,629
+ ,649
+ ,852
+ ,643
+ ,691
+ ,715
+ ,1
+ ,617
+ ,685
+ ,629
+ ,649
+ ,718
+ ,683
+ ,715
+ ,1
+ ,715
+ ,617
+ ,685
+ ,629
+ ,654
+ ,594
+ ,629
+ ,1
+ ,715
+ ,715
+ ,617
+ ,685
+ ,632
+ ,729
+ ,916
+ ,1
+ ,629
+ ,715
+ ,715
+ ,617
+ ,731
+ ,731
+ ,531
+ ,1
+ ,916
+ ,629
+ ,715
+ ,715
+ ,392
+ ,386
+ ,357
+ ,1
+ ,531
+ ,916
+ ,629
+ ,715
+ ,344
+ ,331
+ ,917
+ ,1
+ ,357
+ ,531
+ ,916
+ ,629
+ ,792
+ ,706
+ ,828
+ ,1
+ ,917
+ ,357
+ ,531
+ ,916
+ ,852
+ ,715
+ ,708
+ ,1
+ ,828
+ ,917
+ ,357
+ ,531
+ ,649
+ ,657
+ ,858
+ ,1
+ ,708
+ ,828
+ ,917
+ ,357
+ ,629
+ ,653
+ ,775
+ ,1
+ ,858
+ ,708
+ ,828
+ ,917
+ ,685
+ ,642
+ ,785
+ ,1
+ ,775
+ ,858
+ ,708
+ ,828
+ ,617
+ ,643
+ ,1006
+ ,1
+ ,785
+ ,775
+ ,858
+ ,708
+ ,715
+ ,718
+ ,789
+ ,1
+ ,1006
+ ,785
+ ,775
+ ,858
+ ,715
+ ,654
+ ,734
+ ,1
+ ,789
+ ,1006
+ ,785
+ ,775
+ ,629
+ ,632
+ ,906
+ ,1
+ ,734
+ ,789
+ ,1006
+ ,785
+ ,916
+ ,731
+ ,532
+ ,1
+ ,906
+ ,734
+ ,789
+ ,1006
+ ,531
+ ,392
+ ,387
+ ,1
+ ,532
+ ,906
+ ,734
+ ,789
+ ,357
+ ,344
+ ,991
+ ,1
+ ,387
+ ,532
+ ,906
+ ,734
+ ,917
+ ,792
+ ,841
+ ,1
+ ,991
+ ,387
+ ,532
+ ,906
+ ,828
+ ,852
+ ,892
+ ,1
+ ,841
+ ,991
+ ,387
+ ,532
+ ,708
+ ,649
+ ,782
+ ,1
+ ,892
+ ,841
+ ,991
+ ,387
+ ,858
+ ,629)
+ ,dim=c(8
+ ,60)
+ ,dimnames=list(c('faillissement'
+ ,'crisis'
+ ,'t-1'
+ ,'t-2'
+ ,'t-3'
+ ,'t-4'
+ ,'t-12'
+ ,'t-24')
+ ,1:60))
> y <- array(NA,dim=c(8,60),dimnames=list(c('faillissement','crisis','t-1','t-2','t-3','t-4','t-12','t-24'),1:60))
> 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'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
faillissement crisis t-1 t-2 t-3 t-4 t-12 t-24
1 695 0 641 696 729 839 627 608
2 638 0 695 641 696 729 696 651
3 762 0 638 695 641 696 825 691
4 635 0 762 638 695 641 677 627
5 721 0 635 762 638 695 656 634
6 854 0 721 635 762 638 785 731
7 418 0 854 721 635 762 412 475
8 367 0 418 854 721 635 352 337
9 824 0 367 418 854 721 839 803
10 687 0 824 367 418 854 729 722
11 601 0 687 824 367 418 696 590
12 676 0 601 687 824 367 641 724
13 740 0 676 601 687 824 695 627
14 691 0 740 676 601 687 638 696
15 683 0 691 740 676 601 762 825
16 594 0 683 691 740 676 635 677
17 729 0 594 683 691 740 721 656
18 731 0 729 594 683 691 854 785
19 386 0 731 729 594 683 418 412
20 331 0 386 731 729 594 367 352
21 706 0 331 386 731 729 824 839
22 715 0 706 331 386 731 687 729
23 657 0 715 706 331 386 601 696
24 653 0 657 715 706 331 676 641
25 642 0 653 657 715 706 740 695
26 643 0 642 653 657 715 691 638
27 718 0 643 642 653 657 683 762
28 654 0 718 643 642 653 594 635
29 632 0 654 718 643 642 729 721
30 731 0 632 654 718 643 731 854
31 392 1 731 632 654 718 386 418
32 344 1 392 731 632 654 331 367
33 792 1 344 392 731 632 706 824
34 852 1 792 344 392 731 715 687
35 649 1 852 792 344 392 657 601
36 629 1 649 852 792 344 653 676
37 685 1 629 649 852 792 642 740
38 617 1 685 629 649 852 643 691
39 715 1 617 685 629 649 718 683
40 715 1 715 617 685 629 654 594
41 629 1 715 715 617 685 632 729
42 916 1 629 715 715 617 731 731
43 531 1 916 629 715 715 392 386
44 357 1 531 916 629 715 344 331
45 917 1 357 531 916 629 792 706
46 828 1 917 357 531 916 852 715
47 708 1 828 917 357 531 649 657
48 858 1 708 828 917 357 629 653
49 775 1 858 708 828 917 685 642
50 785 1 775 858 708 828 617 643
51 1006 1 785 775 858 708 715 718
52 789 1 1006 785 775 858 715 654
53 734 1 789 1006 785 775 629 632
54 906 1 734 789 1006 785 916 731
55 532 1 906 734 789 1006 531 392
56 387 1 532 906 734 789 357 344
57 991 1 387 532 906 734 917 792
58 841 1 991 387 532 906 828 852
59 892 1 841 991 387 532 708 649
60 782 1 892 841 991 387 858 629
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) crisis `t-1` `t-2` `t-3` `t-4`
-1.227e+02 7.846e+01 5.238e-02 4.446e-02 6.345e-02 -7.502e-03
`t-12` `t-24`
6.689e-01 3.456e-01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-116.658 -48.727 -7.443 50.785 199.203
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.227e+02 1.186e+02 -1.034 0.305735
crisis 7.846e+01 2.048e+01 3.832 0.000344 ***
`t-1` 5.238e-02 7.240e-02 0.723 0.472633
`t-2` 4.446e-02 8.217e-02 0.541 0.590807
`t-3` 6.345e-02 7.297e-02 0.869 0.388572
`t-4` -7.502e-03 7.504e-02 -0.100 0.920745
`t-12` 6.689e-01 1.432e-01 4.671 2.15e-05 ***
`t-24` 3.456e-01 1.540e-01 2.244 0.029146 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 73.88 on 52 degrees of freedom
Multiple R-squared: 0.81, Adjusted R-squared: 0.7844
F-statistic: 31.66 on 7 and 52 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.36629683 0.73259366 0.6337032
[2,] 0.23551812 0.47103624 0.7644819
[3,] 0.20249549 0.40499097 0.7975045
[4,] 0.11822235 0.23644471 0.8817776
[5,] 0.30626662 0.61253323 0.6937334
[6,] 0.27030215 0.54060431 0.7296978
[7,] 0.19340368 0.38680736 0.8065963
[8,] 0.23894033 0.47788066 0.7610597
[9,] 0.18082212 0.36164423 0.8191779
[10,] 0.14058047 0.28116094 0.8594195
[11,] 0.13692774 0.27385548 0.8630723
[12,] 0.16341215 0.32682430 0.8365879
[13,] 0.16314844 0.32629689 0.8368516
[14,] 0.11299041 0.22598083 0.8870096
[15,] 0.10634885 0.21269769 0.8936512
[16,] 0.07175346 0.14350692 0.9282465
[17,] 0.04872630 0.09745260 0.9512737
[18,] 0.04593059 0.09186119 0.9540694
[19,] 0.04173304 0.08346608 0.9582670
[20,] 0.02585658 0.05171317 0.9741434
[21,] 0.01633876 0.03267752 0.9836612
[22,] 0.01003576 0.02007152 0.9899642
[23,] 0.00843734 0.01687468 0.9915627
[24,] 0.01629443 0.03258886 0.9837056
[25,] 0.01024478 0.02048956 0.9897552
[26,] 0.01733295 0.03466590 0.9826671
[27,] 0.01568979 0.03137957 0.9843102
[28,] 0.02334075 0.04668150 0.9766592
[29,] 0.01890207 0.03780414 0.9810979
[30,] 0.01198951 0.02397903 0.9880105
[31,] 0.07650263 0.15300526 0.9234974
[32,] 0.13511061 0.27022122 0.8648894
[33,] 0.12981659 0.25963317 0.8701834
[34,] 0.10191041 0.20382083 0.8980896
[35,] 0.08266987 0.16533974 0.9173301
[36,] 0.06395201 0.12790402 0.9360480
[37,] 0.06346432 0.12692865 0.9365357
[38,] 0.05213773 0.10427545 0.9478623
[39,] 0.02453573 0.04907145 0.9754643
> postscript(file="/var/www/html/freestat/rcomp/tmp/19fz61292773340.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/2jog91292773340.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3jog91292773340.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4jog91292773340.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5ugfc1292773340.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 = 60
Frequency = 1
1 2 3 4 5 6
83.758830 -33.365883 -5.643042 -19.337167 83.451061 89.496384
7 8 9 10 11 12
-10.361028 36.973105 21.473527 -6.972854 -38.462714 8.237041
13 14 15 16 17 18
81.653576 44.676563 -96.520253 -50.333969 43.008045 -91.500986
19 20 21 22 23 24
-16.508207 -7.913894 -87.761622 55.588808 50.272162 -6.452926
25 26 27 28 29 30
-73.889264 -15.917318 21.839546 57.948148 -84.192660 -33.243010
31 32 33 34 35 36
-68.873711 -48.191978 2.207338 104.447629 -52.598082 -116.657682
37 38 39 40 41 42
-65.789585 -106.240559 -54.821658 12.926260 -104.631665 113.237781
43 44 45 46 47 48
63.721336 -46.304269 91.841603 -35.417215 -11.680455 126.485873
49 50 51 52 53 54
17.156986 76.918266 199.202854 -1.310675 9.097358 -46.491865
55 56 57 58 59 60
-36.978336 -35.206089 52.323488 -59.054592 129.755187 -109.073541
> postscript(file="/var/www/html/freestat/rcomp/tmp/6ugfc1292773340.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 83.758830 NA
1 -33.365883 83.758830
2 -5.643042 -33.365883
3 -19.337167 -5.643042
4 83.451061 -19.337167
5 89.496384 83.451061
6 -10.361028 89.496384
7 36.973105 -10.361028
8 21.473527 36.973105
9 -6.972854 21.473527
10 -38.462714 -6.972854
11 8.237041 -38.462714
12 81.653576 8.237041
13 44.676563 81.653576
14 -96.520253 44.676563
15 -50.333969 -96.520253
16 43.008045 -50.333969
17 -91.500986 43.008045
18 -16.508207 -91.500986
19 -7.913894 -16.508207
20 -87.761622 -7.913894
21 55.588808 -87.761622
22 50.272162 55.588808
23 -6.452926 50.272162
24 -73.889264 -6.452926
25 -15.917318 -73.889264
26 21.839546 -15.917318
27 57.948148 21.839546
28 -84.192660 57.948148
29 -33.243010 -84.192660
30 -68.873711 -33.243010
31 -48.191978 -68.873711
32 2.207338 -48.191978
33 104.447629 2.207338
34 -52.598082 104.447629
35 -116.657682 -52.598082
36 -65.789585 -116.657682
37 -106.240559 -65.789585
38 -54.821658 -106.240559
39 12.926260 -54.821658
40 -104.631665 12.926260
41 113.237781 -104.631665
42 63.721336 113.237781
43 -46.304269 63.721336
44 91.841603 -46.304269
45 -35.417215 91.841603
46 -11.680455 -35.417215
47 126.485873 -11.680455
48 17.156986 126.485873
49 76.918266 17.156986
50 199.202854 76.918266
51 -1.310675 199.202854
52 9.097358 -1.310675
53 -46.491865 9.097358
54 -36.978336 -46.491865
55 -35.206089 -36.978336
56 52.323488 -35.206089
57 -59.054592 52.323488
58 129.755187 -59.054592
59 -109.073541 129.755187
60 NA -109.073541
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -33.365883 83.758830
[2,] -5.643042 -33.365883
[3,] -19.337167 -5.643042
[4,] 83.451061 -19.337167
[5,] 89.496384 83.451061
[6,] -10.361028 89.496384
[7,] 36.973105 -10.361028
[8,] 21.473527 36.973105
[9,] -6.972854 21.473527
[10,] -38.462714 -6.972854
[11,] 8.237041 -38.462714
[12,] 81.653576 8.237041
[13,] 44.676563 81.653576
[14,] -96.520253 44.676563
[15,] -50.333969 -96.520253
[16,] 43.008045 -50.333969
[17,] -91.500986 43.008045
[18,] -16.508207 -91.500986
[19,] -7.913894 -16.508207
[20,] -87.761622 -7.913894
[21,] 55.588808 -87.761622
[22,] 50.272162 55.588808
[23,] -6.452926 50.272162
[24,] -73.889264 -6.452926
[25,] -15.917318 -73.889264
[26,] 21.839546 -15.917318
[27,] 57.948148 21.839546
[28,] -84.192660 57.948148
[29,] -33.243010 -84.192660
[30,] -68.873711 -33.243010
[31,] -48.191978 -68.873711
[32,] 2.207338 -48.191978
[33,] 104.447629 2.207338
[34,] -52.598082 104.447629
[35,] -116.657682 -52.598082
[36,] -65.789585 -116.657682
[37,] -106.240559 -65.789585
[38,] -54.821658 -106.240559
[39,] 12.926260 -54.821658
[40,] -104.631665 12.926260
[41,] 113.237781 -104.631665
[42,] 63.721336 113.237781
[43,] -46.304269 63.721336
[44,] 91.841603 -46.304269
[45,] -35.417215 91.841603
[46,] -11.680455 -35.417215
[47,] 126.485873 -11.680455
[48,] 17.156986 126.485873
[49,] 76.918266 17.156986
[50,] 199.202854 76.918266
[51,] -1.310675 199.202854
[52,] 9.097358 -1.310675
[53,] -46.491865 9.097358
[54,] -36.978336 -46.491865
[55,] -35.206089 -36.978336
[56,] 52.323488 -35.206089
[57,] -59.054592 52.323488
[58,] 129.755187 -59.054592
[59,] -109.073541 129.755187
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -33.365883 83.758830
2 -5.643042 -33.365883
3 -19.337167 -5.643042
4 83.451061 -19.337167
5 89.496384 83.451061
6 -10.361028 89.496384
7 36.973105 -10.361028
8 21.473527 36.973105
9 -6.972854 21.473527
10 -38.462714 -6.972854
11 8.237041 -38.462714
12 81.653576 8.237041
13 44.676563 81.653576
14 -96.520253 44.676563
15 -50.333969 -96.520253
16 43.008045 -50.333969
17 -91.500986 43.008045
18 -16.508207 -91.500986
19 -7.913894 -16.508207
20 -87.761622 -7.913894
21 55.588808 -87.761622
22 50.272162 55.588808
23 -6.452926 50.272162
24 -73.889264 -6.452926
25 -15.917318 -73.889264
26 21.839546 -15.917318
27 57.948148 21.839546
28 -84.192660 57.948148
29 -33.243010 -84.192660
30 -68.873711 -33.243010
31 -48.191978 -68.873711
32 2.207338 -48.191978
33 104.447629 2.207338
34 -52.598082 104.447629
35 -116.657682 -52.598082
36 -65.789585 -116.657682
37 -106.240559 -65.789585
38 -54.821658 -106.240559
39 12.926260 -54.821658
40 -104.631665 12.926260
41 113.237781 -104.631665
42 63.721336 113.237781
43 -46.304269 63.721336
44 91.841603 -46.304269
45 -35.417215 91.841603
46 -11.680455 -35.417215
47 126.485873 -11.680455
48 17.156986 126.485873
49 76.918266 17.156986
50 199.202854 76.918266
51 -1.310675 199.202854
52 9.097358 -1.310675
53 -46.491865 9.097358
54 -36.978336 -46.491865
55 -35.206089 -36.978336
56 52.323488 -35.206089
57 -59.054592 52.323488
58 129.755187 -59.054592
59 -109.073541 129.755187
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/75pwf1292773340.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/85pwf1292773340.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/9fyei1292773340.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10fyei1292773340.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/111hc51292773340.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/124zbt1292773340.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/1309921292773340.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/14mrp81292773340.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/157ane1292773340.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/16as4k1292773340.tab")
+ }
>
> try(system("convert tmp/19fz61292773340.ps tmp/19fz61292773340.png",intern=TRUE))
character(0)
> try(system("convert tmp/2jog91292773340.ps tmp/2jog91292773340.png",intern=TRUE))
character(0)
> try(system("convert tmp/3jog91292773340.ps tmp/3jog91292773340.png",intern=TRUE))
character(0)
> try(system("convert tmp/4jog91292773340.ps tmp/4jog91292773340.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ugfc1292773340.ps tmp/5ugfc1292773340.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ugfc1292773340.ps tmp/6ugfc1292773340.png",intern=TRUE))
character(0)
> try(system("convert tmp/75pwf1292773340.ps tmp/75pwf1292773340.png",intern=TRUE))
character(0)
> try(system("convert tmp/85pwf1292773340.ps tmp/85pwf1292773340.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fyei1292773340.ps tmp/9fyei1292773340.png",intern=TRUE))
character(0)
> try(system("convert tmp/10fyei1292773340.ps tmp/10fyei1292773340.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.968 2.529 4.380