R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(93.7,105.7,109.5,105.3,102.8,100.6,97.6,110.3,107.2,107.2,108.1,97.1,92.2,112.2,111.6,115.7,111.3,104.2,103.2,112.7,106.4,102.6,110.6,95.2,89,112.5,116.8,107.2,113.6,101.8,102.6,122.7,110.3,110.5,121.6,100.3,100.7,123.4,127.1,124.1,131.2,111.6,114.2,130.1,125.9,119,133.8,107.5,113.5,134.4,126.8,135.6,139.9,129.8,131,153.1,134.1,144.1,155.9,123.3,128.1,144.3,153,149.9,150.9,141,138.9,157.4,142.9,151.7,161,138.5,135.9,151.5,164,159.1,157,142.1,144.8,152.1,154.6,148.7,157.7,146.4,136.5),dim=c(1,85),dimnames=list(c('omzet'),1:85))
> y <- array(NA,dim=c(1,85),dimnames=list(c('omzet'),1:85))
> 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
omzet M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 93.7 1 0 0 0 0 0 0 0 0 0 0 1
2 105.7 0 1 0 0 0 0 0 0 0 0 0 2
3 109.5 0 0 1 0 0 0 0 0 0 0 0 3
4 105.3 0 0 0 1 0 0 0 0 0 0 0 4
5 102.8 0 0 0 0 1 0 0 0 0 0 0 5
6 100.6 0 0 0 0 0 1 0 0 0 0 0 6
7 97.6 0 0 0 0 0 0 1 0 0 0 0 7
8 110.3 0 0 0 0 0 0 0 1 0 0 0 8
9 107.2 0 0 0 0 0 0 0 0 1 0 0 9
10 107.2 0 0 0 0 0 0 0 0 0 1 0 10
11 108.1 0 0 0 0 0 0 0 0 0 0 1 11
12 97.1 0 0 0 0 0 0 0 0 0 0 0 12
13 92.2 1 0 0 0 0 0 0 0 0 0 0 13
14 112.2 0 1 0 0 0 0 0 0 0 0 0 14
15 111.6 0 0 1 0 0 0 0 0 0 0 0 15
16 115.7 0 0 0 1 0 0 0 0 0 0 0 16
17 111.3 0 0 0 0 1 0 0 0 0 0 0 17
18 104.2 0 0 0 0 0 1 0 0 0 0 0 18
19 103.2 0 0 0 0 0 0 1 0 0 0 0 19
20 112.7 0 0 0 0 0 0 0 1 0 0 0 20
21 106.4 0 0 0 0 0 0 0 0 1 0 0 21
22 102.6 0 0 0 0 0 0 0 0 0 1 0 22
23 110.6 0 0 0 0 0 0 0 0 0 0 1 23
24 95.2 0 0 0 0 0 0 0 0 0 0 0 24
25 89.0 1 0 0 0 0 0 0 0 0 0 0 25
26 112.5 0 1 0 0 0 0 0 0 0 0 0 26
27 116.8 0 0 1 0 0 0 0 0 0 0 0 27
28 107.2 0 0 0 1 0 0 0 0 0 0 0 28
29 113.6 0 0 0 0 1 0 0 0 0 0 0 29
30 101.8 0 0 0 0 0 1 0 0 0 0 0 30
31 102.6 0 0 0 0 0 0 1 0 0 0 0 31
32 122.7 0 0 0 0 0 0 0 1 0 0 0 32
33 110.3 0 0 0 0 0 0 0 0 1 0 0 33
34 110.5 0 0 0 0 0 0 0 0 0 1 0 34
35 121.6 0 0 0 0 0 0 0 0 0 0 1 35
36 100.3 0 0 0 0 0 0 0 0 0 0 0 36
37 100.7 1 0 0 0 0 0 0 0 0 0 0 37
38 123.4 0 1 0 0 0 0 0 0 0 0 0 38
39 127.1 0 0 1 0 0 0 0 0 0 0 0 39
40 124.1 0 0 0 1 0 0 0 0 0 0 0 40
41 131.2 0 0 0 0 1 0 0 0 0 0 0 41
42 111.6 0 0 0 0 0 1 0 0 0 0 0 42
43 114.2 0 0 0 0 0 0 1 0 0 0 0 43
44 130.1 0 0 0 0 0 0 0 1 0 0 0 44
45 125.9 0 0 0 0 0 0 0 0 1 0 0 45
46 119.0 0 0 0 0 0 0 0 0 0 1 0 46
47 133.8 0 0 0 0 0 0 0 0 0 0 1 47
48 107.5 0 0 0 0 0 0 0 0 0 0 0 48
49 113.5 1 0 0 0 0 0 0 0 0 0 0 49
50 134.4 0 1 0 0 0 0 0 0 0 0 0 50
51 126.8 0 0 1 0 0 0 0 0 0 0 0 51
52 135.6 0 0 0 1 0 0 0 0 0 0 0 52
53 139.9 0 0 0 0 1 0 0 0 0 0 0 53
54 129.8 0 0 0 0 0 1 0 0 0 0 0 54
55 131.0 0 0 0 0 0 0 1 0 0 0 0 55
56 153.1 0 0 0 0 0 0 0 1 0 0 0 56
57 134.1 0 0 0 0 0 0 0 0 1 0 0 57
58 144.1 0 0 0 0 0 0 0 0 0 1 0 58
59 155.9 0 0 0 0 0 0 0 0 0 0 1 59
60 123.3 0 0 0 0 0 0 0 0 0 0 0 60
61 128.1 1 0 0 0 0 0 0 0 0 0 0 61
62 144.3 0 1 0 0 0 0 0 0 0 0 0 62
63 153.0 0 0 1 0 0 0 0 0 0 0 0 63
64 149.9 0 0 0 1 0 0 0 0 0 0 0 64
65 150.9 0 0 0 0 1 0 0 0 0 0 0 65
66 141.0 0 0 0 0 0 1 0 0 0 0 0 66
67 138.9 0 0 0 0 0 0 1 0 0 0 0 67
68 157.4 0 0 0 0 0 0 0 1 0 0 0 68
69 142.9 0 0 0 0 0 0 0 0 1 0 0 69
70 151.7 0 0 0 0 0 0 0 0 0 1 0 70
71 161.0 0 0 0 0 0 0 0 0 0 0 1 71
72 138.5 0 0 0 0 0 0 0 0 0 0 0 72
73 135.9 1 0 0 0 0 0 0 0 0 0 0 73
74 151.5 0 1 0 0 0 0 0 0 0 0 0 74
75 164.0 0 0 1 0 0 0 0 0 0 0 0 75
76 159.1 0 0 0 1 0 0 0 0 0 0 0 76
77 157.0 0 0 0 0 1 0 0 0 0 0 0 77
78 142.1 0 0 0 0 0 1 0 0 0 0 0 78
79 144.8 0 0 0 0 0 0 1 0 0 0 0 79
80 152.1 0 0 0 0 0 0 0 1 0 0 0 80
81 154.6 0 0 0 0 0 0 0 0 1 0 0 81
82 148.7 0 0 0 0 0 0 0 0 0 1 0 82
83 157.7 0 0 0 0 0 0 0 0 0 0 1 83
84 146.4 0 0 0 0 0 0 0 0 0 0 0 84
85 136.5 1 0 0 0 0 0 0 0 0 0 0 85
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
80.2760 -0.6052 18.1467 20.9563 18.5230 19.1898
M6 M7 M8 M9 M10 M11
7.6566 7.0948 21.5187 12.6426 12.2522 20.7904
t
0.7332
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12.130 -4.225 -0.612 4.174 13.296
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 80.27600 2.63295 30.489 < 2e-16 ***
M1 -0.60524 3.14092 -0.193 0.847741
M2 18.14667 3.25217 5.580 3.99e-07 ***
M3 20.95629 3.25003 6.448 1.13e-08 ***
M4 18.52305 3.24813 5.703 2.43e-07 ***
M5 19.18981 3.24644 5.911 1.04e-07 ***
M6 7.65657 3.24498 2.360 0.021014 *
M7 7.09476 3.24374 2.187 0.031976 *
M8 21.51867 3.24273 6.636 5.15e-09 ***
M9 12.64257 3.24194 3.900 0.000214 ***
M10 12.25219 3.24138 3.780 0.000321 ***
M11 20.79038 3.24104 6.415 1.30e-08 ***
t 0.73324 0.02701 27.149 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.063 on 72 degrees of freedom
Multiple R-squared: 0.9226, Adjusted R-squared: 0.9097
F-statistic: 71.5 on 12 and 72 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.31996865 0.63993730 0.68003135
[2,] 0.22016306 0.44032613 0.77983694
[3,] 0.13652771 0.27305542 0.86347229
[4,] 0.07784750 0.15569501 0.92215250
[5,] 0.04461690 0.08923380 0.95538310
[6,] 0.04439489 0.08878978 0.95560511
[7,] 0.08308840 0.16617680 0.91691160
[8,] 0.04869398 0.09738797 0.95130602
[9,] 0.04183173 0.08366346 0.95816827
[10,] 0.06518716 0.13037432 0.93481284
[11,] 0.03961686 0.07923373 0.96038314
[12,] 0.02824961 0.05649922 0.97175039
[13,] 0.03812594 0.07625188 0.96187406
[14,] 0.03147030 0.06294060 0.96852970
[15,] 0.02297035 0.04594069 0.97702965
[16,] 0.01392222 0.02784444 0.98607778
[17,] 0.02544738 0.05089476 0.97455262
[18,] 0.01581731 0.03163462 0.98418269
[19,] 0.01093675 0.02187349 0.98906325
[20,] 0.02077131 0.04154261 0.97922869
[21,] 0.01323520 0.02647041 0.98676480
[22,] 0.00990145 0.01980290 0.99009855
[23,] 0.01253281 0.02506563 0.98746719
[24,] 0.01623512 0.03247023 0.98376488
[25,] 0.01980218 0.03960435 0.98019782
[26,] 0.06448907 0.12897814 0.93551093
[27,] 0.05198527 0.10397053 0.94801473
[28,] 0.04346398 0.08692796 0.95653602
[29,] 0.04084234 0.08168468 0.95915766
[30,] 0.04295775 0.08591550 0.95704225
[31,] 0.05115438 0.10230876 0.94884562
[32,] 0.07836738 0.15673475 0.92163262
[33,] 0.10973100 0.21946199 0.89026900
[34,] 0.10549589 0.21099178 0.89450411
[35,] 0.09744103 0.19488205 0.90255897
[36,] 0.46600625 0.93201251 0.53399375
[37,] 0.60549188 0.78901625 0.39450812
[38,] 0.66675661 0.66648678 0.33324339
[39,] 0.68121678 0.63756644 0.31878322
[40,] 0.69154690 0.61690620 0.30845310
[41,] 0.82250212 0.35499577 0.17749788
[42,] 0.82420956 0.35158088 0.17579044
[43,] 0.84854403 0.30291193 0.15145597
[44,] 0.89129855 0.21740290 0.10870145
[45,] 0.95636740 0.08726521 0.04363260
[46,] 0.93385200 0.13229601 0.06614800
[47,] 0.90153868 0.19692264 0.09846132
[48,] 0.91332627 0.17334746 0.08667373
[49,] 0.90639641 0.18720719 0.09360359
[50,] 0.86613298 0.26773404 0.13386702
[51,] 0.78732252 0.42535495 0.21267748
[52,] 0.70443420 0.59113161 0.29556580
[53,] 0.68307407 0.63385186 0.31692593
[54,] 0.78896589 0.42206822 0.21103411
> postscript(file="/var/www/html/freestat/rcomp/tmp/1u55k1227794852.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/freestat/rcomp/tmp/26dez1227794852.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/freestat/rcomp/tmp/32cxn1227794852.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/freestat/rcomp/tmp/4b1v21227794852.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/freestat/rcomp/tmp/5mgzu1227794852.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 = 85
Frequency = 1
1 2 3 4 5 6
13.29600000 5.81085714 6.06800000 3.56800000 -0.33200000 8.26800000
7 8 9 10 11 12
5.09657143 2.63942857 7.68228571 7.33942857 -1.03200000 8.02514286
13 14 15 16 17 18
2.99714286 3.51200000 -0.63085714 5.16914286 -0.63085714 3.06914286
19 20 21 22 23 24
1.89771429 -3.75942857 -1.91657143 -6.05942857 -7.33085714 -2.67371429
25 26 27 28 29 30
-9.00171429 -4.98685714 -4.22971429 -12.12971429 -7.12971429 -8.12971429
31 32 33 34 35 36
-7.50114286 -2.55828571 -6.81542857 -6.95828571 -5.12971429 -6.37257143
37 38 39 40 41 42
-6.10057143 -2.88571429 -2.72857143 -4.02857143 1.67142857 -7.12857143
43 44 45 46 47 48
-4.70000000 -3.95714286 -0.01428571 -7.25714286 -1.72857143 -7.97142857
49 50 51 52 53 54
-2.09942857 -0.68457143 -11.82742857 -1.32742857 1.57257143 2.27257143
55 56 57 58 59 60
3.30114286 10.24400000 -0.61314286 9.04400000 11.57257143 -0.97028571
61 62 63 64 65 66
3.70171429 0.41657143 5.57371429 4.17371429 3.77371429 4.67371429
67 68 69 70 71 72
2.40228571 5.74514286 -0.61200000 7.84514286 7.87371429 5.43085714
73 74 75 76 77 78
2.70285714 -1.18228571 7.77485714 4.57485714 1.07485714 -3.02514286
79 80 81 82 83 84
-0.49657143 -8.35371429 2.28914286 -3.95371429 -4.22514286 4.53200000
85
-5.49600000
> postscript(file="/var/www/html/freestat/rcomp/tmp/61xdi1227794852.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 = 85
Frequency = 1
lag(myerror, k = 1) myerror
0 13.29600000 NA
1 5.81085714 13.29600000
2 6.06800000 5.81085714
3 3.56800000 6.06800000
4 -0.33200000 3.56800000
5 8.26800000 -0.33200000
6 5.09657143 8.26800000
7 2.63942857 5.09657143
8 7.68228571 2.63942857
9 7.33942857 7.68228571
10 -1.03200000 7.33942857
11 8.02514286 -1.03200000
12 2.99714286 8.02514286
13 3.51200000 2.99714286
14 -0.63085714 3.51200000
15 5.16914286 -0.63085714
16 -0.63085714 5.16914286
17 3.06914286 -0.63085714
18 1.89771429 3.06914286
19 -3.75942857 1.89771429
20 -1.91657143 -3.75942857
21 -6.05942857 -1.91657143
22 -7.33085714 -6.05942857
23 -2.67371429 -7.33085714
24 -9.00171429 -2.67371429
25 -4.98685714 -9.00171429
26 -4.22971429 -4.98685714
27 -12.12971429 -4.22971429
28 -7.12971429 -12.12971429
29 -8.12971429 -7.12971429
30 -7.50114286 -8.12971429
31 -2.55828571 -7.50114286
32 -6.81542857 -2.55828571
33 -6.95828571 -6.81542857
34 -5.12971429 -6.95828571
35 -6.37257143 -5.12971429
36 -6.10057143 -6.37257143
37 -2.88571429 -6.10057143
38 -2.72857143 -2.88571429
39 -4.02857143 -2.72857143
40 1.67142857 -4.02857143
41 -7.12857143 1.67142857
42 -4.70000000 -7.12857143
43 -3.95714286 -4.70000000
44 -0.01428571 -3.95714286
45 -7.25714286 -0.01428571
46 -1.72857143 -7.25714286
47 -7.97142857 -1.72857143
48 -2.09942857 -7.97142857
49 -0.68457143 -2.09942857
50 -11.82742857 -0.68457143
51 -1.32742857 -11.82742857
52 1.57257143 -1.32742857
53 2.27257143 1.57257143
54 3.30114286 2.27257143
55 10.24400000 3.30114286
56 -0.61314286 10.24400000
57 9.04400000 -0.61314286
58 11.57257143 9.04400000
59 -0.97028571 11.57257143
60 3.70171429 -0.97028571
61 0.41657143 3.70171429
62 5.57371429 0.41657143
63 4.17371429 5.57371429
64 3.77371429 4.17371429
65 4.67371429 3.77371429
66 2.40228571 4.67371429
67 5.74514286 2.40228571
68 -0.61200000 5.74514286
69 7.84514286 -0.61200000
70 7.87371429 7.84514286
71 5.43085714 7.87371429
72 2.70285714 5.43085714
73 -1.18228571 2.70285714
74 7.77485714 -1.18228571
75 4.57485714 7.77485714
76 1.07485714 4.57485714
77 -3.02514286 1.07485714
78 -0.49657143 -3.02514286
79 -8.35371429 -0.49657143
80 2.28914286 -8.35371429
81 -3.95371429 2.28914286
82 -4.22514286 -3.95371429
83 4.53200000 -4.22514286
84 -5.49600000 4.53200000
85 NA -5.49600000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 5.81085714 13.29600000
[2,] 6.06800000 5.81085714
[3,] 3.56800000 6.06800000
[4,] -0.33200000 3.56800000
[5,] 8.26800000 -0.33200000
[6,] 5.09657143 8.26800000
[7,] 2.63942857 5.09657143
[8,] 7.68228571 2.63942857
[9,] 7.33942857 7.68228571
[10,] -1.03200000 7.33942857
[11,] 8.02514286 -1.03200000
[12,] 2.99714286 8.02514286
[13,] 3.51200000 2.99714286
[14,] -0.63085714 3.51200000
[15,] 5.16914286 -0.63085714
[16,] -0.63085714 5.16914286
[17,] 3.06914286 -0.63085714
[18,] 1.89771429 3.06914286
[19,] -3.75942857 1.89771429
[20,] -1.91657143 -3.75942857
[21,] -6.05942857 -1.91657143
[22,] -7.33085714 -6.05942857
[23,] -2.67371429 -7.33085714
[24,] -9.00171429 -2.67371429
[25,] -4.98685714 -9.00171429
[26,] -4.22971429 -4.98685714
[27,] -12.12971429 -4.22971429
[28,] -7.12971429 -12.12971429
[29,] -8.12971429 -7.12971429
[30,] -7.50114286 -8.12971429
[31,] -2.55828571 -7.50114286
[32,] -6.81542857 -2.55828571
[33,] -6.95828571 -6.81542857
[34,] -5.12971429 -6.95828571
[35,] -6.37257143 -5.12971429
[36,] -6.10057143 -6.37257143
[37,] -2.88571429 -6.10057143
[38,] -2.72857143 -2.88571429
[39,] -4.02857143 -2.72857143
[40,] 1.67142857 -4.02857143
[41,] -7.12857143 1.67142857
[42,] -4.70000000 -7.12857143
[43,] -3.95714286 -4.70000000
[44,] -0.01428571 -3.95714286
[45,] -7.25714286 -0.01428571
[46,] -1.72857143 -7.25714286
[47,] -7.97142857 -1.72857143
[48,] -2.09942857 -7.97142857
[49,] -0.68457143 -2.09942857
[50,] -11.82742857 -0.68457143
[51,] -1.32742857 -11.82742857
[52,] 1.57257143 -1.32742857
[53,] 2.27257143 1.57257143
[54,] 3.30114286 2.27257143
[55,] 10.24400000 3.30114286
[56,] -0.61314286 10.24400000
[57,] 9.04400000 -0.61314286
[58,] 11.57257143 9.04400000
[59,] -0.97028571 11.57257143
[60,] 3.70171429 -0.97028571
[61,] 0.41657143 3.70171429
[62,] 5.57371429 0.41657143
[63,] 4.17371429 5.57371429
[64,] 3.77371429 4.17371429
[65,] 4.67371429 3.77371429
[66,] 2.40228571 4.67371429
[67,] 5.74514286 2.40228571
[68,] -0.61200000 5.74514286
[69,] 7.84514286 -0.61200000
[70,] 7.87371429 7.84514286
[71,] 5.43085714 7.87371429
[72,] 2.70285714 5.43085714
[73,] -1.18228571 2.70285714
[74,] 7.77485714 -1.18228571
[75,] 4.57485714 7.77485714
[76,] 1.07485714 4.57485714
[77,] -3.02514286 1.07485714
[78,] -0.49657143 -3.02514286
[79,] -8.35371429 -0.49657143
[80,] 2.28914286 -8.35371429
[81,] -3.95371429 2.28914286
[82,] -4.22514286 -3.95371429
[83,] 4.53200000 -4.22514286
[84,] -5.49600000 4.53200000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 5.81085714 13.29600000
2 6.06800000 5.81085714
3 3.56800000 6.06800000
4 -0.33200000 3.56800000
5 8.26800000 -0.33200000
6 5.09657143 8.26800000
7 2.63942857 5.09657143
8 7.68228571 2.63942857
9 7.33942857 7.68228571
10 -1.03200000 7.33942857
11 8.02514286 -1.03200000
12 2.99714286 8.02514286
13 3.51200000 2.99714286
14 -0.63085714 3.51200000
15 5.16914286 -0.63085714
16 -0.63085714 5.16914286
17 3.06914286 -0.63085714
18 1.89771429 3.06914286
19 -3.75942857 1.89771429
20 -1.91657143 -3.75942857
21 -6.05942857 -1.91657143
22 -7.33085714 -6.05942857
23 -2.67371429 -7.33085714
24 -9.00171429 -2.67371429
25 -4.98685714 -9.00171429
26 -4.22971429 -4.98685714
27 -12.12971429 -4.22971429
28 -7.12971429 -12.12971429
29 -8.12971429 -7.12971429
30 -7.50114286 -8.12971429
31 -2.55828571 -7.50114286
32 -6.81542857 -2.55828571
33 -6.95828571 -6.81542857
34 -5.12971429 -6.95828571
35 -6.37257143 -5.12971429
36 -6.10057143 -6.37257143
37 -2.88571429 -6.10057143
38 -2.72857143 -2.88571429
39 -4.02857143 -2.72857143
40 1.67142857 -4.02857143
41 -7.12857143 1.67142857
42 -4.70000000 -7.12857143
43 -3.95714286 -4.70000000
44 -0.01428571 -3.95714286
45 -7.25714286 -0.01428571
46 -1.72857143 -7.25714286
47 -7.97142857 -1.72857143
48 -2.09942857 -7.97142857
49 -0.68457143 -2.09942857
50 -11.82742857 -0.68457143
51 -1.32742857 -11.82742857
52 1.57257143 -1.32742857
53 2.27257143 1.57257143
54 3.30114286 2.27257143
55 10.24400000 3.30114286
56 -0.61314286 10.24400000
57 9.04400000 -0.61314286
58 11.57257143 9.04400000
59 -0.97028571 11.57257143
60 3.70171429 -0.97028571
61 0.41657143 3.70171429
62 5.57371429 0.41657143
63 4.17371429 5.57371429
64 3.77371429 4.17371429
65 4.67371429 3.77371429
66 2.40228571 4.67371429
67 5.74514286 2.40228571
68 -0.61200000 5.74514286
69 7.84514286 -0.61200000
70 7.87371429 7.84514286
71 5.43085714 7.87371429
72 2.70285714 5.43085714
73 -1.18228571 2.70285714
74 7.77485714 -1.18228571
75 4.57485714 7.77485714
76 1.07485714 4.57485714
77 -3.02514286 1.07485714
78 -0.49657143 -3.02514286
79 -8.35371429 -0.49657143
80 2.28914286 -8.35371429
81 -3.95371429 2.28914286
82 -4.22514286 -3.95371429
83 4.53200000 -4.22514286
84 -5.49600000 4.53200000
> 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/freestat/rcomp/tmp/76lx11227794852.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/freestat/rcomp/tmp/8vhqg1227794852.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/freestat/rcomp/tmp/9gjmc1227794852.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/freestat/rcomp/tmp/10v0jd1227794852.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11ykur1227794852.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/freestat/rcomp/tmp/121eac1227794852.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/freestat/rcomp/tmp/13ftqo1227794852.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/freestat/rcomp/tmp/14mppa1227794853.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/freestat/rcomp/tmp/15vrao1227794853.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/freestat/rcomp/tmp/16geqa1227794853.tab")
+ }
>
> system("convert tmp/1u55k1227794852.ps tmp/1u55k1227794852.png")
> system("convert tmp/26dez1227794852.ps tmp/26dez1227794852.png")
> system("convert tmp/32cxn1227794852.ps tmp/32cxn1227794852.png")
> system("convert tmp/4b1v21227794852.ps tmp/4b1v21227794852.png")
> system("convert tmp/5mgzu1227794852.ps tmp/5mgzu1227794852.png")
> system("convert tmp/61xdi1227794852.ps tmp/61xdi1227794852.png")
> system("convert tmp/76lx11227794852.ps tmp/76lx11227794852.png")
> system("convert tmp/8vhqg1227794852.ps tmp/8vhqg1227794852.png")
> system("convert tmp/9gjmc1227794852.ps tmp/9gjmc1227794852.png")
> system("convert tmp/10v0jd1227794852.ps tmp/10v0jd1227794852.png")
>
>
> proc.time()
user system elapsed
3.968 2.495 4.415