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 = 'Linear Trend'
> par2 = 'Include Monthly 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
Faillissementen CPI M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 627 98.60 1 0 0 0 0 0 0 0 0 0 0 1
2 696 98.97 0 1 0 0 0 0 0 0 0 0 0 2
3 825 99.11 0 0 1 0 0 0 0 0 0 0 0 3
4 677 99.64 0 0 0 1 0 0 0 0 0 0 0 4
5 656 100.03 0 0 0 0 1 0 0 0 0 0 0 5
6 785 99.98 0 0 0 0 0 1 0 0 0 0 0 6
7 412 100.32 0 0 0 0 0 0 1 0 0 0 0 7
8 352 100.44 0 0 0 0 0 0 0 1 0 0 0 8
9 839 100.51 0 0 0 0 0 0 0 0 1 0 0 9
10 729 101.00 0 0 0 0 0 0 0 0 0 1 0 10
11 696 100.88 0 0 0 0 0 0 0 0 0 0 1 11
12 641 100.55 0 0 0 0 0 0 0 0 0 0 0 12
13 695 100.83 1 0 0 0 0 0 0 0 0 0 0 13
14 638 101.51 0 1 0 0 0 0 0 0 0 0 0 14
15 762 102.16 0 0 1 0 0 0 0 0 0 0 0 15
16 635 102.39 0 0 0 1 0 0 0 0 0 0 0 16
17 721 102.54 0 0 0 0 1 0 0 0 0 0 0 17
18 854 102.85 0 0 0 0 0 1 0 0 0 0 0 18
19 418 103.47 0 0 0 0 0 0 1 0 0 0 0 19
20 367 103.57 0 0 0 0 0 0 0 1 0 0 0 20
21 824 103.69 0 0 0 0 0 0 0 0 1 0 0 21
22 687 103.50 0 0 0 0 0 0 0 0 0 1 0 22
23 601 103.47 0 0 0 0 0 0 0 0 0 0 1 23
24 676 103.45 0 0 0 0 0 0 0 0 0 0 0 24
25 740 103.48 1 0 0 0 0 0 0 0 0 0 0 25
26 691 103.93 0 1 0 0 0 0 0 0 0 0 0 26
27 683 103.89 0 0 1 0 0 0 0 0 0 0 0 27
28 594 104.40 0 0 0 1 0 0 0 0 0 0 0 28
29 729 104.79 0 0 0 0 1 0 0 0 0 0 0 29
30 731 104.77 0 0 0 0 0 1 0 0 0 0 0 30
31 386 105.13 0 0 0 0 0 0 1 0 0 0 0 31
32 331 105.26 0 0 0 0 0 0 0 1 0 0 0 32
33 707 104.96 0 0 0 0 0 0 0 0 1 0 0 33
34 715 104.75 0 0 0 0 0 0 0 0 0 1 0 34
35 657 105.01 0 0 0 0 0 0 0 0 0 0 1 35
36 653 105.15 0 0 0 0 0 0 0 0 0 0 0 36
37 642 105.20 1 0 0 0 0 0 0 0 0 0 0 37
38 643 105.77 0 1 0 0 0 0 0 0 0 0 0 38
39 718 105.78 0 0 1 0 0 0 0 0 0 0 0 39
40 654 106.26 0 0 0 1 0 0 0 0 0 0 0 40
41 632 106.13 0 0 0 0 1 0 0 0 0 0 0 41
42 731 106.12 0 0 0 0 0 1 0 0 0 0 0 42
43 392 106.57 0 0 0 0 0 0 1 0 0 0 0 43
44 344 106.44 0 0 0 0 0 0 0 1 0 0 0 44
45 792 106.54 0 0 0 0 0 0 0 0 1 0 0 45
46 852 107.10 0 0 0 0 0 0 0 0 0 1 0 46
47 649 108.10 0 0 0 0 0 0 0 0 0 0 1 47
48 629 108.40 0 0 0 0 0 0 0 0 0 0 0 48
49 685 108.84 1 0 0 0 0 0 0 0 0 0 0 49
50 617 109.62 0 1 0 0 0 0 0 0 0 0 0 50
51 715 110.42 0 0 1 0 0 0 0 0 0 0 0 51
52 715 110.67 0 0 0 1 0 0 0 0 0 0 0 52
53 629 111.66 0 0 0 0 1 0 0 0 0 0 0 53
54 916 112.28 0 0 0 0 0 1 0 0 0 0 0 54
55 531 112.87 0 0 0 0 0 0 1 0 0 0 0 55
56 357 112.18 0 0 0 0 0 0 0 1 0 0 0 56
57 917 112.36 0 0 0 0 0 0 0 0 1 0 0 57
58 828 112.16 0 0 0 0 0 0 0 0 0 1 0 58
59 708 111.49 0 0 0 0 0 0 0 0 0 0 1 59
60 858 111.25 0 0 0 0 0 0 0 0 0 0 0 60
61 775 111.36 1 0 0 0 0 0 0 0 0 0 0 61
62 785 111.74 0 1 0 0 0 0 0 0 0 0 0 62
63 1006 111.10 0 0 1 0 0 0 0 0 0 0 0 63
64 789 111.33 0 0 0 1 0 0 0 0 0 0 0 64
65 734 111.25 0 0 0 0 1 0 0 0 0 0 0 65
66 906 111.04 0 0 0 0 0 1 0 0 0 0 0 66
67 532 110.97 0 0 0 0 0 0 1 0 0 0 0 67
68 387 111.31 0 0 0 0 0 0 0 1 0 0 0 68
69 991 111.02 0 0 0 0 0 0 0 0 1 0 0 69
70 841 111.07 0 0 0 0 0 0 0 0 0 1 0 70
71 892 111.36 0 0 0 0 0 0 0 0 0 0 1 71
72 782 111.54 0 0 0 0 0 0 0 0 0 0 0 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) CPI M1 M2 M3 M4
408.959 2.239 7.331 -10.936 93.826 -15.901
M5 M6 M7 M8 M9 M10
-11.767 123.599 -253.983 -344.163 143.154 71.906
M11 t
-4.594 1.395
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-126.154 -47.664 3.389 47.263 166.587
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 408.959 820.447 0.498 0.620046
CPI 2.239 8.336 0.269 0.789186
M1 7.331 39.440 0.186 0.853197
M2 -10.936 39.613 -0.276 0.783470
M3 93.826 39.535 2.373 0.020968 *
M4 -15.901 39.674 -0.401 0.690048
M5 -11.767 39.750 -0.296 0.768264
M6 123.599 39.612 3.120 0.002815 **
M7 -253.983 39.826 -6.377 3.23e-08 ***
M8 -344.163 39.531 -8.706 4.06e-12 ***
M9 143.154 39.325 3.640 0.000581 ***
M10 71.906 39.247 1.832 0.072066 .
M11 -4.594 39.209 -0.117 0.907129
t 1.395 1.708 0.817 0.417509
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 67.85 on 58 degrees of freedom
Multiple R-squared: 0.8465, Adjusted R-squared: 0.8121
F-statistic: 24.61 on 13 and 58 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.149205322 0.298410644 0.8507947
[2,] 0.253624662 0.507249324 0.7463753
[3,] 0.168289289 0.336578578 0.8317107
[4,] 0.116721696 0.233443392 0.8832783
[5,] 0.067604369 0.135208737 0.9323956
[6,] 0.053292268 0.106584536 0.9467077
[7,] 0.074961143 0.149922286 0.9250389
[8,] 0.054226318 0.108452636 0.9457737
[9,] 0.074480151 0.148960303 0.9255198
[10,] 0.059431617 0.118863233 0.9405684
[11,] 0.131364670 0.262729341 0.8686353
[12,] 0.096331760 0.192663520 0.9036682
[13,] 0.176111035 0.352222070 0.8238890
[14,] 0.159762043 0.319524087 0.8402380
[15,] 0.113219806 0.226439611 0.8867802
[16,] 0.104413308 0.208826616 0.8955867
[17,] 0.105628777 0.211257555 0.8943712
[18,] 0.104793504 0.209587008 0.8952065
[19,] 0.088522114 0.177044227 0.9114779
[20,] 0.064195064 0.128390128 0.9358049
[21,] 0.042157279 0.084314557 0.9578427
[22,] 0.029751612 0.059503224 0.9702484
[23,] 0.018233115 0.036466231 0.9817669
[24,] 0.014167237 0.028334473 0.9858328
[25,] 0.011308022 0.022616043 0.9886920
[26,] 0.007292272 0.014584545 0.9927077
[27,] 0.004479977 0.008959954 0.9955200
[28,] 0.005206125 0.010412250 0.9947939
[29,] 0.003230534 0.006461068 0.9967695
[30,] 0.108452087 0.216904173 0.8915479
[31,] 0.082671167 0.165342333 0.9173288
[32,] 0.053576749 0.107153498 0.9464233
[33,] 0.050800158 0.101600316 0.9491998
[34,] 0.033564430 0.067128860 0.9664356
[35,] 0.173234939 0.346469877 0.8267651
[36,] 0.135835153 0.271670306 0.8641648
[37,] 0.117830331 0.235660662 0.8821697
[38,] 0.112766624 0.225533248 0.8872334
[39,] 0.081155007 0.162310014 0.9188450
> postscript(file="/var/www/html/rcomp/tmp/1m8hk1291138261.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/2m8hk1291138261.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/3fhyn1291138261.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/4fhyn1291138261.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/5fhyn1291138261.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
-11.460957 73.582627 96.112482 55.257956 27.856182 20.206976
7 8 9 10 11 12
22.633606 51.149732 49.281545 8.037638 50.412091 -9.837806
13 14 15 16 17 18
34.810060 -6.840481 9.547425 -9.635367 70.500247 66.044960
19 20 21 22 23 24
4.844638 42.405546 10.425403 -56.295906 -67.122973 1.933005
25 26 27 28 29 30
57.140649 24.005104 -90.062000 -71.871745 56.726482 -77.989897
31 32 33 34 35 36
-47.608050 -14.114315 -126.154030 -47.830556 -31.306967 -41.609248
37 38 39 40 41 42
-61.446386 -44.850624 -76.029685 -32.772255 -60.009689 -97.748459
43 44 45 46 47 48
-61.568132 -20.492228 -61.427589 67.171767 -62.961589 -89.622127
49 50 51 52 53 54
-43.332520 -96.206973 -106.154933 1.617492 -92.127750 56.722838
55 56 57 58 59 60
46.589689 -37.080503 33.805008 15.106090 -28.287944 116.260639
61 62 63 64 65 66
24.289154 50.310347 166.586711 57.403919 -2.945471 32.763582
67 68 69 70 71 72
35.108249 -21.868231 94.069663 13.810967 139.267383 22.875538
> postscript(file="/var/www/html/rcomp/tmp/6pqf81291138261.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 -11.460957 NA
1 73.582627 -11.460957
2 96.112482 73.582627
3 55.257956 96.112482
4 27.856182 55.257956
5 20.206976 27.856182
6 22.633606 20.206976
7 51.149732 22.633606
8 49.281545 51.149732
9 8.037638 49.281545
10 50.412091 8.037638
11 -9.837806 50.412091
12 34.810060 -9.837806
13 -6.840481 34.810060
14 9.547425 -6.840481
15 -9.635367 9.547425
16 70.500247 -9.635367
17 66.044960 70.500247
18 4.844638 66.044960
19 42.405546 4.844638
20 10.425403 42.405546
21 -56.295906 10.425403
22 -67.122973 -56.295906
23 1.933005 -67.122973
24 57.140649 1.933005
25 24.005104 57.140649
26 -90.062000 24.005104
27 -71.871745 -90.062000
28 56.726482 -71.871745
29 -77.989897 56.726482
30 -47.608050 -77.989897
31 -14.114315 -47.608050
32 -126.154030 -14.114315
33 -47.830556 -126.154030
34 -31.306967 -47.830556
35 -41.609248 -31.306967
36 -61.446386 -41.609248
37 -44.850624 -61.446386
38 -76.029685 -44.850624
39 -32.772255 -76.029685
40 -60.009689 -32.772255
41 -97.748459 -60.009689
42 -61.568132 -97.748459
43 -20.492228 -61.568132
44 -61.427589 -20.492228
45 67.171767 -61.427589
46 -62.961589 67.171767
47 -89.622127 -62.961589
48 -43.332520 -89.622127
49 -96.206973 -43.332520
50 -106.154933 -96.206973
51 1.617492 -106.154933
52 -92.127750 1.617492
53 56.722838 -92.127750
54 46.589689 56.722838
55 -37.080503 46.589689
56 33.805008 -37.080503
57 15.106090 33.805008
58 -28.287944 15.106090
59 116.260639 -28.287944
60 24.289154 116.260639
61 50.310347 24.289154
62 166.586711 50.310347
63 57.403919 166.586711
64 -2.945471 57.403919
65 32.763582 -2.945471
66 35.108249 32.763582
67 -21.868231 35.108249
68 94.069663 -21.868231
69 13.810967 94.069663
70 139.267383 13.810967
71 22.875538 139.267383
72 NA 22.875538
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 73.582627 -11.460957
[2,] 96.112482 73.582627
[3,] 55.257956 96.112482
[4,] 27.856182 55.257956
[5,] 20.206976 27.856182
[6,] 22.633606 20.206976
[7,] 51.149732 22.633606
[8,] 49.281545 51.149732
[9,] 8.037638 49.281545
[10,] 50.412091 8.037638
[11,] -9.837806 50.412091
[12,] 34.810060 -9.837806
[13,] -6.840481 34.810060
[14,] 9.547425 -6.840481
[15,] -9.635367 9.547425
[16,] 70.500247 -9.635367
[17,] 66.044960 70.500247
[18,] 4.844638 66.044960
[19,] 42.405546 4.844638
[20,] 10.425403 42.405546
[21,] -56.295906 10.425403
[22,] -67.122973 -56.295906
[23,] 1.933005 -67.122973
[24,] 57.140649 1.933005
[25,] 24.005104 57.140649
[26,] -90.062000 24.005104
[27,] -71.871745 -90.062000
[28,] 56.726482 -71.871745
[29,] -77.989897 56.726482
[30,] -47.608050 -77.989897
[31,] -14.114315 -47.608050
[32,] -126.154030 -14.114315
[33,] -47.830556 -126.154030
[34,] -31.306967 -47.830556
[35,] -41.609248 -31.306967
[36,] -61.446386 -41.609248
[37,] -44.850624 -61.446386
[38,] -76.029685 -44.850624
[39,] -32.772255 -76.029685
[40,] -60.009689 -32.772255
[41,] -97.748459 -60.009689
[42,] -61.568132 -97.748459
[43,] -20.492228 -61.568132
[44,] -61.427589 -20.492228
[45,] 67.171767 -61.427589
[46,] -62.961589 67.171767
[47,] -89.622127 -62.961589
[48,] -43.332520 -89.622127
[49,] -96.206973 -43.332520
[50,] -106.154933 -96.206973
[51,] 1.617492 -106.154933
[52,] -92.127750 1.617492
[53,] 56.722838 -92.127750
[54,] 46.589689 56.722838
[55,] -37.080503 46.589689
[56,] 33.805008 -37.080503
[57,] 15.106090 33.805008
[58,] -28.287944 15.106090
[59,] 116.260639 -28.287944
[60,] 24.289154 116.260639
[61,] 50.310347 24.289154
[62,] 166.586711 50.310347
[63,] 57.403919 166.586711
[64,] -2.945471 57.403919
[65,] 32.763582 -2.945471
[66,] 35.108249 32.763582
[67,] -21.868231 35.108249
[68,] 94.069663 -21.868231
[69,] 13.810967 94.069663
[70,] 139.267383 13.810967
[71,] 22.875538 139.267383
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 73.582627 -11.460957
2 96.112482 73.582627
3 55.257956 96.112482
4 27.856182 55.257956
5 20.206976 27.856182
6 22.633606 20.206976
7 51.149732 22.633606
8 49.281545 51.149732
9 8.037638 49.281545
10 50.412091 8.037638
11 -9.837806 50.412091
12 34.810060 -9.837806
13 -6.840481 34.810060
14 9.547425 -6.840481
15 -9.635367 9.547425
16 70.500247 -9.635367
17 66.044960 70.500247
18 4.844638 66.044960
19 42.405546 4.844638
20 10.425403 42.405546
21 -56.295906 10.425403
22 -67.122973 -56.295906
23 1.933005 -67.122973
24 57.140649 1.933005
25 24.005104 57.140649
26 -90.062000 24.005104
27 -71.871745 -90.062000
28 56.726482 -71.871745
29 -77.989897 56.726482
30 -47.608050 -77.989897
31 -14.114315 -47.608050
32 -126.154030 -14.114315
33 -47.830556 -126.154030
34 -31.306967 -47.830556
35 -41.609248 -31.306967
36 -61.446386 -41.609248
37 -44.850624 -61.446386
38 -76.029685 -44.850624
39 -32.772255 -76.029685
40 -60.009689 -32.772255
41 -97.748459 -60.009689
42 -61.568132 -97.748459
43 -20.492228 -61.568132
44 -61.427589 -20.492228
45 67.171767 -61.427589
46 -62.961589 67.171767
47 -89.622127 -62.961589
48 -43.332520 -89.622127
49 -96.206973 -43.332520
50 -106.154933 -96.206973
51 1.617492 -106.154933
52 -92.127750 1.617492
53 56.722838 -92.127750
54 46.589689 56.722838
55 -37.080503 46.589689
56 33.805008 -37.080503
57 15.106090 33.805008
58 -28.287944 15.106090
59 116.260639 -28.287944
60 24.289154 116.260639
61 50.310347 24.289154
62 166.586711 50.310347
63 57.403919 166.586711
64 -2.945471 57.403919
65 32.763582 -2.945471
66 35.108249 32.763582
67 -21.868231 35.108249
68 94.069663 -21.868231
69 13.810967 94.069663
70 139.267383 13.810967
71 22.875538 139.267383
> 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/700wt1291138261.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/800wt1291138261.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/900wt1291138261.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/10bree1291138261.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/11esc21291138261.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/12iabq1291138261.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/13obp11291138261.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/14auo71291138261.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/15dc5v1291138261.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/16942m1291138261.tab")
+ }
>
> try(system("convert tmp/1m8hk1291138261.ps tmp/1m8hk1291138261.png",intern=TRUE))
character(0)
> try(system("convert tmp/2m8hk1291138261.ps tmp/2m8hk1291138261.png",intern=TRUE))
character(0)
> try(system("convert tmp/3fhyn1291138261.ps tmp/3fhyn1291138261.png",intern=TRUE))
character(0)
> try(system("convert tmp/4fhyn1291138261.ps tmp/4fhyn1291138261.png",intern=TRUE))
character(0)
> try(system("convert tmp/5fhyn1291138261.ps tmp/5fhyn1291138261.png",intern=TRUE))
character(0)
> try(system("convert tmp/6pqf81291138261.ps tmp/6pqf81291138261.png",intern=TRUE))
character(0)
> try(system("convert tmp/700wt1291138261.ps tmp/700wt1291138261.png",intern=TRUE))
character(0)
> try(system("convert tmp/800wt1291138261.ps tmp/800wt1291138261.png",intern=TRUE))
character(0)
> try(system("convert tmp/900wt1291138261.ps tmp/900wt1291138261.png",intern=TRUE))
character(0)
> try(system("convert tmp/10bree1291138261.ps tmp/10bree1291138261.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.593 1.618 6.115