R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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.
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(1
+ ,216234.00
+ ,627
+ ,2
+ ,213586.00
+ ,696
+ ,3
+ ,209465.00
+ ,825
+ ,4
+ ,204045.00
+ ,677
+ ,5
+ ,200237.00
+ ,656
+ ,6
+ ,203666.00
+ ,785
+ ,7
+ ,241476.00
+ ,412
+ ,8
+ ,260307.00
+ ,352
+ ,9
+ ,243324.00
+ ,839
+ ,10
+ ,244460.00
+ ,729
+ ,11
+ ,233575.00
+ ,696
+ ,12
+ ,237217.00
+ ,641
+ ,1
+ ,235243.00
+ ,695
+ ,2
+ ,230354.00
+ ,638
+ ,3
+ ,227184.00
+ ,762
+ ,4
+ ,221678.00
+ ,635
+ ,5
+ ,217142.00
+ ,721
+ ,6
+ ,219452.00
+ ,854
+ ,7
+ ,256446.00
+ ,418
+ ,8
+ ,265845.00
+ ,367
+ ,9
+ ,248624.00
+ ,824
+ ,10
+ ,241114.00
+ ,687
+ ,11
+ ,229245.00
+ ,601
+ ,12
+ ,231805.00
+ ,676
+ ,1
+ ,219277.00
+ ,740
+ ,2
+ ,219313.00
+ ,691
+ ,3
+ ,212610.00
+ ,683
+ ,4
+ ,214771.00
+ ,594
+ ,5
+ ,211142.00
+ ,729
+ ,6
+ ,211457.00
+ ,731
+ ,7
+ ,240048.00
+ ,386
+ ,8
+ ,240636.00
+ ,331
+ ,9
+ ,230580.00
+ ,707
+ ,10
+ ,208795.00
+ ,715
+ ,11
+ ,197922.00
+ ,657
+ ,12
+ ,194596.00
+ ,653
+ ,1
+ ,194581.00
+ ,642
+ ,2
+ ,185686.00
+ ,643
+ ,3
+ ,178106.00
+ ,718
+ ,4
+ ,172608.00
+ ,654
+ ,5
+ ,167302.00
+ ,632
+ ,6
+ ,168053.00
+ ,731
+ ,7
+ ,202300.00
+ ,392
+ ,8
+ ,202388.00
+ ,344
+ ,9
+ ,182516.00
+ ,792
+ ,10
+ ,173476.00
+ ,852
+ ,11
+ ,166444.00
+ ,649
+ ,12
+ ,171297.00
+ ,629
+ ,1
+ ,169701.00
+ ,685
+ ,2
+ ,164182.00
+ ,617
+ ,3
+ ,161914.00
+ ,715
+ ,4
+ ,159612.00
+ ,715
+ ,5
+ ,151001.00
+ ,629
+ ,6
+ ,158114.00
+ ,916
+ ,7
+ ,186530.00
+ ,531
+ ,8
+ ,187069.00
+ ,357
+ ,9
+ ,174330.00
+ ,917
+ ,10
+ ,169362.00
+ ,828
+ ,11
+ ,166827.00
+ ,708
+ ,12
+ ,178037.00
+ ,858
+ ,1
+ ,186413.00
+ ,775
+ ,2
+ ,189226.00
+ ,785
+ ,3
+ ,191563.00
+ ,1006
+ ,4
+ ,188906.00
+ ,789
+ ,5
+ ,186005.00
+ ,734
+ ,6
+ ,195309.00
+ ,906
+ ,7
+ ,223532.00
+ ,532
+ ,8
+ ,226899.00
+ ,387
+ ,9
+ ,214126.00
+ ,991
+ ,10
+ ,206903.00
+ ,841
+ ,11
+ ,204442.00
+ ,892
+ ,12
+ ,220375.00
+ ,782)
+ ,dim=c(3
+ ,72)
+ ,dimnames=list(c('month'
+ ,'werklozen'
+ ,'faillissementen')
+ ,1:72))
> y <- array(NA,dim=c(3,72),dimnames=list(c('month','werklozen','faillissementen'),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 = '2'
> #'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
werklozen month faillissementen
1 216234 1 627
2 213586 2 696
3 209465 3 825
4 204045 4 677
5 200237 5 656
6 203666 6 785
7 241476 7 412
8 260307 8 352
9 243324 9 839
10 244460 10 729
11 233575 11 696
12 237217 12 641
13 235243 1 695
14 230354 2 638
15 227184 3 762
16 221678 4 635
17 217142 5 721
18 219452 6 854
19 256446 7 418
20 265845 8 367
21 248624 9 824
22 241114 10 687
23 229245 11 601
24 231805 12 676
25 219277 1 740
26 219313 2 691
27 212610 3 683
28 214771 4 594
29 211142 5 729
30 211457 6 731
31 240048 7 386
32 240636 8 331
33 230580 9 707
34 208795 10 715
35 197922 11 657
36 194596 12 653
37 194581 1 642
38 185686 2 643
39 178106 3 718
40 172608 4 654
41 167302 5 632
42 168053 6 731
43 202300 7 392
44 202388 8 344
45 182516 9 792
46 173476 10 852
47 166444 11 649
48 171297 12 629
49 169701 1 685
50 164182 2 617
51 161914 3 715
52 159612 4 715
53 151001 5 629
54 158114 6 916
55 186530 7 531
56 187069 8 357
57 174330 9 917
58 169362 10 828
59 166827 11 708
60 178037 12 858
61 186413 1 775
62 189226 2 785
63 191563 3 1006
64 188906 4 789
65 186005 5 734
66 195309 6 906
67 223532 7 532
68 226899 8 387
69 214126 9 991
70 206903 10 841
71 204442 11 892
72 220375 12 782
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) month faillissementen
237704.10 1197.29 -59.53
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-55248 -21078 6086 19930 49194
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 237704.10 15295.65 15.541 < 2e-16 ***
month 1197.29 903.78 1.325 0.18962
faillissementen -59.53 20.07 -2.966 0.00415 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 26470 on 69 degrees of freedom
Multiple R-squared: 0.1343, Adjusted R-squared: 0.1092
F-statistic: 5.353 on 2 and 69 DF, p-value: 0.006902
> 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.001218449 0.0024368979 0.9987815510
[2,] 0.030531715 0.0610634295 0.9694682853
[3,] 0.024928934 0.0498578684 0.9750710658
[4,] 0.066367223 0.1327344464 0.9336327768
[5,] 0.036323422 0.0726468449 0.9636765776
[6,] 0.023244959 0.0464899172 0.9767550414
[7,] 0.014753760 0.0295075207 0.9852462397
[8,] 0.032733685 0.0654673704 0.9672663148
[9,] 0.024128138 0.0482562761 0.9758718620
[10,] 0.018624215 0.0372484299 0.9813757850
[11,] 0.010549758 0.0210995155 0.9894502423
[12,] 0.005979433 0.0119588657 0.9940205671
[13,] 0.003474477 0.0069489548 0.9965255226
[14,] 0.003117213 0.0062344260 0.9968827870
[15,] 0.003604169 0.0072083372 0.9963958314
[16,] 0.008843285 0.0176865705 0.9911567148
[17,] 0.007456771 0.0149135419 0.9925432291
[18,] 0.008395827 0.0167916546 0.9916041727
[19,] 0.007820026 0.0156400519 0.9921799740
[20,] 0.006568648 0.0131372954 0.9934313523
[21,] 0.005307589 0.0106151782 0.9946924109
[22,] 0.004398997 0.0087979930 0.9956010035
[23,] 0.004310004 0.0086200075 0.9956899962
[24,] 0.003917409 0.0078348176 0.9960825912
[25,] 0.003776342 0.0075526847 0.9962236576
[26,] 0.004716487 0.0094329733 0.9952835133
[27,] 0.007718094 0.0154361881 0.9922819060
[28,] 0.011758644 0.0235172882 0.9882413559
[29,] 0.020146875 0.0402937492 0.9798531254
[30,] 0.059214221 0.1184284417 0.9407857792
[31,] 0.123079962 0.2461599233 0.8769200383
[32,] 0.145899529 0.2917990571 0.8541004714
[33,] 0.189522542 0.3790450841 0.8104774579
[34,] 0.240497441 0.4809948812 0.7595025594
[35,] 0.352019569 0.7040391380 0.6479804310
[36,] 0.517428650 0.9651426998 0.4825713499
[37,] 0.611015979 0.7779680426 0.3889840213
[38,] 0.619866225 0.7602675503 0.3801337752
[39,] 0.626144181 0.7477116381 0.3738558191
[40,] 0.616833552 0.7663328969 0.3831664485
[41,] 0.637093373 0.7258132545 0.3629066273
[42,] 0.746679775 0.5066404500 0.2533202250
[43,] 0.812069224 0.3758615512 0.1879307756
[44,] 0.794431479 0.4111370427 0.2055685213
[45,] 0.802912987 0.3941740267 0.1970870134
[46,] 0.812437831 0.3751243389 0.1875621695
[47,] 0.843013188 0.3139736233 0.1569868116
[48,] 0.936701285 0.1265974309 0.0632987154
[49,] 0.955473557 0.0890528857 0.0445264429
[50,] 0.943105600 0.1137888002 0.0568944001
[51,] 0.948333449 0.1033331022 0.0516665511
[52,] 0.935630118 0.1287397633 0.0643698817
[53,] 0.952555770 0.0948884600 0.0474442300
[54,] 0.994104501 0.0117909979 0.0058954989
[55,] 0.999765536 0.0004689289 0.0002344645
[56,] 0.999166900 0.0016661992 0.0008330996
[57,] 0.997211860 0.0055762797 0.0027881399
[58,] 0.995385863 0.0092282739 0.0046141370
[59,] 0.985318299 0.0293634024 0.0146817012
[60,] 0.987848428 0.0243031441 0.0121515721
[61,] 0.977457703 0.0450845934 0.0225422967
> postscript(file="/var/www/html/rcomp/tmp/1dxp51292612205.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/rcomp/tmp/2dxp51292612205.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/rcomp/tmp/3op681292612205.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/rcomp/tmp/4op681292612205.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/rcomp/tmp/5op681292612205.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 7
14655.546 14917.561 17278.149 1850.974 -4404.371 5506.218 19915.643
8 9 10 11 12 13 14
33977.776 44786.751 38177.572 24130.913 23301.676 37712.329 28233.040
15 16 17 18 19 20 21
31246.997 16983.873 16369.834 25399.527 35242.801 40408.669 49193.857
22 23 24 25 26 27 28
32331.471 14145.922 19973.094 24425.009 20346.930 11970.426 7636.298
29 30 31 32 33 34 35
10846.043 10082.802 16939.962 13056.725 24185.290 1679.205 -13843.610
36 37 38 39 40 41 42
-18605.009 -6104.560 -16137.329 -20450.156 -30955.129 -38768.000 -33321.198
43 44 45 46 47 48 49
-20450.881 -24417.434 -18818.982 -25484.703 -45797.820 -43332.638 -28424.933
50 51 52 53 54 55 56
-39189.010 -36820.735 -40320.029 -55247.578 -32247.848 -27946.737 -38962.593
57 58 59 60 61 62 63
-19564.204 -31027.332 -41902.773 -22961.134 -6355.573 -4144.605 10150.395
64 65 66 67 68 69 70
-6621.089 -13993.325 4351.890 9114.790 2653.194 24636.736 7287.509
71 72
6665.052 14852.873
> postscript(file="/var/www/html/rcomp/tmp/6hynt1292612205.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 14655.546 NA
1 14917.561 14655.546
2 17278.149 14917.561
3 1850.974 17278.149
4 -4404.371 1850.974
5 5506.218 -4404.371
6 19915.643 5506.218
7 33977.776 19915.643
8 44786.751 33977.776
9 38177.572 44786.751
10 24130.913 38177.572
11 23301.676 24130.913
12 37712.329 23301.676
13 28233.040 37712.329
14 31246.997 28233.040
15 16983.873 31246.997
16 16369.834 16983.873
17 25399.527 16369.834
18 35242.801 25399.527
19 40408.669 35242.801
20 49193.857 40408.669
21 32331.471 49193.857
22 14145.922 32331.471
23 19973.094 14145.922
24 24425.009 19973.094
25 20346.930 24425.009
26 11970.426 20346.930
27 7636.298 11970.426
28 10846.043 7636.298
29 10082.802 10846.043
30 16939.962 10082.802
31 13056.725 16939.962
32 24185.290 13056.725
33 1679.205 24185.290
34 -13843.610 1679.205
35 -18605.009 -13843.610
36 -6104.560 -18605.009
37 -16137.329 -6104.560
38 -20450.156 -16137.329
39 -30955.129 -20450.156
40 -38768.000 -30955.129
41 -33321.198 -38768.000
42 -20450.881 -33321.198
43 -24417.434 -20450.881
44 -18818.982 -24417.434
45 -25484.703 -18818.982
46 -45797.820 -25484.703
47 -43332.638 -45797.820
48 -28424.933 -43332.638
49 -39189.010 -28424.933
50 -36820.735 -39189.010
51 -40320.029 -36820.735
52 -55247.578 -40320.029
53 -32247.848 -55247.578
54 -27946.737 -32247.848
55 -38962.593 -27946.737
56 -19564.204 -38962.593
57 -31027.332 -19564.204
58 -41902.773 -31027.332
59 -22961.134 -41902.773
60 -6355.573 -22961.134
61 -4144.605 -6355.573
62 10150.395 -4144.605
63 -6621.089 10150.395
64 -13993.325 -6621.089
65 4351.890 -13993.325
66 9114.790 4351.890
67 2653.194 9114.790
68 24636.736 2653.194
69 7287.509 24636.736
70 6665.052 7287.509
71 14852.873 6665.052
72 NA 14852.873
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 14917.561 14655.546
[2,] 17278.149 14917.561
[3,] 1850.974 17278.149
[4,] -4404.371 1850.974
[5,] 5506.218 -4404.371
[6,] 19915.643 5506.218
[7,] 33977.776 19915.643
[8,] 44786.751 33977.776
[9,] 38177.572 44786.751
[10,] 24130.913 38177.572
[11,] 23301.676 24130.913
[12,] 37712.329 23301.676
[13,] 28233.040 37712.329
[14,] 31246.997 28233.040
[15,] 16983.873 31246.997
[16,] 16369.834 16983.873
[17,] 25399.527 16369.834
[18,] 35242.801 25399.527
[19,] 40408.669 35242.801
[20,] 49193.857 40408.669
[21,] 32331.471 49193.857
[22,] 14145.922 32331.471
[23,] 19973.094 14145.922
[24,] 24425.009 19973.094
[25,] 20346.930 24425.009
[26,] 11970.426 20346.930
[27,] 7636.298 11970.426
[28,] 10846.043 7636.298
[29,] 10082.802 10846.043
[30,] 16939.962 10082.802
[31,] 13056.725 16939.962
[32,] 24185.290 13056.725
[33,] 1679.205 24185.290
[34,] -13843.610 1679.205
[35,] -18605.009 -13843.610
[36,] -6104.560 -18605.009
[37,] -16137.329 -6104.560
[38,] -20450.156 -16137.329
[39,] -30955.129 -20450.156
[40,] -38768.000 -30955.129
[41,] -33321.198 -38768.000
[42,] -20450.881 -33321.198
[43,] -24417.434 -20450.881
[44,] -18818.982 -24417.434
[45,] -25484.703 -18818.982
[46,] -45797.820 -25484.703
[47,] -43332.638 -45797.820
[48,] -28424.933 -43332.638
[49,] -39189.010 -28424.933
[50,] -36820.735 -39189.010
[51,] -40320.029 -36820.735
[52,] -55247.578 -40320.029
[53,] -32247.848 -55247.578
[54,] -27946.737 -32247.848
[55,] -38962.593 -27946.737
[56,] -19564.204 -38962.593
[57,] -31027.332 -19564.204
[58,] -41902.773 -31027.332
[59,] -22961.134 -41902.773
[60,] -6355.573 -22961.134
[61,] -4144.605 -6355.573
[62,] 10150.395 -4144.605
[63,] -6621.089 10150.395
[64,] -13993.325 -6621.089
[65,] 4351.890 -13993.325
[66,] 9114.790 4351.890
[67,] 2653.194 9114.790
[68,] 24636.736 2653.194
[69,] 7287.509 24636.736
[70,] 6665.052 7287.509
[71,] 14852.873 6665.052
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 14917.561 14655.546
2 17278.149 14917.561
3 1850.974 17278.149
4 -4404.371 1850.974
5 5506.218 -4404.371
6 19915.643 5506.218
7 33977.776 19915.643
8 44786.751 33977.776
9 38177.572 44786.751
10 24130.913 38177.572
11 23301.676 24130.913
12 37712.329 23301.676
13 28233.040 37712.329
14 31246.997 28233.040
15 16983.873 31246.997
16 16369.834 16983.873
17 25399.527 16369.834
18 35242.801 25399.527
19 40408.669 35242.801
20 49193.857 40408.669
21 32331.471 49193.857
22 14145.922 32331.471
23 19973.094 14145.922
24 24425.009 19973.094
25 20346.930 24425.009
26 11970.426 20346.930
27 7636.298 11970.426
28 10846.043 7636.298
29 10082.802 10846.043
30 16939.962 10082.802
31 13056.725 16939.962
32 24185.290 13056.725
33 1679.205 24185.290
34 -13843.610 1679.205
35 -18605.009 -13843.610
36 -6104.560 -18605.009
37 -16137.329 -6104.560
38 -20450.156 -16137.329
39 -30955.129 -20450.156
40 -38768.000 -30955.129
41 -33321.198 -38768.000
42 -20450.881 -33321.198
43 -24417.434 -20450.881
44 -18818.982 -24417.434
45 -25484.703 -18818.982
46 -45797.820 -25484.703
47 -43332.638 -45797.820
48 -28424.933 -43332.638
49 -39189.010 -28424.933
50 -36820.735 -39189.010
51 -40320.029 -36820.735
52 -55247.578 -40320.029
53 -32247.848 -55247.578
54 -27946.737 -32247.848
55 -38962.593 -27946.737
56 -19564.204 -38962.593
57 -31027.332 -19564.204
58 -41902.773 -31027.332
59 -22961.134 -41902.773
60 -6355.573 -22961.134
61 -4144.605 -6355.573
62 10150.395 -4144.605
63 -6621.089 10150.395
64 -13993.325 -6621.089
65 4351.890 -13993.325
66 9114.790 4351.890
67 2653.194 9114.790
68 24636.736 2653.194
69 7287.509 24636.736
70 6665.052 7287.509
71 14852.873 6665.052
> 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/rcomp/tmp/7r75w1292612205.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/rcomp/tmp/8r75w1292612205.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/rcomp/tmp/9r75w1292612205.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/rcomp/tmp/102gmh1292612205.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/rcomp/tmp/115z2m1292612205.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/rcomp/tmp/129h1a1292612205.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/rcomp/tmp/1359y11292612205.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/rcomp/tmp/148axp1292612205.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/rcomp/tmp/15uswd1292612205.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/rcomp/tmp/16xbuj1292612205.tab")
+ }
> try(system("convert tmp/1dxp51292612205.ps tmp/1dxp51292612205.png",intern=TRUE))
character(0)
> try(system("convert tmp/2dxp51292612205.ps tmp/2dxp51292612205.png",intern=TRUE))
character(0)
> try(system("convert tmp/3op681292612205.ps tmp/3op681292612205.png",intern=TRUE))
character(0)
> try(system("convert tmp/4op681292612205.ps tmp/4op681292612205.png",intern=TRUE))
character(0)
> try(system("convert tmp/5op681292612205.ps tmp/5op681292612205.png",intern=TRUE))
character(0)
> try(system("convert tmp/6hynt1292612205.ps tmp/6hynt1292612205.png",intern=TRUE))
character(0)
> try(system("convert tmp/7r75w1292612205.ps tmp/7r75w1292612205.png",intern=TRUE))
character(0)
> try(system("convert tmp/8r75w1292612205.ps tmp/8r75w1292612205.png",intern=TRUE))
character(0)
> try(system("convert tmp/9r75w1292612205.ps tmp/9r75w1292612205.png",intern=TRUE))
character(0)
> try(system("convert tmp/102gmh1292612205.ps tmp/102gmh1292612205.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.645 1.650 6.865