R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(-1.2,23.6,-2.4,25.7,0.8,32.5,-0.1,33.5,-1.5,34.5,-4.4,27.9,-4.2,45.3,3.5,40.8,10,58.5,8.6,32.5,9.5,35.5,9.9,46.7,10.4,53.2,16,36.1,12.7,54,10.2,58.1,8.9,41.8,12.6,43.1,13.6,76,14.8,42.8,9.5,41,13.7,61.4,17,34.2,14.7,53.8,17.4,80.7,9,79.5,9.1,96.5,12.2,108.3,15.9,100.1,12.9,108.5,10.9,127.4,10.6,86.5,13.2,71.4,9.6,88.2,6.4,135.6,5.8,70.5,-1,87.5,-0.2,73.3,2.7,92.2,3.6,61.1,-0.9,45.7,0.3,30.5,-1.1,34.8,-2.5,29.2,-3.4,56.7,-3.5,67.1,-3.9,41.8,-4.6,46.8,-0.1,50.1,4.3,81.9,10.2,115.8,8.7,102.5,13.3,106.6,15,101.4,20.7,136.1,20.7,143.4,26.4,127.5,31.2,113.8,31.4,75.3,26.6,98.5,26.6,113.7,19.2,103.7,6.5,73.9,3.1,52.5,-0.2,63.9,-4,44.9,-12.6,31.3,-13,24.9,-17.6,22.8,-21.7,24.8,-23.2,22.8,-16.8,20.9,-19.8,21.5),dim=c(2,73),dimnames=list(c('werkl','afzetp'),1:73))
> y <- array(NA,dim=c(2,73),dimnames=list(c('werkl','afzetp'),1:73))
> 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
werkl afzetp M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 -1.2 23.6 1 0 0 0 0 0 0 0 0 0 0 1
2 -2.4 25.7 0 1 0 0 0 0 0 0 0 0 0 2
3 0.8 32.5 0 0 1 0 0 0 0 0 0 0 0 3
4 -0.1 33.5 0 0 0 1 0 0 0 0 0 0 0 4
5 -1.5 34.5 0 0 0 0 1 0 0 0 0 0 0 5
6 -4.4 27.9 0 0 0 0 0 1 0 0 0 0 0 6
7 -4.2 45.3 0 0 0 0 0 0 1 0 0 0 0 7
8 3.5 40.8 0 0 0 0 0 0 0 1 0 0 0 8
9 10.0 58.5 0 0 0 0 0 0 0 0 1 0 0 9
10 8.6 32.5 0 0 0 0 0 0 0 0 0 1 0 10
11 9.5 35.5 0 0 0 0 0 0 0 0 0 0 1 11
12 9.9 46.7 0 0 0 0 0 0 0 0 0 0 0 12
13 10.4 53.2 1 0 0 0 0 0 0 0 0 0 0 13
14 16.0 36.1 0 1 0 0 0 0 0 0 0 0 0 14
15 12.7 54.0 0 0 1 0 0 0 0 0 0 0 0 15
16 10.2 58.1 0 0 0 1 0 0 0 0 0 0 0 16
17 8.9 41.8 0 0 0 0 1 0 0 0 0 0 0 17
18 12.6 43.1 0 0 0 0 0 1 0 0 0 0 0 18
19 13.6 76.0 0 0 0 0 0 0 1 0 0 0 0 19
20 14.8 42.8 0 0 0 0 0 0 0 1 0 0 0 20
21 9.5 41.0 0 0 0 0 0 0 0 0 1 0 0 21
22 13.7 61.4 0 0 0 0 0 0 0 0 0 1 0 22
23 17.0 34.2 0 0 0 0 0 0 0 0 0 0 1 23
24 14.7 53.8 0 0 0 0 0 0 0 0 0 0 0 24
25 17.4 80.7 1 0 0 0 0 0 0 0 0 0 0 25
26 9.0 79.5 0 1 0 0 0 0 0 0 0 0 0 26
27 9.1 96.5 0 0 1 0 0 0 0 0 0 0 0 27
28 12.2 108.3 0 0 0 1 0 0 0 0 0 0 0 28
29 15.9 100.1 0 0 0 0 1 0 0 0 0 0 0 29
30 12.9 108.5 0 0 0 0 0 1 0 0 0 0 0 30
31 10.9 127.4 0 0 0 0 0 0 1 0 0 0 0 31
32 10.6 86.5 0 0 0 0 0 0 0 1 0 0 0 32
33 13.2 71.4 0 0 0 0 0 0 0 0 1 0 0 33
34 9.6 88.2 0 0 0 0 0 0 0 0 0 1 0 34
35 6.4 135.6 0 0 0 0 0 0 0 0 0 0 1 35
36 5.8 70.5 0 0 0 0 0 0 0 0 0 0 0 36
37 -1.0 87.5 1 0 0 0 0 0 0 0 0 0 0 37
38 -0.2 73.3 0 1 0 0 0 0 0 0 0 0 0 38
39 2.7 92.2 0 0 1 0 0 0 0 0 0 0 0 39
40 3.6 61.1 0 0 0 1 0 0 0 0 0 0 0 40
41 -0.9 45.7 0 0 0 0 1 0 0 0 0 0 0 41
42 0.3 30.5 0 0 0 0 0 1 0 0 0 0 0 42
43 -1.1 34.8 0 0 0 0 0 0 1 0 0 0 0 43
44 -2.5 29.2 0 0 0 0 0 0 0 1 0 0 0 44
45 -3.4 56.7 0 0 0 0 0 0 0 0 1 0 0 45
46 -3.5 67.1 0 0 0 0 0 0 0 0 0 1 0 46
47 -3.9 41.8 0 0 0 0 0 0 0 0 0 0 1 47
48 -4.6 46.8 0 0 0 0 0 0 0 0 0 0 0 48
49 -0.1 50.1 1 0 0 0 0 0 0 0 0 0 0 49
50 4.3 81.9 0 1 0 0 0 0 0 0 0 0 0 50
51 10.2 115.8 0 0 1 0 0 0 0 0 0 0 0 51
52 8.7 102.5 0 0 0 1 0 0 0 0 0 0 0 52
53 13.3 106.6 0 0 0 0 1 0 0 0 0 0 0 53
54 15.0 101.4 0 0 0 0 0 1 0 0 0 0 0 54
55 20.7 136.1 0 0 0 0 0 0 1 0 0 0 0 55
56 20.7 143.4 0 0 0 0 0 0 0 1 0 0 0 56
57 26.4 127.5 0 0 0 0 0 0 0 0 1 0 0 57
58 31.2 113.8 0 0 0 0 0 0 0 0 0 1 0 58
59 31.4 75.3 0 0 0 0 0 0 0 0 0 0 1 59
60 26.6 98.5 0 0 0 0 0 0 0 0 0 0 0 60
61 26.6 113.7 1 0 0 0 0 0 0 0 0 0 0 61
62 19.2 103.7 0 1 0 0 0 0 0 0 0 0 0 62
63 6.5 73.9 0 0 1 0 0 0 0 0 0 0 0 63
64 3.1 52.5 0 0 0 1 0 0 0 0 0 0 0 64
65 -0.2 63.9 0 0 0 0 1 0 0 0 0 0 0 65
66 -4.0 44.9 0 0 0 0 0 1 0 0 0 0 0 66
67 -12.6 31.3 0 0 0 0 0 0 1 0 0 0 0 67
68 -13.0 24.9 0 0 0 0 0 0 0 1 0 0 0 68
69 -17.6 22.8 0 0 0 0 0 0 0 0 1 0 0 69
70 -21.7 24.8 0 0 0 0 0 0 0 0 0 1 0 70
71 -23.2 22.8 0 0 0 0 0 0 0 0 0 0 1 71
72 -16.8 20.9 0 0 0 0 0 0 0 0 0 0 0 72
73 -19.8 21.5 1 0 0 0 0 0 0 0 0 0 0 73
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) afzetp M1 M2 M3 M4
-0.3957 0.2696 -3.7901 -3.2140 -6.5608 -4.8705
M5 M6 M7 M8 M9 M10
-3.9758 -2.6516 -7.5417 -2.4559 -2.0420 -2.3101
M11 t
-0.3028 -0.2100
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-22.1044 -5.2982 -0.2371 5.0406 24.1910
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.39571 4.14503 -0.095 0.924
afzetp 0.26957 0.03108 8.674 4.01e-12 ***
M1 -3.79014 4.67620 -0.811 0.421
M2 -3.21396 4.88417 -0.658 0.513
M3 -6.56078 4.91988 -1.334 0.187
M4 -4.87045 4.88165 -0.998 0.322
M5 -3.97579 4.86632 -0.817 0.417
M6 -2.65155 4.85224 -0.546 0.587
M7 -7.54173 4.88991 -1.542 0.128
M8 -2.45587 4.84902 -0.506 0.614
M9 -2.04195 4.84936 -0.421 0.675
M10 -2.31005 4.85038 -0.476 0.636
M11 -0.30277 4.84156 -0.063 0.950
t -0.21002 0.04853 -4.328 5.91e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.385 on 59 degrees of freedom
Multiple R-squared: 0.5791, Adjusted R-squared: 0.4863
F-statistic: 6.244 on 13 and 59 DF, p-value: 3.467e-07
> 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.026955489 0.053910979 0.9730445
[2,] 0.010480089 0.020960178 0.9895199
[3,] 0.006653067 0.013306134 0.9933469
[4,] 0.002294576 0.004589153 0.9977054
[5,] 0.003466280 0.006932559 0.9965337
[6,] 0.006329758 0.012659515 0.9936702
[7,] 0.003781685 0.007563369 0.9962183
[8,] 0.002406009 0.004812019 0.9975940
[9,] 0.002252572 0.004505143 0.9977474
[10,] 0.012558721 0.025117441 0.9874413
[11,] 0.014465818 0.028931635 0.9855342
[12,] 0.007764797 0.015529594 0.9922352
[13,] 0.004143905 0.008287809 0.9958561
[14,] 0.002074347 0.004148694 0.9979257
[15,] 0.001298579 0.002597157 0.9987014
[16,] 0.001338812 0.002677625 0.9986612
[17,] 0.001334097 0.002668193 0.9986659
[18,] 0.001737135 0.003474270 0.9982629
[19,] 0.009387243 0.018774487 0.9906128
[20,] 0.022998569 0.045997138 0.9770014
[21,] 0.099002523 0.198005046 0.9009975
[22,] 0.128507798 0.257015596 0.8714922
[23,] 0.122944689 0.245889379 0.8770553
[24,] 0.091594032 0.183188063 0.9084060
[25,] 0.074262790 0.148525579 0.9257372
[26,] 0.056614480 0.113228960 0.9433855
[27,] 0.050645734 0.101291469 0.9493543
[28,] 0.073358037 0.146716074 0.9266420
[29,] 0.065832846 0.131665691 0.9341672
[30,] 0.059480183 0.118960367 0.9405198
[31,] 0.042663187 0.085326374 0.9573368
[32,] 0.031384577 0.062769153 0.9686154
[33,] 0.018884670 0.037769340 0.9811153
[34,] 0.011233389 0.022466778 0.9887666
[35,] 0.013536887 0.027073774 0.9864631
[36,] 0.035793261 0.071586522 0.9642067
[37,] 0.105008939 0.210017878 0.8949911
[38,] 0.857667767 0.284664466 0.1423322
[39,] 0.873848261 0.252303479 0.1261517
[40,] 0.854302040 0.291395920 0.1456980
> postscript(file="/var/www/html/rcomp/tmp/11ieu1261311809.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/2j24p1261311809.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/3ndtj1261311809.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/4a1pb1261311809.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/51xzq1261311809.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 = 73
Frequency = 1
1 2 3 4 5 6
-3.16593859 -5.29819440 -0.37441545 -3.02429539 -5.37851161 -7.61358258
7 8 9 10 11 12
-7.00387687 -2.96665599 -1.44191682 4.64496992 2.93900044 0.22708232
13 14 15 16 17 18
2.97504242 12.81849097 8.25006584 3.16452530 5.57383436 7.80917670
19 20 21 22 23 24
5.04057942 10.31440003 5.29571393 4.47464849 13.30963065 5.63334188
25 26 27 28 29 30
5.08211610 -3.36056551 -4.28637951 -5.84759314 -0.62178435 -7.00037435
31 32 33 34 35 36
-8.99502054 -3.14552684 3.32104060 -4.32958028 -22.10436606 -5.24825075
37 38 39 40 41 42
-12.63075396 -8.36905243 -7.00704550 0.79620528 -0.23709681 3.94611646
43 44 45 46 47 48
6.48716212 1.72090772 -6.79611886 -9.22150496 -4.59870188 -6.73929880
49 50 51 52 53 54
0.87127869 -3.66714477 -3.34865687 -2.74371534 0.06640785 2.05394177
55 56 57 58 59 60
3.50012220 -3.34355855 6.43866324 15.40986437 24.19096419 13.04423086
61 62 63 64 65 66
12.94694992 7.87646613 6.76643149 7.65487330 0.59715056 0.80472198
67 68 69 70 71 72
0.97103367 -2.57956638 -6.81738209 -10.97839753 -13.73652734 -6.91710551
73
-6.07869459
> postscript(file="/var/www/html/rcomp/tmp/66q3i1261311809.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 -3.16593859 NA
1 -5.29819440 -3.16593859
2 -0.37441545 -5.29819440
3 -3.02429539 -0.37441545
4 -5.37851161 -3.02429539
5 -7.61358258 -5.37851161
6 -7.00387687 -7.61358258
7 -2.96665599 -7.00387687
8 -1.44191682 -2.96665599
9 4.64496992 -1.44191682
10 2.93900044 4.64496992
11 0.22708232 2.93900044
12 2.97504242 0.22708232
13 12.81849097 2.97504242
14 8.25006584 12.81849097
15 3.16452530 8.25006584
16 5.57383436 3.16452530
17 7.80917670 5.57383436
18 5.04057942 7.80917670
19 10.31440003 5.04057942
20 5.29571393 10.31440003
21 4.47464849 5.29571393
22 13.30963065 4.47464849
23 5.63334188 13.30963065
24 5.08211610 5.63334188
25 -3.36056551 5.08211610
26 -4.28637951 -3.36056551
27 -5.84759314 -4.28637951
28 -0.62178435 -5.84759314
29 -7.00037435 -0.62178435
30 -8.99502054 -7.00037435
31 -3.14552684 -8.99502054
32 3.32104060 -3.14552684
33 -4.32958028 3.32104060
34 -22.10436606 -4.32958028
35 -5.24825075 -22.10436606
36 -12.63075396 -5.24825075
37 -8.36905243 -12.63075396
38 -7.00704550 -8.36905243
39 0.79620528 -7.00704550
40 -0.23709681 0.79620528
41 3.94611646 -0.23709681
42 6.48716212 3.94611646
43 1.72090772 6.48716212
44 -6.79611886 1.72090772
45 -9.22150496 -6.79611886
46 -4.59870188 -9.22150496
47 -6.73929880 -4.59870188
48 0.87127869 -6.73929880
49 -3.66714477 0.87127869
50 -3.34865687 -3.66714477
51 -2.74371534 -3.34865687
52 0.06640785 -2.74371534
53 2.05394177 0.06640785
54 3.50012220 2.05394177
55 -3.34355855 3.50012220
56 6.43866324 -3.34355855
57 15.40986437 6.43866324
58 24.19096419 15.40986437
59 13.04423086 24.19096419
60 12.94694992 13.04423086
61 7.87646613 12.94694992
62 6.76643149 7.87646613
63 7.65487330 6.76643149
64 0.59715056 7.65487330
65 0.80472198 0.59715056
66 0.97103367 0.80472198
67 -2.57956638 0.97103367
68 -6.81738209 -2.57956638
69 -10.97839753 -6.81738209
70 -13.73652734 -10.97839753
71 -6.91710551 -13.73652734
72 -6.07869459 -6.91710551
73 NA -6.07869459
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -5.29819440 -3.16593859
[2,] -0.37441545 -5.29819440
[3,] -3.02429539 -0.37441545
[4,] -5.37851161 -3.02429539
[5,] -7.61358258 -5.37851161
[6,] -7.00387687 -7.61358258
[7,] -2.96665599 -7.00387687
[8,] -1.44191682 -2.96665599
[9,] 4.64496992 -1.44191682
[10,] 2.93900044 4.64496992
[11,] 0.22708232 2.93900044
[12,] 2.97504242 0.22708232
[13,] 12.81849097 2.97504242
[14,] 8.25006584 12.81849097
[15,] 3.16452530 8.25006584
[16,] 5.57383436 3.16452530
[17,] 7.80917670 5.57383436
[18,] 5.04057942 7.80917670
[19,] 10.31440003 5.04057942
[20,] 5.29571393 10.31440003
[21,] 4.47464849 5.29571393
[22,] 13.30963065 4.47464849
[23,] 5.63334188 13.30963065
[24,] 5.08211610 5.63334188
[25,] -3.36056551 5.08211610
[26,] -4.28637951 -3.36056551
[27,] -5.84759314 -4.28637951
[28,] -0.62178435 -5.84759314
[29,] -7.00037435 -0.62178435
[30,] -8.99502054 -7.00037435
[31,] -3.14552684 -8.99502054
[32,] 3.32104060 -3.14552684
[33,] -4.32958028 3.32104060
[34,] -22.10436606 -4.32958028
[35,] -5.24825075 -22.10436606
[36,] -12.63075396 -5.24825075
[37,] -8.36905243 -12.63075396
[38,] -7.00704550 -8.36905243
[39,] 0.79620528 -7.00704550
[40,] -0.23709681 0.79620528
[41,] 3.94611646 -0.23709681
[42,] 6.48716212 3.94611646
[43,] 1.72090772 6.48716212
[44,] -6.79611886 1.72090772
[45,] -9.22150496 -6.79611886
[46,] -4.59870188 -9.22150496
[47,] -6.73929880 -4.59870188
[48,] 0.87127869 -6.73929880
[49,] -3.66714477 0.87127869
[50,] -3.34865687 -3.66714477
[51,] -2.74371534 -3.34865687
[52,] 0.06640785 -2.74371534
[53,] 2.05394177 0.06640785
[54,] 3.50012220 2.05394177
[55,] -3.34355855 3.50012220
[56,] 6.43866324 -3.34355855
[57,] 15.40986437 6.43866324
[58,] 24.19096419 15.40986437
[59,] 13.04423086 24.19096419
[60,] 12.94694992 13.04423086
[61,] 7.87646613 12.94694992
[62,] 6.76643149 7.87646613
[63,] 7.65487330 6.76643149
[64,] 0.59715056 7.65487330
[65,] 0.80472198 0.59715056
[66,] 0.97103367 0.80472198
[67,] -2.57956638 0.97103367
[68,] -6.81738209 -2.57956638
[69,] -10.97839753 -6.81738209
[70,] -13.73652734 -10.97839753
[71,] -6.91710551 -13.73652734
[72,] -6.07869459 -6.91710551
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -5.29819440 -3.16593859
2 -0.37441545 -5.29819440
3 -3.02429539 -0.37441545
4 -5.37851161 -3.02429539
5 -7.61358258 -5.37851161
6 -7.00387687 -7.61358258
7 -2.96665599 -7.00387687
8 -1.44191682 -2.96665599
9 4.64496992 -1.44191682
10 2.93900044 4.64496992
11 0.22708232 2.93900044
12 2.97504242 0.22708232
13 12.81849097 2.97504242
14 8.25006584 12.81849097
15 3.16452530 8.25006584
16 5.57383436 3.16452530
17 7.80917670 5.57383436
18 5.04057942 7.80917670
19 10.31440003 5.04057942
20 5.29571393 10.31440003
21 4.47464849 5.29571393
22 13.30963065 4.47464849
23 5.63334188 13.30963065
24 5.08211610 5.63334188
25 -3.36056551 5.08211610
26 -4.28637951 -3.36056551
27 -5.84759314 -4.28637951
28 -0.62178435 -5.84759314
29 -7.00037435 -0.62178435
30 -8.99502054 -7.00037435
31 -3.14552684 -8.99502054
32 3.32104060 -3.14552684
33 -4.32958028 3.32104060
34 -22.10436606 -4.32958028
35 -5.24825075 -22.10436606
36 -12.63075396 -5.24825075
37 -8.36905243 -12.63075396
38 -7.00704550 -8.36905243
39 0.79620528 -7.00704550
40 -0.23709681 0.79620528
41 3.94611646 -0.23709681
42 6.48716212 3.94611646
43 1.72090772 6.48716212
44 -6.79611886 1.72090772
45 -9.22150496 -6.79611886
46 -4.59870188 -9.22150496
47 -6.73929880 -4.59870188
48 0.87127869 -6.73929880
49 -3.66714477 0.87127869
50 -3.34865687 -3.66714477
51 -2.74371534 -3.34865687
52 0.06640785 -2.74371534
53 2.05394177 0.06640785
54 3.50012220 2.05394177
55 -3.34355855 3.50012220
56 6.43866324 -3.34355855
57 15.40986437 6.43866324
58 24.19096419 15.40986437
59 13.04423086 24.19096419
60 12.94694992 13.04423086
61 7.87646613 12.94694992
62 6.76643149 7.87646613
63 7.65487330 6.76643149
64 0.59715056 7.65487330
65 0.80472198 0.59715056
66 0.97103367 0.80472198
67 -2.57956638 0.97103367
68 -6.81738209 -2.57956638
69 -10.97839753 -6.81738209
70 -13.73652734 -10.97839753
71 -6.91710551 -13.73652734
72 -6.07869459 -6.91710551
> 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/70r8f1261311809.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/8vy5e1261311809.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/9n8px1261311809.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/10ml4a1261311809.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/110cbk1261311809.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/12ktcg1261311809.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/13kihv1261311809.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/143qff1261311809.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/151jlk1261311809.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/16vcao1261311809.tab")
+ }
>
> try(system("convert tmp/11ieu1261311809.ps tmp/11ieu1261311809.png",intern=TRUE))
character(0)
> try(system("convert tmp/2j24p1261311809.ps tmp/2j24p1261311809.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ndtj1261311809.ps tmp/3ndtj1261311809.png",intern=TRUE))
character(0)
> try(system("convert tmp/4a1pb1261311809.ps tmp/4a1pb1261311809.png",intern=TRUE))
character(0)
> try(system("convert tmp/51xzq1261311809.ps tmp/51xzq1261311809.png",intern=TRUE))
character(0)
> try(system("convert tmp/66q3i1261311809.ps tmp/66q3i1261311809.png",intern=TRUE))
character(0)
> try(system("convert tmp/70r8f1261311809.ps tmp/70r8f1261311809.png",intern=TRUE))
character(0)
> try(system("convert tmp/8vy5e1261311809.ps tmp/8vy5e1261311809.png",intern=TRUE))
character(0)
> try(system("convert tmp/9n8px1261311809.ps tmp/9n8px1261311809.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ml4a1261311809.ps tmp/10ml4a1261311809.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.576 1.558 3.154