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(98.60,627,98.97,696,99.11,825,99.64,677,100.03,656,99.98,785,100.32,412,100.44,352,100.51,839,101.00,729,100.88,696,100.55,641,100.83,695,101.51,638,102.16,762,102.39,635,102.54,721,102.85,854,103.47,418,103.57,367,103.69,824,103.50,687,103.47,601,103.45,676,103.48,740,103.93,691,103.89,683,104.40,594,104.79,729,104.77,731,105.13,386,105.26,331,104.96,707,104.75,715,105.01,657,105.15,653,105.20,642,105.77,643,105.78,718,106.26,654,106.13,632,106.12,731,106.57,392,106.44,344,106.54,792,107.10,852,108.10,649,108.40,629,108.84,685,109.62,617,110.42,715,110.67,715,111.66,629,112.28,916,112.87,531,112.18,357,112.36,917,112.16,828,111.49,708,111.25,858,111.36,775,111.74,785,111.10,1006,111.33,789,111.25,734,111.04,906,110.97,532,111.31,387,111.02,991,111.07,841,111.36,892,111.54,782),dim=c(2,72),dimnames=list(c('CPI','Faillissementen'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('CPI','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 = '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
CPI Faillissementen
1 98.60 627
2 98.97 696
3 99.11 825
4 99.64 677
5 100.03 656
6 99.98 785
7 100.32 412
8 100.44 352
9 100.51 839
10 101.00 729
11 100.88 696
12 100.55 641
13 100.83 695
14 101.51 638
15 102.16 762
16 102.39 635
17 102.54 721
18 102.85 854
19 103.47 418
20 103.57 367
21 103.69 824
22 103.50 687
23 103.47 601
24 103.45 676
25 103.48 740
26 103.93 691
27 103.89 683
28 104.40 594
29 104.79 729
30 104.77 731
31 105.13 386
32 105.26 331
33 104.96 707
34 104.75 715
35 105.01 657
36 105.15 653
37 105.20 642
38 105.77 643
39 105.78 718
40 106.26 654
41 106.13 632
42 106.12 731
43 106.57 392
44 106.44 344
45 106.54 792
46 107.10 852
47 108.10 649
48 108.40 629
49 108.84 685
50 109.62 617
51 110.42 715
52 110.67 715
53 111.66 629
54 112.28 916
55 112.87 531
56 112.18 357
57 112.36 917
58 112.16 828
59 111.49 708
60 111.25 858
61 111.36 775
62 111.74 785
63 111.10 1006
64 111.33 789
65 111.25 734
66 111.04 906
67 110.97 532
68 111.31 387
69 111.02 991
70 111.07 841
71 111.36 892
72 111.54 782
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Faillissementen
1.021e+02 5.844e-03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-7.8210 -3.2840 -0.3937 4.0395 7.9838
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.021e+02 2.227e+00 45.859 <2e-16 ***
Faillissementen 5.844e-03 3.189e-03 1.832 0.0712 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.207 on 70 degrees of freedom
Multiple R-squared: 0.04576, Adjusted R-squared: 0.03213
F-statistic: 3.357 on 1 and 70 DF, p-value: 0.07118
> 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,] 6.430621e-03 1.286124e-02 9.935694e-01
[2,] 1.599632e-03 3.199264e-03 9.984004e-01
[3,] 4.249176e-04 8.498352e-04 9.995751e-01
[4,] 7.193266e-05 1.438653e-04 9.999281e-01
[5,] 6.768015e-05 1.353603e-04 9.999323e-01
[6,] 6.125174e-05 1.225035e-04 9.999387e-01
[7,] 3.275713e-05 6.551427e-05 9.999672e-01
[8,] 1.143919e-05 2.287838e-05 9.999886e-01
[9,] 5.607627e-06 1.121525e-05 9.999944e-01
[10,] 6.087969e-06 1.217594e-05 9.999939e-01
[11,] 1.711101e-05 3.422202e-05 9.999829e-01
[12,] 3.505982e-05 7.011964e-05 9.999649e-01
[13,] 6.579898e-05 1.315980e-04 9.999342e-01
[14,] 1.271025e-04 2.542049e-04 9.998729e-01
[15,] 3.356305e-04 6.712611e-04 9.996644e-01
[16,] 3.989066e-04 7.978133e-04 9.996011e-01
[17,] 1.099879e-03 2.199757e-03 9.989001e-01
[18,] 1.628292e-03 3.256583e-03 9.983717e-01
[19,] 1.976384e-03 3.952767e-03 9.980236e-01
[20,] 2.616949e-03 5.233897e-03 9.973831e-01
[21,] 3.864606e-03 7.729213e-03 9.961354e-01
[22,] 5.948453e-03 1.189691e-02 9.940515e-01
[23,] 8.809429e-03 1.761886e-02 9.911906e-01
[24,] 1.279956e-02 2.559913e-02 9.872004e-01
[25,] 2.314719e-02 4.629437e-02 9.768528e-01
[26,] 3.872540e-02 7.745080e-02 9.612746e-01
[27,] 4.536979e-02 9.073959e-02 9.546302e-01
[28,] 4.568805e-02 9.137611e-02 9.543119e-01
[29,] 7.122505e-02 1.424501e-01 9.287750e-01
[30,] 1.094339e-01 2.188678e-01 8.905661e-01
[31,] 1.553226e-01 3.106452e-01 8.446774e-01
[32,] 2.183832e-01 4.367664e-01 7.816168e-01
[33,] 3.016678e-01 6.033357e-01 6.983322e-01
[34,] 4.023746e-01 8.047492e-01 5.976254e-01
[35,] 5.452852e-01 9.094296e-01 4.547148e-01
[36,] 6.652421e-01 6.695159e-01 3.347579e-01
[37,] 7.778675e-01 4.442649e-01 2.221325e-01
[38,] 8.949713e-01 2.100575e-01 1.050287e-01
[39,] 9.321026e-01 1.357948e-01 6.789738e-02
[40,] 9.713242e-01 5.735163e-02 2.867581e-02
[41,] 9.963345e-01 7.330964e-03 3.665482e-03
[42,] 9.998115e-01 3.770893e-04 1.885446e-04
[43,] 9.999851e-01 2.977891e-05 1.488945e-05
[44,] 9.999994e-01 1.144578e-06 5.722888e-07
[45,] 1.000000e+00 2.022733e-08 1.011366e-08
[46,] 1.000000e+00 7.784208e-10 3.892104e-10
[47,] 1.000000e+00 2.366548e-10 1.183274e-10
[48,] 1.000000e+00 1.376624e-10 6.883121e-11
[49,] 1.000000e+00 3.057455e-10 1.528727e-10
[50,] 1.000000e+00 2.088430e-10 1.044215e-10
[51,] 1.000000e+00 2.533021e-11 1.266510e-11
[52,] 1.000000e+00 3.296943e-11 1.648472e-11
[53,] 1.000000e+00 5.334865e-12 2.667432e-12
[54,] 1.000000e+00 3.847896e-13 1.923948e-13
[55,] 1.000000e+00 3.392435e-12 1.696218e-12
[56,] 1.000000e+00 4.879611e-11 2.439806e-11
[57,] 1.000000e+00 6.370827e-10 3.185413e-10
[58,] 1.000000e+00 7.165132e-10 3.582566e-10
[59,] 1.000000e+00 1.356740e-08 6.783701e-09
[60,] 9.999999e-01 2.139005e-07 1.069503e-07
[61,] 9.999981e-01 3.839594e-06 1.919797e-06
[62,] 9.999708e-01 5.834027e-05 2.917014e-05
[63,] 9.997547e-01 4.906748e-04 2.453374e-04
> postscript(file="/var/www/html/rcomp/tmp/1r0pq1291122100.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/2r0pq1291122100.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/3r0pq1291122100.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/42r6t1291122100.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/52r6t1291122100.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.17393974 -7.20714890 -7.82097473 -6.42612029 -5.91340446 -6.71723029
7 8 9 10 11 12
-4.19756335 -3.72694669 -6.50278529 -5.36998807 -5.29714890 -5.30575029
13 14 15 16 17 18
-5.34130529 -4.32821946 -4.40282724 -3.43068863 -3.78323918 -4.25043946
19 20 21 22 23 24
-1.08262502 -0.68460086 -3.23513112 -2.62455640 -2.15200585 -2.61027668
25 26 27 28 29 30
-2.95426779 -2.21793085 -2.21118196 -1.18110057 -1.57998807 -1.61167529
31 32 33 34 35 36
0.76437053 1.21576914 -1.28142863 -1.53817751 -0.93924807 -0.77587363
37 38 39 40 41 42
-0.66159391 -0.09743752 -0.52570835 0.32828276 0.32684221 -0.26167529
43 44 45 46 47 48
2.16930887 2.31980220 -0.19813557 0.01124777 2.19750082 2.61437304
49 50 51 52 53 54
2.72713082 3.90449637 4.13182249 4.38182249 5.87437304 4.81725666
55 56 57 58 59 60
7.65704693 7.98383526 4.89141305 5.21149443 5.24272776 4.12618610
61 62 63 64 65 66
4.72120582 5.04276971 3.11133166 4.60939526 4.85079387 3.63569277
67 68 69 70 71 72
5.75120331 6.93852692 3.11898582 4.04552749 4.03750332 4.86030054
> postscript(file="/var/www/html/rcomp/tmp/62r6t1291122100.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 -7.17393974 NA
1 -7.20714890 -7.17393974
2 -7.82097473 -7.20714890
3 -6.42612029 -7.82097473
4 -5.91340446 -6.42612029
5 -6.71723029 -5.91340446
6 -4.19756335 -6.71723029
7 -3.72694669 -4.19756335
8 -6.50278529 -3.72694669
9 -5.36998807 -6.50278529
10 -5.29714890 -5.36998807
11 -5.30575029 -5.29714890
12 -5.34130529 -5.30575029
13 -4.32821946 -5.34130529
14 -4.40282724 -4.32821946
15 -3.43068863 -4.40282724
16 -3.78323918 -3.43068863
17 -4.25043946 -3.78323918
18 -1.08262502 -4.25043946
19 -0.68460086 -1.08262502
20 -3.23513112 -0.68460086
21 -2.62455640 -3.23513112
22 -2.15200585 -2.62455640
23 -2.61027668 -2.15200585
24 -2.95426779 -2.61027668
25 -2.21793085 -2.95426779
26 -2.21118196 -2.21793085
27 -1.18110057 -2.21118196
28 -1.57998807 -1.18110057
29 -1.61167529 -1.57998807
30 0.76437053 -1.61167529
31 1.21576914 0.76437053
32 -1.28142863 1.21576914
33 -1.53817751 -1.28142863
34 -0.93924807 -1.53817751
35 -0.77587363 -0.93924807
36 -0.66159391 -0.77587363
37 -0.09743752 -0.66159391
38 -0.52570835 -0.09743752
39 0.32828276 -0.52570835
40 0.32684221 0.32828276
41 -0.26167529 0.32684221
42 2.16930887 -0.26167529
43 2.31980220 2.16930887
44 -0.19813557 2.31980220
45 0.01124777 -0.19813557
46 2.19750082 0.01124777
47 2.61437304 2.19750082
48 2.72713082 2.61437304
49 3.90449637 2.72713082
50 4.13182249 3.90449637
51 4.38182249 4.13182249
52 5.87437304 4.38182249
53 4.81725666 5.87437304
54 7.65704693 4.81725666
55 7.98383526 7.65704693
56 4.89141305 7.98383526
57 5.21149443 4.89141305
58 5.24272776 5.21149443
59 4.12618610 5.24272776
60 4.72120582 4.12618610
61 5.04276971 4.72120582
62 3.11133166 5.04276971
63 4.60939526 3.11133166
64 4.85079387 4.60939526
65 3.63569277 4.85079387
66 5.75120331 3.63569277
67 6.93852692 5.75120331
68 3.11898582 6.93852692
69 4.04552749 3.11898582
70 4.03750332 4.04552749
71 4.86030054 4.03750332
72 NA 4.86030054
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -7.20714890 -7.17393974
[2,] -7.82097473 -7.20714890
[3,] -6.42612029 -7.82097473
[4,] -5.91340446 -6.42612029
[5,] -6.71723029 -5.91340446
[6,] -4.19756335 -6.71723029
[7,] -3.72694669 -4.19756335
[8,] -6.50278529 -3.72694669
[9,] -5.36998807 -6.50278529
[10,] -5.29714890 -5.36998807
[11,] -5.30575029 -5.29714890
[12,] -5.34130529 -5.30575029
[13,] -4.32821946 -5.34130529
[14,] -4.40282724 -4.32821946
[15,] -3.43068863 -4.40282724
[16,] -3.78323918 -3.43068863
[17,] -4.25043946 -3.78323918
[18,] -1.08262502 -4.25043946
[19,] -0.68460086 -1.08262502
[20,] -3.23513112 -0.68460086
[21,] -2.62455640 -3.23513112
[22,] -2.15200585 -2.62455640
[23,] -2.61027668 -2.15200585
[24,] -2.95426779 -2.61027668
[25,] -2.21793085 -2.95426779
[26,] -2.21118196 -2.21793085
[27,] -1.18110057 -2.21118196
[28,] -1.57998807 -1.18110057
[29,] -1.61167529 -1.57998807
[30,] 0.76437053 -1.61167529
[31,] 1.21576914 0.76437053
[32,] -1.28142863 1.21576914
[33,] -1.53817751 -1.28142863
[34,] -0.93924807 -1.53817751
[35,] -0.77587363 -0.93924807
[36,] -0.66159391 -0.77587363
[37,] -0.09743752 -0.66159391
[38,] -0.52570835 -0.09743752
[39,] 0.32828276 -0.52570835
[40,] 0.32684221 0.32828276
[41,] -0.26167529 0.32684221
[42,] 2.16930887 -0.26167529
[43,] 2.31980220 2.16930887
[44,] -0.19813557 2.31980220
[45,] 0.01124777 -0.19813557
[46,] 2.19750082 0.01124777
[47,] 2.61437304 2.19750082
[48,] 2.72713082 2.61437304
[49,] 3.90449637 2.72713082
[50,] 4.13182249 3.90449637
[51,] 4.38182249 4.13182249
[52,] 5.87437304 4.38182249
[53,] 4.81725666 5.87437304
[54,] 7.65704693 4.81725666
[55,] 7.98383526 7.65704693
[56,] 4.89141305 7.98383526
[57,] 5.21149443 4.89141305
[58,] 5.24272776 5.21149443
[59,] 4.12618610 5.24272776
[60,] 4.72120582 4.12618610
[61,] 5.04276971 4.72120582
[62,] 3.11133166 5.04276971
[63,] 4.60939526 3.11133166
[64,] 4.85079387 4.60939526
[65,] 3.63569277 4.85079387
[66,] 5.75120331 3.63569277
[67,] 6.93852692 5.75120331
[68,] 3.11898582 6.93852692
[69,] 4.04552749 3.11898582
[70,] 4.03750332 4.04552749
[71,] 4.86030054 4.03750332
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -7.20714890 -7.17393974
2 -7.82097473 -7.20714890
3 -6.42612029 -7.82097473
4 -5.91340446 -6.42612029
5 -6.71723029 -5.91340446
6 -4.19756335 -6.71723029
7 -3.72694669 -4.19756335
8 -6.50278529 -3.72694669
9 -5.36998807 -6.50278529
10 -5.29714890 -5.36998807
11 -5.30575029 -5.29714890
12 -5.34130529 -5.30575029
13 -4.32821946 -5.34130529
14 -4.40282724 -4.32821946
15 -3.43068863 -4.40282724
16 -3.78323918 -3.43068863
17 -4.25043946 -3.78323918
18 -1.08262502 -4.25043946
19 -0.68460086 -1.08262502
20 -3.23513112 -0.68460086
21 -2.62455640 -3.23513112
22 -2.15200585 -2.62455640
23 -2.61027668 -2.15200585
24 -2.95426779 -2.61027668
25 -2.21793085 -2.95426779
26 -2.21118196 -2.21793085
27 -1.18110057 -2.21118196
28 -1.57998807 -1.18110057
29 -1.61167529 -1.57998807
30 0.76437053 -1.61167529
31 1.21576914 0.76437053
32 -1.28142863 1.21576914
33 -1.53817751 -1.28142863
34 -0.93924807 -1.53817751
35 -0.77587363 -0.93924807
36 -0.66159391 -0.77587363
37 -0.09743752 -0.66159391
38 -0.52570835 -0.09743752
39 0.32828276 -0.52570835
40 0.32684221 0.32828276
41 -0.26167529 0.32684221
42 2.16930887 -0.26167529
43 2.31980220 2.16930887
44 -0.19813557 2.31980220
45 0.01124777 -0.19813557
46 2.19750082 0.01124777
47 2.61437304 2.19750082
48 2.72713082 2.61437304
49 3.90449637 2.72713082
50 4.13182249 3.90449637
51 4.38182249 4.13182249
52 5.87437304 4.38182249
53 4.81725666 5.87437304
54 7.65704693 4.81725666
55 7.98383526 7.65704693
56 4.89141305 7.98383526
57 5.21149443 4.89141305
58 5.24272776 5.21149443
59 4.12618610 5.24272776
60 4.72120582 4.12618610
61 5.04276971 4.72120582
62 3.11133166 5.04276971
63 4.60939526 3.11133166
64 4.85079387 4.60939526
65 3.63569277 4.85079387
66 5.75120331 3.63569277
67 6.93852692 5.75120331
68 3.11898582 6.93852692
69 4.04552749 3.11898582
70 4.03750332 4.04552749
71 4.86030054 4.03750332
> 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/7ujoe1291122100.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/85s5h1291122100.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/95s5h1291122100.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/105s5h1291122100.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/118a3n1291122100.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/12cbkt1291122100.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/131uh51291122100.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/144cxt1291122100.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/15x4xv1291122100.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/16teum1291122100.tab")
+ }
>
> try(system("convert tmp/1r0pq1291122100.ps tmp/1r0pq1291122100.png",intern=TRUE))
character(0)
> try(system("convert tmp/2r0pq1291122100.ps tmp/2r0pq1291122100.png",intern=TRUE))
character(0)
> try(system("convert tmp/3r0pq1291122100.ps tmp/3r0pq1291122100.png",intern=TRUE))
character(0)
> try(system("convert tmp/42r6t1291122100.ps tmp/42r6t1291122100.png",intern=TRUE))
character(0)
> try(system("convert tmp/52r6t1291122100.ps tmp/52r6t1291122100.png",intern=TRUE))
character(0)
> try(system("convert tmp/62r6t1291122100.ps tmp/62r6t1291122100.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ujoe1291122100.ps tmp/7ujoe1291122100.png",intern=TRUE))
character(0)
> try(system("convert tmp/85s5h1291122100.ps tmp/85s5h1291122100.png",intern=TRUE))
character(0)
> try(system("convert tmp/95s5h1291122100.ps tmp/95s5h1291122100.png",intern=TRUE))
character(0)
> try(system("convert tmp/105s5h1291122100.ps tmp/105s5h1291122100.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.574 1.631 6.577