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(267413,294912,267366,293488,264777,290555,258863,284736,254844,281818,254868,287854,277267,316263,285351,325412,286602,326011,283042,328282,276687,317480,277915,317539,277128,313737,277103,312276,275037,309391,270150,302950,267140,300316,264993,304035,287259,333476,291186,337698,292300,335932,288186,323931,281477,313927,282656,314485,280190,313218,280408,309664,276836,302963,275216,298989,274352,298423,271311,301631,289802,329765,290726,335083,292300,327616,278506,309119,269826,295916,265861,291413,269034,291542,264176,284678,255198,276475,253353,272566,246057,264981,235372,263290,258556,296806,260993,303598,254663,286994,250643,276427,243422,266424,247105,267153,248541,268381,245039,262522,237080,255542,237085,253158,225554,243803,226839,250741,247934,280445,248333,285257,246969,270976,245098,261076,246263,255603,255765,260376,264319,263903,268347,264291,273046,263276,273963,262572,267430,256167,271993,264221,292710,293860),dim=c(2,67),dimnames=list(c('Y','X'),1:67))
> y <- array(NA,dim=c(2,67),dimnames=list(c('Y','X'),1:67))
> 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 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 267413 294912 1 0 0 0 0 0 0 0 0 0 0 1
2 267366 293488 0 1 0 0 0 0 0 0 0 0 0 2
3 264777 290555 0 0 1 0 0 0 0 0 0 0 0 3
4 258863 284736 0 0 0 1 0 0 0 0 0 0 0 4
5 254844 281818 0 0 0 0 1 0 0 0 0 0 0 5
6 254868 287854 0 0 0 0 0 1 0 0 0 0 0 6
7 277267 316263 0 0 0 0 0 0 1 0 0 0 0 7
8 285351 325412 0 0 0 0 0 0 0 1 0 0 0 8
9 286602 326011 0 0 0 0 0 0 0 0 1 0 0 9
10 283042 328282 0 0 0 0 0 0 0 0 0 1 0 10
11 276687 317480 0 0 0 0 0 0 0 0 0 0 1 11
12 277915 317539 0 0 0 0 0 0 0 0 0 0 0 12
13 277128 313737 1 0 0 0 0 0 0 0 0 0 0 13
14 277103 312276 0 1 0 0 0 0 0 0 0 0 0 14
15 275037 309391 0 0 1 0 0 0 0 0 0 0 0 15
16 270150 302950 0 0 0 1 0 0 0 0 0 0 0 16
17 267140 300316 0 0 0 0 1 0 0 0 0 0 0 17
18 264993 304035 0 0 0 0 0 1 0 0 0 0 0 18
19 287259 333476 0 0 0 0 0 0 1 0 0 0 0 19
20 291186 337698 0 0 0 0 0 0 0 1 0 0 0 20
21 292300 335932 0 0 0 0 0 0 0 0 1 0 0 21
22 288186 323931 0 0 0 0 0 0 0 0 0 1 0 22
23 281477 313927 0 0 0 0 0 0 0 0 0 0 1 23
24 282656 314485 0 0 0 0 0 0 0 0 0 0 0 24
25 280190 313218 1 0 0 0 0 0 0 0 0 0 0 25
26 280408 309664 0 1 0 0 0 0 0 0 0 0 0 26
27 276836 302963 0 0 1 0 0 0 0 0 0 0 0 27
28 275216 298989 0 0 0 1 0 0 0 0 0 0 0 28
29 274352 298423 0 0 0 0 1 0 0 0 0 0 0 29
30 271311 301631 0 0 0 0 0 1 0 0 0 0 0 30
31 289802 329765 0 0 0 0 0 0 1 0 0 0 0 31
32 290726 335083 0 0 0 0 0 0 0 1 0 0 0 32
33 292300 327616 0 0 0 0 0 0 0 0 1 0 0 33
34 278506 309119 0 0 0 0 0 0 0 0 0 1 0 34
35 269826 295916 0 0 0 0 0 0 0 0 0 0 1 35
36 265861 291413 0 0 0 0 0 0 0 0 0 0 0 36
37 269034 291542 1 0 0 0 0 0 0 0 0 0 0 37
38 264176 284678 0 1 0 0 0 0 0 0 0 0 0 38
39 255198 276475 0 0 1 0 0 0 0 0 0 0 0 39
40 253353 272566 0 0 0 1 0 0 0 0 0 0 0 40
41 246057 264981 0 0 0 0 1 0 0 0 0 0 0 41
42 235372 263290 0 0 0 0 0 1 0 0 0 0 0 42
43 258556 296806 0 0 0 0 0 0 1 0 0 0 0 43
44 260993 303598 0 0 0 0 0 0 0 1 0 0 0 44
45 254663 286994 0 0 0 0 0 0 0 0 1 0 0 45
46 250643 276427 0 0 0 0 0 0 0 0 0 1 0 46
47 243422 266424 0 0 0 0 0 0 0 0 0 0 1 47
48 247105 267153 0 0 0 0 0 0 0 0 0 0 0 48
49 248541 268381 1 0 0 0 0 0 0 0 0 0 0 49
50 245039 262522 0 1 0 0 0 0 0 0 0 0 0 50
51 237080 255542 0 0 1 0 0 0 0 0 0 0 0 51
52 237085 253158 0 0 0 1 0 0 0 0 0 0 0 52
53 225554 243803 0 0 0 0 1 0 0 0 0 0 0 53
54 226839 250741 0 0 0 0 0 1 0 0 0 0 0 54
55 247934 280445 0 0 0 0 0 0 1 0 0 0 0 55
56 248333 285257 0 0 0 0 0 0 0 1 0 0 0 56
57 246969 270976 0 0 0 0 0 0 0 0 1 0 0 57
58 245098 261076 0 0 0 0 0 0 0 0 0 1 0 58
59 246263 255603 0 0 0 0 0 0 0 0 0 0 1 59
60 255765 260376 0 0 0 0 0 0 0 0 0 0 0 60
61 264319 263903 1 0 0 0 0 0 0 0 0 0 0 61
62 268347 264291 0 1 0 0 0 0 0 0 0 0 0 62
63 273046 263276 0 0 1 0 0 0 0 0 0 0 0 63
64 273963 262572 0 0 0 1 0 0 0 0 0 0 0 64
65 267430 256167 0 0 0 0 1 0 0 0 0 0 0 65
66 271993 264221 0 0 0 0 0 1 0 0 0 0 0 66
67 292710 293860 0 0 0 0 0 0 1 0 0 0 0 67
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
-19851.073 0.928 3488.774 5238.899 5813.703 6726.842
M5 M6 M7 M8 M9 M10
5285.694 -899.207 -7657.433 -13975.701 -7847.906 -4738.037
M11 t
-1569.558 455.913
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-10293 -4653 -1886 4595 17457
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.985e+04 2.499e+04 -0.794 0.4305
X 9.280e-01 7.641e-02 12.145 < 2e-16 ***
M1 3.489e+03 4.731e+03 0.737 0.4641
M2 5.239e+03 4.742e+03 1.105 0.2743
M3 5.814e+03 4.778e+03 1.217 0.2291
M4 6.727e+03 4.818e+03 1.396 0.1685
M5 5.286e+03 4.889e+03 1.081 0.2845
M6 -8.992e+02 4.798e+03 -0.187 0.8520
M7 -7.657e+03 4.938e+03 -1.551 0.1269
M8 -1.398e+04 5.248e+03 -2.663 0.0102 *
M9 -7.848e+03 5.088e+03 -1.542 0.1289
M10 -4.738e+03 4.962e+03 -0.955 0.3440
M11 -1.570e+03 4.927e+03 -0.319 0.7513
t 4.559e+02 8.538e+01 5.340 1.98e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7788 on 53 degrees of freedom
Multiple R-squared: 0.8307, Adjusted R-squared: 0.7892
F-statistic: 20.01 on 13 and 53 DF, p-value: 5.95e-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.007262019 1.452404e-02 9.927380e-01
[2,] 0.005546337 1.109267e-02 9.944537e-01
[3,] 0.001750749 3.501499e-03 9.982493e-01
[4,] 0.007060267 1.412053e-02 9.929397e-01
[5,] 0.002505557 5.011114e-03 9.974944e-01
[6,] 0.017547102 3.509420e-02 9.824529e-01
[7,] 0.008424222 1.684844e-02 9.915758e-01
[8,] 0.003408026 6.816051e-03 9.965920e-01
[9,] 0.004353054 8.706107e-03 9.956469e-01
[10,] 0.002256327 4.512653e-03 9.977437e-01
[11,] 0.001047448 2.094896e-03 9.989526e-01
[12,] 0.001480809 2.961618e-03 9.985192e-01
[13,] 0.007425248 1.485050e-02 9.925748e-01
[14,] 0.006003241 1.200648e-02 9.939968e-01
[15,] 0.004181185 8.362370e-03 9.958188e-01
[16,] 0.038755030 7.751006e-02 9.612450e-01
[17,] 0.113675972 2.273519e-01 8.863240e-01
[18,] 0.552578638 8.948427e-01 4.474214e-01
[19,] 0.798966973 4.020661e-01 2.010330e-01
[20,] 0.932690654 1.346187e-01 6.730935e-02
[21,] 0.942016082 1.159678e-01 5.798392e-02
[22,] 0.936026294 1.279474e-01 6.397371e-02
[23,] 0.954289156 9.142169e-02 4.571084e-02
[24,] 0.943208974 1.135821e-01 5.679103e-02
[25,] 0.949441181 1.011176e-01 5.055882e-02
[26,] 0.999722884 5.542321e-04 2.771160e-04
[27,] 0.999951529 9.694270e-05 4.847135e-05
[28,] 0.999963656 7.268766e-05 3.634383e-05
[29,] 0.999918092 1.638159e-04 8.190796e-05
[30,] 0.999651688 6.966233e-04 3.483116e-04
[31,] 0.998504598 2.990805e-03 1.495402e-03
[32,] 0.997196533 5.606935e-03 2.803467e-03
[33,] 0.996474420 7.051161e-03 3.525580e-03
[34,] 0.987433638 2.513272e-02 1.256636e-02
> postscript(file="/var/www/html/rcomp/tmp/1kau81259606617.ps",horizontal=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/2sxbi1259606617.ps",horizontal=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/3yt691259606617.ps",horizontal=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/4bds71259606617.ps",horizontal=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/57a7w1259606617.ps",horizontal=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 = 67
Frequency = 1
1 2 3 4 5 6
9641.7577 8710.1885 7812.2876 5929.2536 5603.3850 5754.9797
7 8 9 10 11 12
8092.8083 13548.9122 7660.3342 -1572.9314 -1528.0936 -2380.3166
13 14 15 16 17 18
-3583.7572 -4458.9905 -4878.4353 -5157.2548 -4737.6747 -4606.9096
19 20 21 22 23 24
-3359.7745 2511.5736 -1319.2901 2137.8261 1088.1218 -276.1720
25 26 27 28 29 30
-5511.0865 -4201.0209 -2585.2268 -1886.4164 -1239.9353 -1528.9634
31 32 33 34 35 36
-2843.9354 -992.6728 926.9779 732.3664 680.3265 -1131.3715
37 38 39 40 41 42
-2022.7706 -2717.0329 -5113.3864 -4699.8959 -3971.7996 -7358.5675
43 44 45 46 47 48
-8975.0226 -6978.6284 -4483.8637 -2263.4961 -3826.1284 -2845.1098
49 50 51 52 53 54
-6493.3783 -6764.2782 -9276.5728 -8428.2785 -10292.6265 -9717.0857
55 56 57 58 59 60
-9885.0140 -8089.1845 -2784.1582 966.2350 3585.7736 6632.9699
61 62 63 64 65 66
7969.2349 9431.1340 14041.3337 14242.5920 14638.6511 17456.5466
67
16970.9382
> postscript(file="/var/www/html/rcomp/tmp/60elt1259606617.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 9641.7577 NA
1 8710.1885 9641.7577
2 7812.2876 8710.1885
3 5929.2536 7812.2876
4 5603.3850 5929.2536
5 5754.9797 5603.3850
6 8092.8083 5754.9797
7 13548.9122 8092.8083
8 7660.3342 13548.9122
9 -1572.9314 7660.3342
10 -1528.0936 -1572.9314
11 -2380.3166 -1528.0936
12 -3583.7572 -2380.3166
13 -4458.9905 -3583.7572
14 -4878.4353 -4458.9905
15 -5157.2548 -4878.4353
16 -4737.6747 -5157.2548
17 -4606.9096 -4737.6747
18 -3359.7745 -4606.9096
19 2511.5736 -3359.7745
20 -1319.2901 2511.5736
21 2137.8261 -1319.2901
22 1088.1218 2137.8261
23 -276.1720 1088.1218
24 -5511.0865 -276.1720
25 -4201.0209 -5511.0865
26 -2585.2268 -4201.0209
27 -1886.4164 -2585.2268
28 -1239.9353 -1886.4164
29 -1528.9634 -1239.9353
30 -2843.9354 -1528.9634
31 -992.6728 -2843.9354
32 926.9779 -992.6728
33 732.3664 926.9779
34 680.3265 732.3664
35 -1131.3715 680.3265
36 -2022.7706 -1131.3715
37 -2717.0329 -2022.7706
38 -5113.3864 -2717.0329
39 -4699.8959 -5113.3864
40 -3971.7996 -4699.8959
41 -7358.5675 -3971.7996
42 -8975.0226 -7358.5675
43 -6978.6284 -8975.0226
44 -4483.8637 -6978.6284
45 -2263.4961 -4483.8637
46 -3826.1284 -2263.4961
47 -2845.1098 -3826.1284
48 -6493.3783 -2845.1098
49 -6764.2782 -6493.3783
50 -9276.5728 -6764.2782
51 -8428.2785 -9276.5728
52 -10292.6265 -8428.2785
53 -9717.0857 -10292.6265
54 -9885.0140 -9717.0857
55 -8089.1845 -9885.0140
56 -2784.1582 -8089.1845
57 966.2350 -2784.1582
58 3585.7736 966.2350
59 6632.9699 3585.7736
60 7969.2349 6632.9699
61 9431.1340 7969.2349
62 14041.3337 9431.1340
63 14242.5920 14041.3337
64 14638.6511 14242.5920
65 17456.5466 14638.6511
66 16970.9382 17456.5466
67 NA 16970.9382
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 8710.1885 9641.7577
[2,] 7812.2876 8710.1885
[3,] 5929.2536 7812.2876
[4,] 5603.3850 5929.2536
[5,] 5754.9797 5603.3850
[6,] 8092.8083 5754.9797
[7,] 13548.9122 8092.8083
[8,] 7660.3342 13548.9122
[9,] -1572.9314 7660.3342
[10,] -1528.0936 -1572.9314
[11,] -2380.3166 -1528.0936
[12,] -3583.7572 -2380.3166
[13,] -4458.9905 -3583.7572
[14,] -4878.4353 -4458.9905
[15,] -5157.2548 -4878.4353
[16,] -4737.6747 -5157.2548
[17,] -4606.9096 -4737.6747
[18,] -3359.7745 -4606.9096
[19,] 2511.5736 -3359.7745
[20,] -1319.2901 2511.5736
[21,] 2137.8261 -1319.2901
[22,] 1088.1218 2137.8261
[23,] -276.1720 1088.1218
[24,] -5511.0865 -276.1720
[25,] -4201.0209 -5511.0865
[26,] -2585.2268 -4201.0209
[27,] -1886.4164 -2585.2268
[28,] -1239.9353 -1886.4164
[29,] -1528.9634 -1239.9353
[30,] -2843.9354 -1528.9634
[31,] -992.6728 -2843.9354
[32,] 926.9779 -992.6728
[33,] 732.3664 926.9779
[34,] 680.3265 732.3664
[35,] -1131.3715 680.3265
[36,] -2022.7706 -1131.3715
[37,] -2717.0329 -2022.7706
[38,] -5113.3864 -2717.0329
[39,] -4699.8959 -5113.3864
[40,] -3971.7996 -4699.8959
[41,] -7358.5675 -3971.7996
[42,] -8975.0226 -7358.5675
[43,] -6978.6284 -8975.0226
[44,] -4483.8637 -6978.6284
[45,] -2263.4961 -4483.8637
[46,] -3826.1284 -2263.4961
[47,] -2845.1098 -3826.1284
[48,] -6493.3783 -2845.1098
[49,] -6764.2782 -6493.3783
[50,] -9276.5728 -6764.2782
[51,] -8428.2785 -9276.5728
[52,] -10292.6265 -8428.2785
[53,] -9717.0857 -10292.6265
[54,] -9885.0140 -9717.0857
[55,] -8089.1845 -9885.0140
[56,] -2784.1582 -8089.1845
[57,] 966.2350 -2784.1582
[58,] 3585.7736 966.2350
[59,] 6632.9699 3585.7736
[60,] 7969.2349 6632.9699
[61,] 9431.1340 7969.2349
[62,] 14041.3337 9431.1340
[63,] 14242.5920 14041.3337
[64,] 14638.6511 14242.5920
[65,] 17456.5466 14638.6511
[66,] 16970.9382 17456.5466
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 8710.1885 9641.7577
2 7812.2876 8710.1885
3 5929.2536 7812.2876
4 5603.3850 5929.2536
5 5754.9797 5603.3850
6 8092.8083 5754.9797
7 13548.9122 8092.8083
8 7660.3342 13548.9122
9 -1572.9314 7660.3342
10 -1528.0936 -1572.9314
11 -2380.3166 -1528.0936
12 -3583.7572 -2380.3166
13 -4458.9905 -3583.7572
14 -4878.4353 -4458.9905
15 -5157.2548 -4878.4353
16 -4737.6747 -5157.2548
17 -4606.9096 -4737.6747
18 -3359.7745 -4606.9096
19 2511.5736 -3359.7745
20 -1319.2901 2511.5736
21 2137.8261 -1319.2901
22 1088.1218 2137.8261
23 -276.1720 1088.1218
24 -5511.0865 -276.1720
25 -4201.0209 -5511.0865
26 -2585.2268 -4201.0209
27 -1886.4164 -2585.2268
28 -1239.9353 -1886.4164
29 -1528.9634 -1239.9353
30 -2843.9354 -1528.9634
31 -992.6728 -2843.9354
32 926.9779 -992.6728
33 732.3664 926.9779
34 680.3265 732.3664
35 -1131.3715 680.3265
36 -2022.7706 -1131.3715
37 -2717.0329 -2022.7706
38 -5113.3864 -2717.0329
39 -4699.8959 -5113.3864
40 -3971.7996 -4699.8959
41 -7358.5675 -3971.7996
42 -8975.0226 -7358.5675
43 -6978.6284 -8975.0226
44 -4483.8637 -6978.6284
45 -2263.4961 -4483.8637
46 -3826.1284 -2263.4961
47 -2845.1098 -3826.1284
48 -6493.3783 -2845.1098
49 -6764.2782 -6493.3783
50 -9276.5728 -6764.2782
51 -8428.2785 -9276.5728
52 -10292.6265 -8428.2785
53 -9717.0857 -10292.6265
54 -9885.0140 -9717.0857
55 -8089.1845 -9885.0140
56 -2784.1582 -8089.1845
57 966.2350 -2784.1582
58 3585.7736 966.2350
59 6632.9699 3585.7736
60 7969.2349 6632.9699
61 9431.1340 7969.2349
62 14041.3337 9431.1340
63 14242.5920 14041.3337
64 14638.6511 14242.5920
65 17456.5466 14638.6511
66 16970.9382 17456.5466
> 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/759mm1259606617.ps",horizontal=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/84a2q1259606617.ps",horizontal=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/9xzj81259606617.ps",horizontal=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/10thbi1259606617.ps",horizontal=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/1129pu1259606617.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/122q541259606617.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/139cbw1259606617.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/144kq31259606617.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/1518ee1259606617.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/16t1tb1259606617.tab")
+ }
> system("convert tmp/1kau81259606617.ps tmp/1kau81259606617.png")
> system("convert tmp/2sxbi1259606617.ps tmp/2sxbi1259606617.png")
> system("convert tmp/3yt691259606617.ps tmp/3yt691259606617.png")
> system("convert tmp/4bds71259606617.ps tmp/4bds71259606617.png")
> system("convert tmp/57a7w1259606617.ps tmp/57a7w1259606617.png")
> system("convert tmp/60elt1259606617.ps tmp/60elt1259606617.png")
> system("convert tmp/759mm1259606617.ps tmp/759mm1259606617.png")
> system("convert tmp/84a2q1259606617.ps tmp/84a2q1259606617.png")
> system("convert tmp/9xzj81259606617.ps tmp/9xzj81259606617.png")
> system("convert tmp/10thbi1259606617.ps tmp/10thbi1259606617.png")
>
>
> proc.time()
user system elapsed
2.479 1.596 4.773