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(13.193,15.234,14.718,16.961,13.945,15.876,16.226,18.316,16.748,17.904,17.209,18.950,17.225,18.710,17.236,18.687,17.580,19.568,17.381,19.580,17.260,18.661,15.658,18.674,15.908,17.475,17.725,19.562,16.368,19.555,17.743,19.867,15.703,19.324,18.162,19.074,15.323,19.704,18.375,18.352,13.927,17.795,16.761,18.902,16.239,19.158,18.279,15.698,16.239,18.431,18.414,19.801,14.995,18.706,18.232,19.409,16.263,19.017,20.298,19.891,15.203,17.845,17.502,18.532,15.737,17.770,17.224,17.601,14.940,18.507,17.635,19.392,15.699,17.661,18.243,19.643,15.770,17.344,17.229,17.322,16.152,17.919,16.918,18.114,16.308,17.759,16.021,17.952,15.954,17.762,16.610,17.751,15.458,18.106,15.990,15.349,13.185,15.409,16.007,16.633,14.800,15.974,15.693),dim=c(1,103),dimnames=list(c('aantal'),1:103))
> y <- array(NA,dim=c(1,103),dimnames=list(c('aantal'),1:103))
> 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 = 'Include Quarterly 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
aantal Q1 Q2 Q3
1 13.193 1 0 0
2 15.234 0 1 0
3 14.718 0 0 1
4 16.961 0 0 0
5 13.945 1 0 0
6 15.876 0 1 0
7 16.226 0 0 1
8 18.316 0 0 0
9 16.748 1 0 0
10 17.904 0 1 0
11 17.209 0 0 1
12 18.950 0 0 0
13 17.225 1 0 0
14 18.710 0 1 0
15 17.236 0 0 1
16 18.687 0 0 0
17 17.580 1 0 0
18 19.568 0 1 0
19 17.381 0 0 1
20 19.580 0 0 0
21 17.260 1 0 0
22 18.661 0 1 0
23 15.658 0 0 1
24 18.674 0 0 0
25 15.908 1 0 0
26 17.475 0 1 0
27 17.725 0 0 1
28 19.562 0 0 0
29 16.368 1 0 0
30 19.555 0 1 0
31 17.743 0 0 1
32 19.867 0 0 0
33 15.703 1 0 0
34 19.324 0 1 0
35 18.162 0 0 1
36 19.074 0 0 0
37 15.323 1 0 0
38 19.704 0 1 0
39 18.375 0 0 1
40 18.352 0 0 0
41 13.927 1 0 0
42 17.795 0 1 0
43 16.761 0 0 1
44 18.902 0 0 0
45 16.239 1 0 0
46 19.158 0 1 0
47 18.279 0 0 1
48 15.698 0 0 0
49 16.239 1 0 0
50 18.431 0 1 0
51 18.414 0 0 1
52 19.801 0 0 0
53 14.995 1 0 0
54 18.706 0 1 0
55 18.232 0 0 1
56 19.409 0 0 0
57 16.263 1 0 0
58 19.017 0 1 0
59 20.298 0 0 1
60 19.891 0 0 0
61 15.203 1 0 0
62 17.845 0 1 0
63 17.502 0 0 1
64 18.532 0 0 0
65 15.737 1 0 0
66 17.770 0 1 0
67 17.224 0 0 1
68 17.601 0 0 0
69 14.940 1 0 0
70 18.507 0 1 0
71 17.635 0 0 1
72 19.392 0 0 0
73 15.699 1 0 0
74 17.661 0 1 0
75 18.243 0 0 1
76 19.643 0 0 0
77 15.770 1 0 0
78 17.344 0 1 0
79 17.229 0 0 1
80 17.322 0 0 0
81 16.152 1 0 0
82 17.919 0 1 0
83 16.918 0 0 1
84 18.114 0 0 0
85 16.308 1 0 0
86 17.759 0 1 0
87 16.021 0 0 1
88 17.952 0 0 0
89 15.954 1 0 0
90 17.762 0 1 0
91 16.610 0 0 1
92 17.751 0 0 0
93 15.458 1 0 0
94 18.106 0 1 0
95 15.990 0 0 1
96 15.349 0 0 0
97 13.185 1 0 0
98 15.409 0 1 0
99 16.007 0 0 1
100 16.633 0 0 0
101 14.800 1 0 0
102 15.974 0 1 0
103 15.693 0 0 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Q1 Q2 Q3
18.4005 -2.7804 -0.4323 -1.1894
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3.0515 -0.6127 0.1169 0.7448 3.0869
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.4005 0.2400 76.655 < 2e-16 ***
Q1 -2.7804 0.3362 -8.270 6.3e-13 ***
Q2 -0.4323 0.3362 -1.286 0.201499
Q3 -1.1894 0.3362 -3.538 0.000616 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.2 on 99 degrees of freedom
Multiple R-squared: 0.4474, Adjusted R-squared: 0.4306
F-statistic: 26.72 on 3 and 99 DF, p-value: 9.622e-13
> 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.2828383 0.56567653 0.71716173
[2,] 0.2814237 0.56284734 0.71857633
[3,] 0.8099903 0.38001947 0.19000973
[4,] 0.8825769 0.23484616 0.11742308
[5,] 0.8853018 0.22939640 0.11469820
[6,] 0.8649436 0.27011271 0.13505635
[7,] 0.9362098 0.12758034 0.06379017
[8,] 0.9595211 0.08095776 0.04047888
[9,] 0.9484977 0.10300465 0.05150232
[10,] 0.9261139 0.14777220 0.07388610
[11,] 0.9581674 0.08366513 0.04183257
[12,] 0.9798485 0.04030292 0.02015146
[13,] 0.9729957 0.05400866 0.02700433
[14,] 0.9707464 0.05850721 0.02925361
[15,] 0.9738735 0.05225298 0.02612649
[16,] 0.9685941 0.06281184 0.03140592
[17,] 0.9664780 0.06704402 0.03352201
[18,] 0.9518823 0.09623541 0.04811771
[19,] 0.9329956 0.13400889 0.06700445
[20,] 0.9104628 0.17907449 0.08953724
[21,] 0.9007935 0.19841300 0.09920650
[22,] 0.8920767 0.21584668 0.10792334
[23,] 0.8669624 0.26607523 0.13303761
[24,] 0.8947397 0.21052064 0.10526032
[25,] 0.8794279 0.24114419 0.12057210
[26,] 0.8835431 0.23291377 0.11645688
[27,] 0.8519808 0.29603848 0.14801924
[28,] 0.8597336 0.28053284 0.14026642
[29,] 0.8551814 0.28963728 0.14481864
[30,] 0.8268144 0.34637115 0.17318557
[31,] 0.7921616 0.41567673 0.20783836
[32,] 0.8291065 0.34178697 0.17089349
[33,] 0.8311804 0.33763916 0.16881958
[34,] 0.7949161 0.41016787 0.20508393
[35,] 0.8367619 0.32647614 0.16323807
[36,] 0.7997314 0.40053712 0.20026856
[37,] 0.7615424 0.47691514 0.23845757
[38,] 0.7234114 0.55317727 0.27658863
[39,] 0.6852700 0.62945993 0.31472996
[40,] 0.6819766 0.63604683 0.31802341
[41,] 0.6724590 0.65508201 0.32754100
[42,] 0.8510531 0.29789380 0.14894690
[43,] 0.8254623 0.34907530 0.17453765
[44,] 0.7938732 0.41225360 0.20612680
[45,] 0.7931520 0.41369592 0.20684796
[46,] 0.8103018 0.37939631 0.18969815
[47,] 0.7796323 0.44073535 0.22036768
[48,] 0.7576954 0.48460926 0.24230463
[49,] 0.7455309 0.50893814 0.25446907
[50,] 0.7400413 0.51991731 0.25995865
[51,] 0.7086669 0.58266613 0.29133307
[52,] 0.7130941 0.57381178 0.28690589
[53,] 0.9370642 0.12587165 0.06293583
[54,] 0.9596659 0.08066811 0.04033406
[55,] 0.9462001 0.10759983 0.05379992
[56,] 0.9303322 0.13933562 0.06966781
[57,] 0.9144002 0.17119968 0.08559984
[58,] 0.8968857 0.20622858 0.10311429
[59,] 0.8685906 0.26281880 0.13140940
[60,] 0.8376571 0.32468574 0.16234287
[61,] 0.8027453 0.39450941 0.19725470
[62,] 0.7699063 0.46018736 0.23009368
[63,] 0.7316472 0.53670558 0.26835279
[64,] 0.7210602 0.55787962 0.27893981
[65,] 0.6963121 0.60737588 0.30368794
[66,] 0.7434384 0.51312325 0.25656163
[67,] 0.6903660 0.61926806 0.30963403
[68,] 0.6409127 0.71817451 0.35908726
[69,] 0.7031306 0.59373877 0.29686939
[70,] 0.8399742 0.32005164 0.16002582
[71,] 0.7993910 0.40121792 0.20060896
[72,] 0.7532863 0.49342734 0.24671367
[73,] 0.7324826 0.53503487 0.26751743
[74,] 0.6883443 0.62331148 0.31165574
[75,] 0.6622743 0.67545149 0.33772574
[76,] 0.6262320 0.74753600 0.37376800
[77,] 0.5856417 0.82871652 0.41435826
[78,] 0.5728493 0.85430144 0.42715072
[79,] 0.5939794 0.81204128 0.40602064
[80,] 0.5554646 0.88907083 0.44453542
[81,] 0.4856755 0.97135099 0.51432451
[82,] 0.4901026 0.98020513 0.50989743
[83,] 0.5079500 0.98409997 0.49204999
[84,] 0.4989261 0.99785219 0.50107391
[85,] 0.4284041 0.85680822 0.57159589
[86,] 0.4945915 0.98918310 0.50540845
[87,] 0.5110605 0.97787902 0.48893951
[88,] 0.8451184 0.30976323 0.15488161
[89,] 0.7368605 0.52627906 0.26313953
[90,] 0.7554254 0.48914918 0.24457459
> postscript(file="/var/www/html/rcomp/tmp/18z7c1292941512.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/21rof1292941512.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/31rof1292941512.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/41rof1292941512.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/5ci601292941512.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 = 103
Frequency = 1
1 2 3 4 5 6
-2.427076923 -2.734230769 -2.493115385 -1.439520000 -1.675076923 -2.092230769
7 8 9 10 11 12
-0.985115385 -0.084520000 1.127923077 -0.064230769 -0.002115385 0.549480000
13 14 15 16 17 18
1.604923077 0.741769231 0.024884615 0.286480000 1.959923077 1.599769231
19 20 21 22 23 24
0.169884615 1.179480000 1.639923077 0.692769231 -1.553115385 0.273480000
25 26 27 28 29 30
0.287923077 -0.493230769 0.513884615 1.161480000 0.747923077 1.586769231
31 32 33 34 35 36
0.531884615 1.466480000 0.082923077 1.355769231 0.950884615 0.673480000
37 38 39 40 41 42
-0.297076923 1.735769231 1.163884615 -0.048520000 -1.693076923 -0.173230769
43 44 45 46 47 48
-0.450115385 0.501480000 0.618923077 1.189769231 1.067884615 -2.702520000
49 50 51 52 53 54
0.618923077 0.462769231 1.202884615 1.400480000 -0.625076923 0.737769231
55 56 57 58 59 60
1.020884615 1.008480000 0.642923077 1.048769231 3.086884615 1.490480000
61 62 63 64 65 66
-0.417076923 -0.123230769 0.290884615 0.131480000 0.116923077 -0.198230769
67 68 69 70 71 72
0.012884615 -0.799520000 -0.680076923 0.538769231 0.423884615 0.991480000
73 74 75 76 77 78
0.078923077 -0.307230769 1.031884615 1.242480000 0.149923077 -0.624230769
79 80 81 82 83 84
0.017884615 -1.078520000 0.531923077 -0.049230769 -0.293115385 -0.286520000
85 86 87 88 89 90
0.687923077 -0.209230769 -1.190115385 -0.448520000 0.333923077 -0.206230769
91 92 93 94 95 96
-0.601115385 -0.649520000 -0.162076923 0.137769231 -1.221115385 -3.051520000
97 98 99 100 101 102
-2.435076923 -2.559230769 -1.204115385 -1.767520000 -0.820076923 -1.994230769
103
-1.518115385
> postscript(file="/var/www/html/rcomp/tmp/6ci601292941512.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 = 103
Frequency = 1
lag(myerror, k = 1) myerror
0 -2.427076923 NA
1 -2.734230769 -2.427076923
2 -2.493115385 -2.734230769
3 -1.439520000 -2.493115385
4 -1.675076923 -1.439520000
5 -2.092230769 -1.675076923
6 -0.985115385 -2.092230769
7 -0.084520000 -0.985115385
8 1.127923077 -0.084520000
9 -0.064230769 1.127923077
10 -0.002115385 -0.064230769
11 0.549480000 -0.002115385
12 1.604923077 0.549480000
13 0.741769231 1.604923077
14 0.024884615 0.741769231
15 0.286480000 0.024884615
16 1.959923077 0.286480000
17 1.599769231 1.959923077
18 0.169884615 1.599769231
19 1.179480000 0.169884615
20 1.639923077 1.179480000
21 0.692769231 1.639923077
22 -1.553115385 0.692769231
23 0.273480000 -1.553115385
24 0.287923077 0.273480000
25 -0.493230769 0.287923077
26 0.513884615 -0.493230769
27 1.161480000 0.513884615
28 0.747923077 1.161480000
29 1.586769231 0.747923077
30 0.531884615 1.586769231
31 1.466480000 0.531884615
32 0.082923077 1.466480000
33 1.355769231 0.082923077
34 0.950884615 1.355769231
35 0.673480000 0.950884615
36 -0.297076923 0.673480000
37 1.735769231 -0.297076923
38 1.163884615 1.735769231
39 -0.048520000 1.163884615
40 -1.693076923 -0.048520000
41 -0.173230769 -1.693076923
42 -0.450115385 -0.173230769
43 0.501480000 -0.450115385
44 0.618923077 0.501480000
45 1.189769231 0.618923077
46 1.067884615 1.189769231
47 -2.702520000 1.067884615
48 0.618923077 -2.702520000
49 0.462769231 0.618923077
50 1.202884615 0.462769231
51 1.400480000 1.202884615
52 -0.625076923 1.400480000
53 0.737769231 -0.625076923
54 1.020884615 0.737769231
55 1.008480000 1.020884615
56 0.642923077 1.008480000
57 1.048769231 0.642923077
58 3.086884615 1.048769231
59 1.490480000 3.086884615
60 -0.417076923 1.490480000
61 -0.123230769 -0.417076923
62 0.290884615 -0.123230769
63 0.131480000 0.290884615
64 0.116923077 0.131480000
65 -0.198230769 0.116923077
66 0.012884615 -0.198230769
67 -0.799520000 0.012884615
68 -0.680076923 -0.799520000
69 0.538769231 -0.680076923
70 0.423884615 0.538769231
71 0.991480000 0.423884615
72 0.078923077 0.991480000
73 -0.307230769 0.078923077
74 1.031884615 -0.307230769
75 1.242480000 1.031884615
76 0.149923077 1.242480000
77 -0.624230769 0.149923077
78 0.017884615 -0.624230769
79 -1.078520000 0.017884615
80 0.531923077 -1.078520000
81 -0.049230769 0.531923077
82 -0.293115385 -0.049230769
83 -0.286520000 -0.293115385
84 0.687923077 -0.286520000
85 -0.209230769 0.687923077
86 -1.190115385 -0.209230769
87 -0.448520000 -1.190115385
88 0.333923077 -0.448520000
89 -0.206230769 0.333923077
90 -0.601115385 -0.206230769
91 -0.649520000 -0.601115385
92 -0.162076923 -0.649520000
93 0.137769231 -0.162076923
94 -1.221115385 0.137769231
95 -3.051520000 -1.221115385
96 -2.435076923 -3.051520000
97 -2.559230769 -2.435076923
98 -1.204115385 -2.559230769
99 -1.767520000 -1.204115385
100 -0.820076923 -1.767520000
101 -1.994230769 -0.820076923
102 -1.518115385 -1.994230769
103 NA -1.518115385
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.734230769 -2.427076923
[2,] -2.493115385 -2.734230769
[3,] -1.439520000 -2.493115385
[4,] -1.675076923 -1.439520000
[5,] -2.092230769 -1.675076923
[6,] -0.985115385 -2.092230769
[7,] -0.084520000 -0.985115385
[8,] 1.127923077 -0.084520000
[9,] -0.064230769 1.127923077
[10,] -0.002115385 -0.064230769
[11,] 0.549480000 -0.002115385
[12,] 1.604923077 0.549480000
[13,] 0.741769231 1.604923077
[14,] 0.024884615 0.741769231
[15,] 0.286480000 0.024884615
[16,] 1.959923077 0.286480000
[17,] 1.599769231 1.959923077
[18,] 0.169884615 1.599769231
[19,] 1.179480000 0.169884615
[20,] 1.639923077 1.179480000
[21,] 0.692769231 1.639923077
[22,] -1.553115385 0.692769231
[23,] 0.273480000 -1.553115385
[24,] 0.287923077 0.273480000
[25,] -0.493230769 0.287923077
[26,] 0.513884615 -0.493230769
[27,] 1.161480000 0.513884615
[28,] 0.747923077 1.161480000
[29,] 1.586769231 0.747923077
[30,] 0.531884615 1.586769231
[31,] 1.466480000 0.531884615
[32,] 0.082923077 1.466480000
[33,] 1.355769231 0.082923077
[34,] 0.950884615 1.355769231
[35,] 0.673480000 0.950884615
[36,] -0.297076923 0.673480000
[37,] 1.735769231 -0.297076923
[38,] 1.163884615 1.735769231
[39,] -0.048520000 1.163884615
[40,] -1.693076923 -0.048520000
[41,] -0.173230769 -1.693076923
[42,] -0.450115385 -0.173230769
[43,] 0.501480000 -0.450115385
[44,] 0.618923077 0.501480000
[45,] 1.189769231 0.618923077
[46,] 1.067884615 1.189769231
[47,] -2.702520000 1.067884615
[48,] 0.618923077 -2.702520000
[49,] 0.462769231 0.618923077
[50,] 1.202884615 0.462769231
[51,] 1.400480000 1.202884615
[52,] -0.625076923 1.400480000
[53,] 0.737769231 -0.625076923
[54,] 1.020884615 0.737769231
[55,] 1.008480000 1.020884615
[56,] 0.642923077 1.008480000
[57,] 1.048769231 0.642923077
[58,] 3.086884615 1.048769231
[59,] 1.490480000 3.086884615
[60,] -0.417076923 1.490480000
[61,] -0.123230769 -0.417076923
[62,] 0.290884615 -0.123230769
[63,] 0.131480000 0.290884615
[64,] 0.116923077 0.131480000
[65,] -0.198230769 0.116923077
[66,] 0.012884615 -0.198230769
[67,] -0.799520000 0.012884615
[68,] -0.680076923 -0.799520000
[69,] 0.538769231 -0.680076923
[70,] 0.423884615 0.538769231
[71,] 0.991480000 0.423884615
[72,] 0.078923077 0.991480000
[73,] -0.307230769 0.078923077
[74,] 1.031884615 -0.307230769
[75,] 1.242480000 1.031884615
[76,] 0.149923077 1.242480000
[77,] -0.624230769 0.149923077
[78,] 0.017884615 -0.624230769
[79,] -1.078520000 0.017884615
[80,] 0.531923077 -1.078520000
[81,] -0.049230769 0.531923077
[82,] -0.293115385 -0.049230769
[83,] -0.286520000 -0.293115385
[84,] 0.687923077 -0.286520000
[85,] -0.209230769 0.687923077
[86,] -1.190115385 -0.209230769
[87,] -0.448520000 -1.190115385
[88,] 0.333923077 -0.448520000
[89,] -0.206230769 0.333923077
[90,] -0.601115385 -0.206230769
[91,] -0.649520000 -0.601115385
[92,] -0.162076923 -0.649520000
[93,] 0.137769231 -0.162076923
[94,] -1.221115385 0.137769231
[95,] -3.051520000 -1.221115385
[96,] -2.435076923 -3.051520000
[97,] -2.559230769 -2.435076923
[98,] -1.204115385 -2.559230769
[99,] -1.767520000 -1.204115385
[100,] -0.820076923 -1.767520000
[101,] -1.994230769 -0.820076923
[102,] -1.518115385 -1.994230769
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.734230769 -2.427076923
2 -2.493115385 -2.734230769
3 -1.439520000 -2.493115385
4 -1.675076923 -1.439520000
5 -2.092230769 -1.675076923
6 -0.985115385 -2.092230769
7 -0.084520000 -0.985115385
8 1.127923077 -0.084520000
9 -0.064230769 1.127923077
10 -0.002115385 -0.064230769
11 0.549480000 -0.002115385
12 1.604923077 0.549480000
13 0.741769231 1.604923077
14 0.024884615 0.741769231
15 0.286480000 0.024884615
16 1.959923077 0.286480000
17 1.599769231 1.959923077
18 0.169884615 1.599769231
19 1.179480000 0.169884615
20 1.639923077 1.179480000
21 0.692769231 1.639923077
22 -1.553115385 0.692769231
23 0.273480000 -1.553115385
24 0.287923077 0.273480000
25 -0.493230769 0.287923077
26 0.513884615 -0.493230769
27 1.161480000 0.513884615
28 0.747923077 1.161480000
29 1.586769231 0.747923077
30 0.531884615 1.586769231
31 1.466480000 0.531884615
32 0.082923077 1.466480000
33 1.355769231 0.082923077
34 0.950884615 1.355769231
35 0.673480000 0.950884615
36 -0.297076923 0.673480000
37 1.735769231 -0.297076923
38 1.163884615 1.735769231
39 -0.048520000 1.163884615
40 -1.693076923 -0.048520000
41 -0.173230769 -1.693076923
42 -0.450115385 -0.173230769
43 0.501480000 -0.450115385
44 0.618923077 0.501480000
45 1.189769231 0.618923077
46 1.067884615 1.189769231
47 -2.702520000 1.067884615
48 0.618923077 -2.702520000
49 0.462769231 0.618923077
50 1.202884615 0.462769231
51 1.400480000 1.202884615
52 -0.625076923 1.400480000
53 0.737769231 -0.625076923
54 1.020884615 0.737769231
55 1.008480000 1.020884615
56 0.642923077 1.008480000
57 1.048769231 0.642923077
58 3.086884615 1.048769231
59 1.490480000 3.086884615
60 -0.417076923 1.490480000
61 -0.123230769 -0.417076923
62 0.290884615 -0.123230769
63 0.131480000 0.290884615
64 0.116923077 0.131480000
65 -0.198230769 0.116923077
66 0.012884615 -0.198230769
67 -0.799520000 0.012884615
68 -0.680076923 -0.799520000
69 0.538769231 -0.680076923
70 0.423884615 0.538769231
71 0.991480000 0.423884615
72 0.078923077 0.991480000
73 -0.307230769 0.078923077
74 1.031884615 -0.307230769
75 1.242480000 1.031884615
76 0.149923077 1.242480000
77 -0.624230769 0.149923077
78 0.017884615 -0.624230769
79 -1.078520000 0.017884615
80 0.531923077 -1.078520000
81 -0.049230769 0.531923077
82 -0.293115385 -0.049230769
83 -0.286520000 -0.293115385
84 0.687923077 -0.286520000
85 -0.209230769 0.687923077
86 -1.190115385 -0.209230769
87 -0.448520000 -1.190115385
88 0.333923077 -0.448520000
89 -0.206230769 0.333923077
90 -0.601115385 -0.206230769
91 -0.649520000 -0.601115385
92 -0.162076923 -0.649520000
93 0.137769231 -0.162076923
94 -1.221115385 0.137769231
95 -3.051520000 -1.221115385
96 -2.435076923 -3.051520000
97 -2.559230769 -2.435076923
98 -1.204115385 -2.559230769
99 -1.767520000 -1.204115385
100 -0.820076923 -1.767520000
101 -1.994230769 -0.820076923
102 -1.518115385 -1.994230769
> 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/7m95k1292941512.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/8m95k1292941512.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/9fjmn1292941512.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/10fjmn1292941512.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/1101kt1292941512.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/12m21h1292941512.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/135fof1292941512.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/143ufw1292941512.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/15pcwk1292941512.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/16sduq1292941512.tab")
+ }
>
> try(system("convert tmp/18z7c1292941512.ps tmp/18z7c1292941512.png",intern=TRUE))
character(0)
> try(system("convert tmp/21rof1292941512.ps tmp/21rof1292941512.png",intern=TRUE))
character(0)
> try(system("convert tmp/31rof1292941512.ps tmp/31rof1292941512.png",intern=TRUE))
character(0)
> try(system("convert tmp/41rof1292941512.ps tmp/41rof1292941512.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ci601292941512.ps tmp/5ci601292941512.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ci601292941512.ps tmp/6ci601292941512.png",intern=TRUE))
character(0)
> try(system("convert tmp/7m95k1292941512.ps tmp/7m95k1292941512.png",intern=TRUE))
character(0)
> try(system("convert tmp/8m95k1292941512.ps tmp/8m95k1292941512.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fjmn1292941512.ps tmp/9fjmn1292941512.png",intern=TRUE))
character(0)
> try(system("convert tmp/10fjmn1292941512.ps tmp/10fjmn1292941512.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.035 1.686 14.845