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(627
+ ,0
+ ,724
+ ,590
+ ,722
+ ,803
+ ,608
+ ,696
+ ,0
+ ,627
+ ,724
+ ,590
+ ,722
+ ,651
+ ,825
+ ,0
+ ,696
+ ,627
+ ,724
+ ,590
+ ,691
+ ,677
+ ,0
+ ,825
+ ,696
+ ,627
+ ,724
+ ,627
+ ,656
+ ,0
+ ,677
+ ,825
+ ,696
+ ,627
+ ,634
+ ,785
+ ,0
+ ,656
+ ,677
+ ,825
+ ,696
+ ,731
+ ,412
+ ,0
+ ,785
+ ,656
+ ,677
+ ,825
+ ,475
+ ,352
+ ,0
+ ,412
+ ,785
+ ,656
+ ,677
+ ,337
+ ,839
+ ,0
+ ,352
+ ,412
+ ,785
+ ,656
+ ,803
+ ,729
+ ,0
+ ,839
+ ,352
+ ,412
+ ,785
+ ,722
+ ,696
+ ,0
+ ,729
+ ,839
+ ,352
+ ,412
+ ,590
+ ,641
+ ,0
+ ,696
+ ,729
+ ,839
+ ,352
+ ,724
+ ,695
+ ,0
+ ,641
+ ,696
+ ,729
+ ,839
+ ,627
+ ,638
+ ,0
+ ,695
+ ,641
+ ,696
+ ,729
+ ,696
+ ,762
+ ,0
+ ,638
+ ,695
+ ,641
+ ,696
+ ,825
+ ,635
+ ,0
+ ,762
+ ,638
+ ,695
+ ,641
+ ,677
+ ,721
+ ,0
+ ,635
+ ,762
+ ,638
+ ,695
+ ,656
+ ,854
+ ,0
+ ,721
+ ,635
+ ,762
+ ,638
+ ,785
+ ,418
+ ,0
+ ,854
+ ,721
+ ,635
+ ,762
+ ,412
+ ,367
+ ,0
+ ,418
+ ,854
+ ,721
+ ,635
+ ,352
+ ,824
+ ,0
+ ,367
+ ,418
+ ,854
+ ,721
+ ,839
+ ,687
+ ,0
+ ,824
+ ,367
+ ,418
+ ,854
+ ,729
+ ,601
+ ,0
+ ,687
+ ,824
+ ,367
+ ,418
+ ,696
+ ,676
+ ,0
+ ,601
+ ,687
+ ,824
+ ,367
+ ,641
+ ,740
+ ,0
+ ,676
+ ,601
+ ,687
+ ,824
+ ,695
+ ,691
+ ,0
+ ,740
+ ,676
+ ,601
+ ,687
+ ,638
+ ,683
+ ,0
+ ,691
+ ,740
+ ,676
+ ,601
+ ,762
+ ,594
+ ,0
+ ,683
+ ,691
+ ,740
+ ,676
+ ,635
+ ,729
+ ,0
+ ,594
+ ,683
+ ,691
+ ,740
+ ,721
+ ,731
+ ,0
+ ,729
+ ,594
+ ,683
+ ,691
+ ,854
+ ,386
+ ,0
+ ,731
+ ,729
+ ,594
+ ,683
+ ,418
+ ,331
+ ,0
+ ,386
+ ,731
+ ,729
+ ,594
+ ,367
+ ,706
+ ,0
+ ,331
+ ,386
+ ,731
+ ,729
+ ,824
+ ,715
+ ,0
+ ,706
+ ,331
+ ,386
+ ,731
+ ,687
+ ,657
+ ,0
+ ,715
+ ,706
+ ,331
+ ,386
+ ,601
+ ,653
+ ,0
+ ,657
+ ,715
+ ,706
+ ,331
+ ,676
+ ,642
+ ,0
+ ,653
+ ,657
+ ,715
+ ,706
+ ,740
+ ,643
+ ,0
+ ,642
+ ,653
+ ,657
+ ,715
+ ,691
+ ,718
+ ,0
+ ,643
+ ,642
+ ,653
+ ,657
+ ,683
+ ,654
+ ,0
+ ,718
+ ,643
+ ,642
+ ,653
+ ,594
+ ,632
+ ,0
+ ,654
+ ,718
+ ,643
+ ,642
+ ,729
+ ,731
+ ,0
+ ,632
+ ,654
+ ,718
+ ,643
+ ,731
+ ,392
+ ,1
+ ,731
+ ,632
+ ,654
+ ,718
+ ,386
+ ,344
+ ,1
+ ,392
+ ,731
+ ,632
+ ,654
+ ,331
+ ,792
+ ,1
+ ,344
+ ,392
+ ,731
+ ,632
+ ,706
+ ,852
+ ,1
+ ,792
+ ,344
+ ,392
+ ,731
+ ,715
+ ,649
+ ,1
+ ,852
+ ,792
+ ,344
+ ,392
+ ,657
+ ,629
+ ,1
+ ,649
+ ,852
+ ,792
+ ,344
+ ,653
+ ,685
+ ,1
+ ,629
+ ,649
+ ,852
+ ,792
+ ,642
+ ,617
+ ,1
+ ,685
+ ,629
+ ,649
+ ,852
+ ,643
+ ,715
+ ,1
+ ,617
+ ,685
+ ,629
+ ,649
+ ,718
+ ,715
+ ,1
+ ,715
+ ,617
+ ,685
+ ,629
+ ,654
+ ,629
+ ,1
+ ,715
+ ,715
+ ,617
+ ,685
+ ,632
+ ,916
+ ,1
+ ,629
+ ,715
+ ,715
+ ,617
+ ,731
+ ,531
+ ,1
+ ,916
+ ,629
+ ,715
+ ,715
+ ,392
+ ,357
+ ,1
+ ,531
+ ,916
+ ,629
+ ,715
+ ,344
+ ,917
+ ,1
+ ,357
+ ,531
+ ,916
+ ,629
+ ,792
+ ,828
+ ,1
+ ,917
+ ,357
+ ,531
+ ,916
+ ,852
+ ,708
+ ,1
+ ,828
+ ,917
+ ,357
+ ,531
+ ,649
+ ,858
+ ,1
+ ,708
+ ,828
+ ,917
+ ,357
+ ,629
+ ,775
+ ,1
+ ,858
+ ,708
+ ,828
+ ,917
+ ,685
+ ,785
+ ,1
+ ,775
+ ,858
+ ,708
+ ,828
+ ,617
+ ,1.006
+ ,1
+ ,785
+ ,775
+ ,858
+ ,708
+ ,715
+ ,789
+ ,1
+ ,1006
+ ,785
+ ,775
+ ,858
+ ,715
+ ,734
+ ,1
+ ,789
+ ,1006
+ ,785
+ ,775
+ ,629
+ ,906
+ ,1
+ ,734
+ ,789
+ ,1006
+ ,785
+ ,916
+ ,532
+ ,1
+ ,906
+ ,734
+ ,789
+ ,1006
+ ,531
+ ,387
+ ,1
+ ,532
+ ,906
+ ,734
+ ,789
+ ,357
+ ,991
+ ,1
+ ,387
+ ,532
+ ,906
+ ,734
+ ,917
+ ,841
+ ,1
+ ,991
+ ,387
+ ,532
+ ,906
+ ,828
+ ,892
+ ,1
+ ,841
+ ,991
+ ,387
+ ,532
+ ,708
+ ,782
+ ,1
+ ,892
+ ,841
+ ,991
+ ,387
+ ,858)
+ ,dim=c(7
+ ,72)
+ ,dimnames=list(c('faillissement'
+ ,'crisis'
+ ,'t-1'
+ ,'t-2'
+ ,'t-3'
+ ,'t-4'
+ ,'t-12')
+ ,1:72))
> y <- array(NA,dim=c(7,72),dimnames=list(c('faillissement','crisis','t-1','t-2','t-3','t-4','t-12'),1:72))
> 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
1 627.000 0 724 590 722 803 608
2 696.000 0 627 724 590 722 651
3 825.000 0 696 627 724 590 691
4 677.000 0 825 696 627 724 627
5 656.000 0 677 825 696 627 634
6 785.000 0 656 677 825 696 731
7 412.000 0 785 656 677 825 475
8 352.000 0 412 785 656 677 337
9 839.000 0 352 412 785 656 803
10 729.000 0 839 352 412 785 722
11 696.000 0 729 839 352 412 590
12 641.000 0 696 729 839 352 724
13 695.000 0 641 696 729 839 627
14 638.000 0 695 641 696 729 696
15 762.000 0 638 695 641 696 825
16 635.000 0 762 638 695 641 677
17 721.000 0 635 762 638 695 656
18 854.000 0 721 635 762 638 785
19 418.000 0 854 721 635 762 412
20 367.000 0 418 854 721 635 352
21 824.000 0 367 418 854 721 839
22 687.000 0 824 367 418 854 729
23 601.000 0 687 824 367 418 696
24 676.000 0 601 687 824 367 641
25 740.000 0 676 601 687 824 695
26 691.000 0 740 676 601 687 638
27 683.000 0 691 740 676 601 762
28 594.000 0 683 691 740 676 635
29 729.000 0 594 683 691 740 721
30 731.000 0 729 594 683 691 854
31 386.000 0 731 729 594 683 418
32 331.000 0 386 731 729 594 367
33 706.000 0 331 386 731 729 824
34 715.000 0 706 331 386 731 687
35 657.000 0 715 706 331 386 601
36 653.000 0 657 715 706 331 676
37 642.000 0 653 657 715 706 740
38 643.000 0 642 653 657 715 691
39 718.000 0 643 642 653 657 683
40 654.000 0 718 643 642 653 594
41 632.000 0 654 718 643 642 729
42 731.000 0 632 654 718 643 731
43 392.000 1 731 632 654 718 386
44 344.000 1 392 731 632 654 331
45 792.000 1 344 392 731 632 706
46 852.000 1 792 344 392 731 715
47 649.000 1 852 792 344 392 657
48 629.000 1 649 852 792 344 653
49 685.000 1 629 649 852 792 642
50 617.000 1 685 629 649 852 643
51 715.000 1 617 685 629 649 718
52 715.000 1 715 617 685 629 654
53 629.000 1 715 715 617 685 632
54 916.000 1 629 715 715 617 731
55 531.000 1 916 629 715 715 392
56 357.000 1 531 916 629 715 344
57 917.000 1 357 531 916 629 792
58 828.000 1 917 357 531 916 852
59 708.000 1 828 917 357 531 649
60 858.000 1 708 828 917 357 629
61 775.000 1 858 708 828 917 685
62 785.000 1 775 858 708 828 617
63 1.006 1 785 775 858 708 715
64 789.000 1 1006 785 775 858 715
65 734.000 1 789 1006 785 775 629
66 906.000 1 734 789 1006 785 916
67 532.000 1 906 734 789 1006 531
68 387.000 1 532 906 734 789 357
69 991.000 1 387 532 906 734 917
70 841.000 1 991 387 532 906 828
71 892.000 1 841 991 387 532 708
72 782.000 1 892 841 991 387 858
> 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`
179.31295 44.39697 -0.04541 -0.02566 -0.09497 -0.03389
`t-12`
0.92150
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-720.566 -48.570 3.459 59.523 207.248
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 179.31295 160.67113 1.116 0.269
crisis 44.39697 29.45785 1.507 0.137
`t-1` -0.04541 0.10377 -0.438 0.663
`t-2` -0.02566 0.11718 -0.219 0.827
`t-3` -0.09497 0.10423 -0.911 0.366
`t-4` -0.03389 0.10826 -0.313 0.755
`t-12` 0.92150 0.11662 7.902 4.36e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 115.6 on 65 degrees of freedom
Multiple R-squared: 0.5825, Adjusted R-squared: 0.544
F-statistic: 15.11 on 6 and 65 DF, p-value: 9.602e-11
> 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,] 7.350324e-02 1.470065e-01 0.9264968
[2,] 3.001349e-02 6.002698e-02 0.9699865
[3,] 1.161643e-01 2.323286e-01 0.8838357
[4,] 5.944055e-02 1.188811e-01 0.9405595
[5,] 6.178142e-02 1.235628e-01 0.9382186
[6,] 8.363758e-02 1.672752e-01 0.9163624
[7,] 5.170780e-02 1.034156e-01 0.9482922
[8,] 3.029782e-02 6.059563e-02 0.9697022
[9,] 2.227477e-02 4.454955e-02 0.9777252
[10,] 1.139488e-02 2.278975e-02 0.9886051
[11,] 5.624393e-03 1.124879e-02 0.9943756
[12,] 2.773232e-03 5.546465e-03 0.9972268
[13,] 1.679824e-03 3.359648e-03 0.9983202
[14,] 3.149184e-03 6.298368e-03 0.9968508
[15,] 1.673361e-03 3.346721e-03 0.9983266
[16,] 8.900830e-04 1.780166e-03 0.9991099
[17,] 4.938550e-04 9.877099e-04 0.9995061
[18,] 4.362226e-04 8.724451e-04 0.9995638
[19,] 2.475878e-04 4.951756e-04 0.9997524
[20,] 1.121234e-04 2.242468e-04 0.9998879
[21,] 1.468576e-04 2.937152e-04 0.9998531
[22,] 8.105464e-05 1.621093e-04 0.9999189
[23,] 4.886717e-05 9.773434e-05 0.9999511
[24,] 5.134104e-05 1.026821e-04 0.9999487
[25,] 2.599821e-05 5.199642e-05 0.9999740
[26,] 1.322827e-05 2.645654e-05 0.9999868
[27,] 5.588153e-06 1.117631e-05 0.9999944
[28,] 4.568305e-06 9.136611e-06 0.9999954
[29,] 2.122823e-06 4.245646e-06 0.9999979
[30,] 9.541866e-07 1.908373e-06 0.9999990
[31,] 5.113581e-07 1.022716e-06 0.9999995
[32,] 3.855242e-07 7.710484e-07 0.9999996
[33,] 1.419046e-07 2.838091e-07 0.9999999
[34,] 5.138926e-08 1.027785e-07 0.9999999
[35,] 2.084688e-08 4.169376e-08 1.0000000
[36,] 1.419981e-08 2.839962e-08 1.0000000
[37,] 9.218279e-09 1.843656e-08 1.0000000
[38,] 4.817276e-09 9.634552e-09 1.0000000
[39,] 2.010211e-09 4.020421e-09 1.0000000
[40,] 6.487089e-10 1.297418e-09 1.0000000
[41,] 2.586214e-10 5.172428e-10 1.0000000
[42,] 8.445094e-11 1.689019e-10 1.0000000
[43,] 2.535459e-11 5.070918e-11 1.0000000
[44,] 7.893972e-12 1.578794e-11 1.0000000
[45,] 4.081007e-11 8.162015e-11 1.0000000
[46,] 2.119871e-11 4.239743e-11 1.0000000
[47,] 6.864604e-12 1.372921e-11 1.0000000
[48,] 3.872048e-12 7.744097e-12 1.0000000
[49,] 1.250222e-12 2.500444e-12 1.0000000
[50,] 3.905898e-13 7.811796e-13 1.0000000
[51,] 3.097954e-11 6.195908e-11 1.0000000
[52,] 1.082432e-11 2.164864e-11 1.0000000
[53,] 7.671970e-12 1.534394e-11 1.0000000
> postscript(file="/var/www/html/freestat/rcomp/tmp/1s2jq1292770153.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/2s2jq1292770153.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/3lb0b1292770153.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/4lb0b1292770153.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/5lb0b1292770153.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 = 72
Frequency = 1
1 2 3 4 5 6
31.212279 44.340152 145.376958 59.310279 31.715027 81.168023
7 8 9 10 11 12
-60.292804 -13.763123 43.061668 -2.774858 75.025004 -63.559549
13 14 15 16 17 18
82.539602 -43.865157 -46.283243 -29.468584 69.714542 94.332003
19 20 21 22 23 24
2.438892 -5.792966 4.478605 -48.613337 -118.318216 41.617330
25 26 27 28 29 30
59.532061 55.077682 -63.562975 -28.532999 20.486706 -98.646854
31 32 33 34 35 36
-47.041224 -60.854367 -113.565013 4.600189 18.964555 -22.800530
37 38 39 40 41 42
-80.882823 -40.534725 39.254912 59.519528 -88.142602 13.529952
43 44 45 46 47 48
-51.554679 -65.983643 34.232014 76.209706 -75.170863 -58.242836
49 50 51 52 53 54
22.657843 -61.479708 -43.022293 23.299532 -44.473025 154.395905
55 56 57 58 59 60
95.931452 -52.122038 101.607883 -48.555584 -0.735864 207.248430
61 62 63 64 65 66
86.902696 145.231961 -720.565922 74.920859 93.123826 13.915106
67 68 69 70 71 72
-12.027064 -21.832870 64.416961 -9.553502 134.267703 -63.041985
> postscript(file="/var/www/html/freestat/rcomp/tmp/6ekiw1292770153.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 31.212279 NA
1 44.340152 31.212279
2 145.376958 44.340152
3 59.310279 145.376958
4 31.715027 59.310279
5 81.168023 31.715027
6 -60.292804 81.168023
7 -13.763123 -60.292804
8 43.061668 -13.763123
9 -2.774858 43.061668
10 75.025004 -2.774858
11 -63.559549 75.025004
12 82.539602 -63.559549
13 -43.865157 82.539602
14 -46.283243 -43.865157
15 -29.468584 -46.283243
16 69.714542 -29.468584
17 94.332003 69.714542
18 2.438892 94.332003
19 -5.792966 2.438892
20 4.478605 -5.792966
21 -48.613337 4.478605
22 -118.318216 -48.613337
23 41.617330 -118.318216
24 59.532061 41.617330
25 55.077682 59.532061
26 -63.562975 55.077682
27 -28.532999 -63.562975
28 20.486706 -28.532999
29 -98.646854 20.486706
30 -47.041224 -98.646854
31 -60.854367 -47.041224
32 -113.565013 -60.854367
33 4.600189 -113.565013
34 18.964555 4.600189
35 -22.800530 18.964555
36 -80.882823 -22.800530
37 -40.534725 -80.882823
38 39.254912 -40.534725
39 59.519528 39.254912
40 -88.142602 59.519528
41 13.529952 -88.142602
42 -51.554679 13.529952
43 -65.983643 -51.554679
44 34.232014 -65.983643
45 76.209706 34.232014
46 -75.170863 76.209706
47 -58.242836 -75.170863
48 22.657843 -58.242836
49 -61.479708 22.657843
50 -43.022293 -61.479708
51 23.299532 -43.022293
52 -44.473025 23.299532
53 154.395905 -44.473025
54 95.931452 154.395905
55 -52.122038 95.931452
56 101.607883 -52.122038
57 -48.555584 101.607883
58 -0.735864 -48.555584
59 207.248430 -0.735864
60 86.902696 207.248430
61 145.231961 86.902696
62 -720.565922 145.231961
63 74.920859 -720.565922
64 93.123826 74.920859
65 13.915106 93.123826
66 -12.027064 13.915106
67 -21.832870 -12.027064
68 64.416961 -21.832870
69 -9.553502 64.416961
70 134.267703 -9.553502
71 -63.041985 134.267703
72 NA -63.041985
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 44.340152 31.212279
[2,] 145.376958 44.340152
[3,] 59.310279 145.376958
[4,] 31.715027 59.310279
[5,] 81.168023 31.715027
[6,] -60.292804 81.168023
[7,] -13.763123 -60.292804
[8,] 43.061668 -13.763123
[9,] -2.774858 43.061668
[10,] 75.025004 -2.774858
[11,] -63.559549 75.025004
[12,] 82.539602 -63.559549
[13,] -43.865157 82.539602
[14,] -46.283243 -43.865157
[15,] -29.468584 -46.283243
[16,] 69.714542 -29.468584
[17,] 94.332003 69.714542
[18,] 2.438892 94.332003
[19,] -5.792966 2.438892
[20,] 4.478605 -5.792966
[21,] -48.613337 4.478605
[22,] -118.318216 -48.613337
[23,] 41.617330 -118.318216
[24,] 59.532061 41.617330
[25,] 55.077682 59.532061
[26,] -63.562975 55.077682
[27,] -28.532999 -63.562975
[28,] 20.486706 -28.532999
[29,] -98.646854 20.486706
[30,] -47.041224 -98.646854
[31,] -60.854367 -47.041224
[32,] -113.565013 -60.854367
[33,] 4.600189 -113.565013
[34,] 18.964555 4.600189
[35,] -22.800530 18.964555
[36,] -80.882823 -22.800530
[37,] -40.534725 -80.882823
[38,] 39.254912 -40.534725
[39,] 59.519528 39.254912
[40,] -88.142602 59.519528
[41,] 13.529952 -88.142602
[42,] -51.554679 13.529952
[43,] -65.983643 -51.554679
[44,] 34.232014 -65.983643
[45,] 76.209706 34.232014
[46,] -75.170863 76.209706
[47,] -58.242836 -75.170863
[48,] 22.657843 -58.242836
[49,] -61.479708 22.657843
[50,] -43.022293 -61.479708
[51,] 23.299532 -43.022293
[52,] -44.473025 23.299532
[53,] 154.395905 -44.473025
[54,] 95.931452 154.395905
[55,] -52.122038 95.931452
[56,] 101.607883 -52.122038
[57,] -48.555584 101.607883
[58,] -0.735864 -48.555584
[59,] 207.248430 -0.735864
[60,] 86.902696 207.248430
[61,] 145.231961 86.902696
[62,] -720.565922 145.231961
[63,] 74.920859 -720.565922
[64,] 93.123826 74.920859
[65,] 13.915106 93.123826
[66,] -12.027064 13.915106
[67,] -21.832870 -12.027064
[68,] 64.416961 -21.832870
[69,] -9.553502 64.416961
[70,] 134.267703 -9.553502
[71,] -63.041985 134.267703
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 44.340152 31.212279
2 145.376958 44.340152
3 59.310279 145.376958
4 31.715027 59.310279
5 81.168023 31.715027
6 -60.292804 81.168023
7 -13.763123 -60.292804
8 43.061668 -13.763123
9 -2.774858 43.061668
10 75.025004 -2.774858
11 -63.559549 75.025004
12 82.539602 -63.559549
13 -43.865157 82.539602
14 -46.283243 -43.865157
15 -29.468584 -46.283243
16 69.714542 -29.468584
17 94.332003 69.714542
18 2.438892 94.332003
19 -5.792966 2.438892
20 4.478605 -5.792966
21 -48.613337 4.478605
22 -118.318216 -48.613337
23 41.617330 -118.318216
24 59.532061 41.617330
25 55.077682 59.532061
26 -63.562975 55.077682
27 -28.532999 -63.562975
28 20.486706 -28.532999
29 -98.646854 20.486706
30 -47.041224 -98.646854
31 -60.854367 -47.041224
32 -113.565013 -60.854367
33 4.600189 -113.565013
34 18.964555 4.600189
35 -22.800530 18.964555
36 -80.882823 -22.800530
37 -40.534725 -80.882823
38 39.254912 -40.534725
39 59.519528 39.254912
40 -88.142602 59.519528
41 13.529952 -88.142602
42 -51.554679 13.529952
43 -65.983643 -51.554679
44 34.232014 -65.983643
45 76.209706 34.232014
46 -75.170863 76.209706
47 -58.242836 -75.170863
48 22.657843 -58.242836
49 -61.479708 22.657843
50 -43.022293 -61.479708
51 23.299532 -43.022293
52 -44.473025 23.299532
53 154.395905 -44.473025
54 95.931452 154.395905
55 -52.122038 95.931452
56 101.607883 -52.122038
57 -48.555584 101.607883
58 -0.735864 -48.555584
59 207.248430 -0.735864
60 86.902696 207.248430
61 145.231961 86.902696
62 -720.565922 145.231961
63 74.920859 -720.565922
64 93.123826 74.920859
65 13.915106 93.123826
66 -12.027064 13.915106
67 -21.832870 -12.027064
68 64.416961 -21.832870
69 -9.553502 64.416961
70 134.267703 -9.553502
71 -63.041985 134.267703
> 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/76uzh1292770153.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/86uzh1292770153.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/96uzh1292770153.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/10h3gk1292770153.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/11lmxq1292770153.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/1264dw1292770153.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/13kebm1292770153.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/145wrs1292770153.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/15rfqg1292770153.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/16uf641292770153.tab")
+ }
>
> try(system("convert tmp/1s2jq1292770153.ps tmp/1s2jq1292770153.png",intern=TRUE))
character(0)
> try(system("convert tmp/2s2jq1292770153.ps tmp/2s2jq1292770153.png",intern=TRUE))
character(0)
> try(system("convert tmp/3lb0b1292770153.ps tmp/3lb0b1292770153.png",intern=TRUE))
character(0)
> try(system("convert tmp/4lb0b1292770153.ps tmp/4lb0b1292770153.png",intern=TRUE))
character(0)
> try(system("convert tmp/5lb0b1292770153.ps tmp/5lb0b1292770153.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ekiw1292770153.ps tmp/6ekiw1292770153.png",intern=TRUE))
character(0)
> try(system("convert tmp/76uzh1292770153.ps tmp/76uzh1292770153.png",intern=TRUE))
character(0)
> try(system("convert tmp/86uzh1292770153.ps tmp/86uzh1292770153.png",intern=TRUE))
character(0)
> try(system("convert tmp/96uzh1292770153.ps tmp/96uzh1292770153.png",intern=TRUE))
character(0)
> try(system("convert tmp/10h3gk1292770153.ps tmp/10h3gk1292770153.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.107 2.514 5.381