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(16643
+ ,16196.7
+ ,18252.1
+ ,17570.4
+ ,15836.8
+ ,89.1
+ ,17729
+ ,16643
+ ,16196.7
+ ,18252.1
+ ,17570.4
+ ,82.6
+ ,16446.1
+ ,17729
+ ,16643
+ ,16196.7
+ ,18252.1
+ ,102.7
+ ,15993.8
+ ,16446.1
+ ,17729
+ ,16643
+ ,16196.7
+ ,91.8
+ ,16373.5
+ ,15993.8
+ ,16446.1
+ ,17729
+ ,16643
+ ,94.1
+ ,17842.2
+ ,16373.5
+ ,15993.8
+ ,16446.1
+ ,17729
+ ,103.1
+ ,22321.5
+ ,17842.2
+ ,16373.5
+ ,15993.8
+ ,16446.1
+ ,93.2
+ ,22786.7
+ ,22321.5
+ ,17842.2
+ ,16373.5
+ ,15993.8
+ ,91
+ ,18274.1
+ ,22786.7
+ ,22321.5
+ ,17842.2
+ ,16373.5
+ ,94.3
+ ,22392.9
+ ,18274.1
+ ,22786.7
+ ,22321.5
+ ,17842.2
+ ,99.4
+ ,23899.3
+ ,22392.9
+ ,18274.1
+ ,22786.7
+ ,22321.5
+ ,115.7
+ ,21343.5
+ ,23899.3
+ ,22392.9
+ ,18274.1
+ ,22786.7
+ ,116.8
+ ,22952.3
+ ,21343.5
+ ,23899.3
+ ,22392.9
+ ,18274.1
+ ,99.8
+ ,21374.4
+ ,22952.3
+ ,21343.5
+ ,23899.3
+ ,22392.9
+ ,96
+ ,21164.1
+ ,21374.4
+ ,22952.3
+ ,21343.5
+ ,23899.3
+ ,115.9
+ ,20906.5
+ ,21164.1
+ ,21374.4
+ ,22952.3
+ ,21343.5
+ ,109.1
+ ,17877.4
+ ,20906.5
+ ,21164.1
+ ,21374.4
+ ,22952.3
+ ,117.3
+ ,20664.3
+ ,17877.4
+ ,20906.5
+ ,21164.1
+ ,21374.4
+ ,109.8
+ ,22160
+ ,20664.3
+ ,17877.4
+ ,20906.5
+ ,21164.1
+ ,112.8
+ ,19813.6
+ ,22160
+ ,20664.3
+ ,17877.4
+ ,20906.5
+ ,110.7
+ ,17735.4
+ ,19813.6
+ ,22160
+ ,20664.3
+ ,17877.4
+ ,100
+ ,19640.2
+ ,17735.4
+ ,19813.6
+ ,22160
+ ,20664.3
+ ,113.3
+ ,20844.4
+ ,19640.2
+ ,17735.4
+ ,19813.6
+ ,22160
+ ,122.4
+ ,19823.1
+ ,20844.4
+ ,19640.2
+ ,17735.4
+ ,19813.6
+ ,112.5
+ ,18594.6
+ ,19823.1
+ ,20844.4
+ ,19640.2
+ ,17735.4
+ ,104.2
+ ,21350.6
+ ,18594.6
+ ,19823.1
+ ,20844.4
+ ,19640.2
+ ,92.5
+ ,18574.1
+ ,21350.6
+ ,18594.6
+ ,19823.1
+ ,20844.4
+ ,117.2
+ ,18924.2
+ ,18574.1
+ ,21350.6
+ ,18594.6
+ ,19823.1
+ ,109.3
+ ,17343.4
+ ,18924.2
+ ,18574.1
+ ,21350.6
+ ,18594.6
+ ,106.1
+ ,19961.2
+ ,17343.4
+ ,18924.2
+ ,18574.1
+ ,21350.6
+ ,118.8
+ ,19932.1
+ ,19961.2
+ ,17343.4
+ ,18924.2
+ ,18574.1
+ ,105.3
+ ,19464.6
+ ,19932.1
+ ,19961.2
+ ,17343.4
+ ,18924.2
+ ,106
+ ,16165.4
+ ,19464.6
+ ,19932.1
+ ,19961.2
+ ,17343.4
+ ,102
+ ,17574.9
+ ,16165.4
+ ,19464.6
+ ,19932.1
+ ,19961.2
+ ,112.9
+ ,19795.4
+ ,17574.9
+ ,16165.4
+ ,19464.6
+ ,19932.1
+ ,116.5
+ ,19439.5
+ ,19795.4
+ ,17574.9
+ ,16165.4
+ ,19464.6
+ ,114.8
+ ,17170
+ ,19439.5
+ ,19795.4
+ ,17574.9
+ ,16165.4
+ ,100.5
+ ,21072.4
+ ,17170
+ ,19439.5
+ ,19795.4
+ ,17574.9
+ ,85.4
+ ,17751.8
+ ,21072.4
+ ,17170
+ ,19439.5
+ ,19795.4
+ ,114.6
+ ,17515.5
+ ,17751.8
+ ,21072.4
+ ,17170
+ ,19439.5
+ ,109.9
+ ,18040.3
+ ,17515.5
+ ,17751.8
+ ,21072.4
+ ,17170
+ ,100.7
+ ,19090.1
+ ,18040.3
+ ,17515.5
+ ,17751.8
+ ,21072.4
+ ,115.5
+ ,17746.5
+ ,19090.1
+ ,18040.3
+ ,17515.5
+ ,17751.8
+ ,100.7
+ ,19202.1
+ ,17746.5
+ ,19090.1
+ ,18040.3
+ ,17515.5
+ ,99
+ ,15141.6
+ ,19202.1
+ ,17746.5
+ ,19090.1
+ ,18040.3
+ ,102.3
+ ,16258.1
+ ,15141.6
+ ,19202.1
+ ,17746.5
+ ,19090.1
+ ,108.8
+ ,18586.5
+ ,16258.1
+ ,15141.6
+ ,19202.1
+ ,17746.5
+ ,105.9
+ ,17209.4
+ ,18586.5
+ ,16258.1
+ ,15141.6
+ ,19202.1
+ ,113.2
+ ,17838.7
+ ,17209.4
+ ,18586.5
+ ,16258.1
+ ,15141.6
+ ,95.7
+ ,19123.5
+ ,17838.7
+ ,17209.4
+ ,18586.5
+ ,16258.1
+ ,80.9
+ ,16583.6
+ ,19123.5
+ ,17838.7
+ ,17209.4
+ ,18586.5
+ ,113.9
+ ,15991.2
+ ,16583.6
+ ,19123.5
+ ,17838.7
+ ,17209.4
+ ,98.1
+ ,16704.4
+ ,15991.2
+ ,16583.6
+ ,19123.5
+ ,17838.7
+ ,102.8
+ ,17420.4
+ ,16704.4
+ ,15991.2
+ ,16583.6
+ ,19123.5
+ ,104.7
+ ,17872
+ ,17420.4
+ ,16704.4
+ ,15991.2
+ ,16583.6
+ ,95.9
+ ,17823.2
+ ,17872
+ ,17420.4
+ ,16704.4
+ ,15991.2
+ ,94.6)
+ ,dim=c(6
+ ,56)
+ ,dimnames=list(c('uitvoer'
+ ,'uitvoer1'
+ ,'uitvoer2'
+ ,'uitvoer3'
+ ,'uitvoer4'
+ ,'indprod')
+ ,1:56))
> y <- array(NA,dim=c(6,56),dimnames=list(c('uitvoer','uitvoer1','uitvoer2','uitvoer3','uitvoer4','indprod'),1:56))
> 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
uitvoer uitvoer1 uitvoer2 uitvoer3 uitvoer4 indprod M1 M2 M3 M4 M5 M6 M7 M8
1 16643.0 16196.7 18252.1 17570.4 15836.8 89.1 1 0 0 0 0 0 0 0
2 17729.0 16643.0 16196.7 18252.1 17570.4 82.6 0 1 0 0 0 0 0 0
3 16446.1 17729.0 16643.0 16196.7 18252.1 102.7 0 0 1 0 0 0 0 0
4 15993.8 16446.1 17729.0 16643.0 16196.7 91.8 0 0 0 1 0 0 0 0
5 16373.5 15993.8 16446.1 17729.0 16643.0 94.1 0 0 0 0 1 0 0 0
6 17842.2 16373.5 15993.8 16446.1 17729.0 103.1 0 0 0 0 0 1 0 0
7 22321.5 17842.2 16373.5 15993.8 16446.1 93.2 0 0 0 0 0 0 1 0
8 22786.7 22321.5 17842.2 16373.5 15993.8 91.0 0 0 0 0 0 0 0 1
9 18274.1 22786.7 22321.5 17842.2 16373.5 94.3 0 0 0 0 0 0 0 0
10 22392.9 18274.1 22786.7 22321.5 17842.2 99.4 0 0 0 0 0 0 0 0
11 23899.3 22392.9 18274.1 22786.7 22321.5 115.7 0 0 0 0 0 0 0 0
12 21343.5 23899.3 22392.9 18274.1 22786.7 116.8 0 0 0 0 0 0 0 0
13 22952.3 21343.5 23899.3 22392.9 18274.1 99.8 1 0 0 0 0 0 0 0
14 21374.4 22952.3 21343.5 23899.3 22392.9 96.0 0 1 0 0 0 0 0 0
15 21164.1 21374.4 22952.3 21343.5 23899.3 115.9 0 0 1 0 0 0 0 0
16 20906.5 21164.1 21374.4 22952.3 21343.5 109.1 0 0 0 1 0 0 0 0
17 17877.4 20906.5 21164.1 21374.4 22952.3 117.3 0 0 0 0 1 0 0 0
18 20664.3 17877.4 20906.5 21164.1 21374.4 109.8 0 0 0 0 0 1 0 0
19 22160.0 20664.3 17877.4 20906.5 21164.1 112.8 0 0 0 0 0 0 1 0
20 19813.6 22160.0 20664.3 17877.4 20906.5 110.7 0 0 0 0 0 0 0 1
21 17735.4 19813.6 22160.0 20664.3 17877.4 100.0 0 0 0 0 0 0 0 0
22 19640.2 17735.4 19813.6 22160.0 20664.3 113.3 0 0 0 0 0 0 0 0
23 20844.4 19640.2 17735.4 19813.6 22160.0 122.4 0 0 0 0 0 0 0 0
24 19823.1 20844.4 19640.2 17735.4 19813.6 112.5 0 0 0 0 0 0 0 0
25 18594.6 19823.1 20844.4 19640.2 17735.4 104.2 1 0 0 0 0 0 0 0
26 21350.6 18594.6 19823.1 20844.4 19640.2 92.5 0 1 0 0 0 0 0 0
27 18574.1 21350.6 18594.6 19823.1 20844.4 117.2 0 0 1 0 0 0 0 0
28 18924.2 18574.1 21350.6 18594.6 19823.1 109.3 0 0 0 1 0 0 0 0
29 17343.4 18924.2 18574.1 21350.6 18594.6 106.1 0 0 0 0 1 0 0 0
30 19961.2 17343.4 18924.2 18574.1 21350.6 118.8 0 0 0 0 0 1 0 0
31 19932.1 19961.2 17343.4 18924.2 18574.1 105.3 0 0 0 0 0 0 1 0
32 19464.6 19932.1 19961.2 17343.4 18924.2 106.0 0 0 0 0 0 0 0 1
33 16165.4 19464.6 19932.1 19961.2 17343.4 102.0 0 0 0 0 0 0 0 0
34 17574.9 16165.4 19464.6 19932.1 19961.2 112.9 0 0 0 0 0 0 0 0
35 19795.4 17574.9 16165.4 19464.6 19932.1 116.5 0 0 0 0 0 0 0 0
36 19439.5 19795.4 17574.9 16165.4 19464.6 114.8 0 0 0 0 0 0 0 0
37 17170.0 19439.5 19795.4 17574.9 16165.4 100.5 1 0 0 0 0 0 0 0
38 21072.4 17170.0 19439.5 19795.4 17574.9 85.4 0 1 0 0 0 0 0 0
39 17751.8 21072.4 17170.0 19439.5 19795.4 114.6 0 0 1 0 0 0 0 0
40 17515.5 17751.8 21072.4 17170.0 19439.5 109.9 0 0 0 1 0 0 0 0
41 18040.3 17515.5 17751.8 21072.4 17170.0 100.7 0 0 0 0 1 0 0 0
42 19090.1 18040.3 17515.5 17751.8 21072.4 115.5 0 0 0 0 0 1 0 0
43 17746.5 19090.1 18040.3 17515.5 17751.8 100.7 0 0 0 0 0 0 1 0
44 19202.1 17746.5 19090.1 18040.3 17515.5 99.0 0 0 0 0 0 0 0 1
45 15141.6 19202.1 17746.5 19090.1 18040.3 102.3 0 0 0 0 0 0 0 0
46 16258.1 15141.6 19202.1 17746.5 19090.1 108.8 0 0 0 0 0 0 0 0
47 18586.5 16258.1 15141.6 19202.1 17746.5 105.9 0 0 0 0 0 0 0 0
48 17209.4 18586.5 16258.1 15141.6 19202.1 113.2 0 0 0 0 0 0 0 0
49 17838.7 17209.4 18586.5 16258.1 15141.6 95.7 1 0 0 0 0 0 0 0
50 19123.5 17838.7 17209.4 18586.5 16258.1 80.9 0 1 0 0 0 0 0 0
51 16583.6 19123.5 17838.7 17209.4 18586.5 113.9 0 0 1 0 0 0 0 0
52 15991.2 16583.6 19123.5 17838.7 17209.4 98.1 0 0 0 1 0 0 0 0
53 16704.4 15991.2 16583.6 19123.5 17838.7 102.8 0 0 0 0 1 0 0 0
54 17420.4 16704.4 15991.2 16583.6 19123.5 104.7 0 0 0 0 0 1 0 0
55 17872.0 17420.4 16704.4 15991.2 16583.6 95.9 0 0 0 0 0 0 1 0
56 17823.2 17872.0 17420.4 16704.4 15991.2 94.6 0 0 0 0 0 0 0 1
M9 M10 M11 t
1 0 0 0 1
2 0 0 0 2
3 0 0 0 3
4 0 0 0 4
5 0 0 0 5
6 0 0 0 6
7 0 0 0 7
8 0 0 0 8
9 1 0 0 9
10 0 1 0 10
11 0 0 1 11
12 0 0 0 12
13 0 0 0 13
14 0 0 0 14
15 0 0 0 15
16 0 0 0 16
17 0 0 0 17
18 0 0 0 18
19 0 0 0 19
20 0 0 0 20
21 1 0 0 21
22 0 1 0 22
23 0 0 1 23
24 0 0 0 24
25 0 0 0 25
26 0 0 0 26
27 0 0 0 27
28 0 0 0 28
29 0 0 0 29
30 0 0 0 30
31 0 0 0 31
32 0 0 0 32
33 1 0 0 33
34 0 1 0 34
35 0 0 1 35
36 0 0 0 36
37 0 0 0 37
38 0 0 0 38
39 0 0 0 39
40 0 0 0 40
41 0 0 0 41
42 0 0 0 42
43 0 0 0 43
44 0 0 0 44
45 1 0 0 45
46 0 1 0 46
47 0 0 1 47
48 0 0 0 48
49 0 0 0 49
50 0 0 0 50
51 0 0 0 51
52 0 0 0 52
53 0 0 0 53
54 0 0 0 54
55 0 0 0 55
56 0 0 0 56
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) uitvoer1 uitvoer2 uitvoer3 uitvoer4 indprod
6470.9462 0.3224 0.3836 0.4096 -0.3752 3.9111
M1 M2 M3 M4 M5 M6
-2876.5130 -583.2869 -1903.6624 -2486.1245 -2880.6311 519.8767
M7 M8 M9 M10 M11 t
598.5011 -393.8065 -5232.1043 -1624.0076 1361.0215 -23.7051
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2287.3 -635.7 -104.9 561.2 2638.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6470.9462 3486.2938 1.856 0.07120 .
uitvoer1 0.3224 0.1436 2.245 0.03066 *
uitvoer2 0.3836 0.1428 2.687 0.01065 *
uitvoer3 0.4096 0.1411 2.903 0.00612 **
uitvoer4 -0.3752 0.2620 -1.432 0.16036
indprod 3.9111 60.8227 0.064 0.94907
M1 -2876.5130 1019.9424 -2.820 0.00758 **
M2 -583.2869 1496.6099 -0.390 0.69891
M3 -1903.6624 777.9012 -2.447 0.01913 *
M4 -2486.1245 937.0171 -2.653 0.01157 *
M5 -2880.6311 1034.5733 -2.784 0.00831 **
M6 519.8767 911.7697 0.570 0.57191
M7 598.5011 868.6644 0.689 0.49502
M8 -393.8065 843.5878 -0.467 0.64329
M9 -5232.1043 1003.4002 -5.214 6.78e-06 ***
M10 -1624.0076 1174.3880 -1.383 0.17478
M11 1361.0215 1029.8005 1.322 0.19419
t -23.7051 11.8782 -1.996 0.05317 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1050 on 38 degrees of freedom
Multiple R-squared: 0.822, Adjusted R-squared: 0.7423
F-statistic: 10.32 on 17 and 38 DF, p-value: 1.751e-09
> 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.8773987 0.24520267 0.12260134
[2,] 0.9676789 0.06464220 0.03232110
[3,] 0.9474269 0.10514614 0.05257307
[4,] 0.9402030 0.11959404 0.05979702
[5,] 0.9002188 0.19956241 0.09978120
[6,] 0.9235787 0.15284251 0.07642125
[7,] 0.8846693 0.23066143 0.11533071
[8,] 0.8501595 0.29968105 0.14984052
[9,] 0.9157553 0.16848935 0.08424467
[10,] 0.9432656 0.11346878 0.05673439
[11,] 0.9669618 0.06607644 0.03303822
[12,] 0.9520027 0.09599464 0.04799732
[13,] 0.9556841 0.08863176 0.04431588
[14,] 0.9644754 0.07104923 0.03552461
[15,] 0.9887954 0.02240917 0.01120459
> postscript(file="/var/www/html/rcomp/tmp/1wmz11258481536.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/24lpa1258481536.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/3fm0p1258481536.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/47spz1258481536.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/5ycsv1258481536.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 = 56
Frequency = 1
1 2 3 4 5 6
-755.37121 -897.65229 -338.77195 -1099.24071 50.27649 -908.99559
7 8 9 10 11 12
2638.85120 1795.83195 -195.17386 312.03618 686.60974 -531.62237
13 14 15 16 17 18
909.85578 -1532.62273 1026.92505 457.34394 -771.96919 -762.96912
19 20 21 22 23 24
956.21362 -773.24432 -42.88716 228.65216 141.32910 -604.61054
25 26 27 28 29 30
-593.04641 948.48769 -127.68326 817.39334 -970.19532 767.77226
31 32 33 34 35 36
-686.12514 -356.37754 -281.39697 -261.78819 -24.88954 929.91127
37 38 39 40 41 42
-935.70012 1243.78409 -255.58810 502.25186 381.38414 742.21023
43 44 45 46 47 48
-2287.27956 -82.17716 519.45799 -278.90015 -803.04929 206.32164
49 50 51 52 53 54
1374.26196 238.00323 -304.88175 -677.74844 1310.50389 161.98222
55 56
-621.66011 -584.03293
> postscript(file="/var/www/html/rcomp/tmp/6smbz1258481536.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 = 56
Frequency = 1
lag(myerror, k = 1) myerror
0 -755.37121 NA
1 -897.65229 -755.37121
2 -338.77195 -897.65229
3 -1099.24071 -338.77195
4 50.27649 -1099.24071
5 -908.99559 50.27649
6 2638.85120 -908.99559
7 1795.83195 2638.85120
8 -195.17386 1795.83195
9 312.03618 -195.17386
10 686.60974 312.03618
11 -531.62237 686.60974
12 909.85578 -531.62237
13 -1532.62273 909.85578
14 1026.92505 -1532.62273
15 457.34394 1026.92505
16 -771.96919 457.34394
17 -762.96912 -771.96919
18 956.21362 -762.96912
19 -773.24432 956.21362
20 -42.88716 -773.24432
21 228.65216 -42.88716
22 141.32910 228.65216
23 -604.61054 141.32910
24 -593.04641 -604.61054
25 948.48769 -593.04641
26 -127.68326 948.48769
27 817.39334 -127.68326
28 -970.19532 817.39334
29 767.77226 -970.19532
30 -686.12514 767.77226
31 -356.37754 -686.12514
32 -281.39697 -356.37754
33 -261.78819 -281.39697
34 -24.88954 -261.78819
35 929.91127 -24.88954
36 -935.70012 929.91127
37 1243.78409 -935.70012
38 -255.58810 1243.78409
39 502.25186 -255.58810
40 381.38414 502.25186
41 742.21023 381.38414
42 -2287.27956 742.21023
43 -82.17716 -2287.27956
44 519.45799 -82.17716
45 -278.90015 519.45799
46 -803.04929 -278.90015
47 206.32164 -803.04929
48 1374.26196 206.32164
49 238.00323 1374.26196
50 -304.88175 238.00323
51 -677.74844 -304.88175
52 1310.50389 -677.74844
53 161.98222 1310.50389
54 -621.66011 161.98222
55 -584.03293 -621.66011
56 NA -584.03293
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -897.65229 -755.37121
[2,] -338.77195 -897.65229
[3,] -1099.24071 -338.77195
[4,] 50.27649 -1099.24071
[5,] -908.99559 50.27649
[6,] 2638.85120 -908.99559
[7,] 1795.83195 2638.85120
[8,] -195.17386 1795.83195
[9,] 312.03618 -195.17386
[10,] 686.60974 312.03618
[11,] -531.62237 686.60974
[12,] 909.85578 -531.62237
[13,] -1532.62273 909.85578
[14,] 1026.92505 -1532.62273
[15,] 457.34394 1026.92505
[16,] -771.96919 457.34394
[17,] -762.96912 -771.96919
[18,] 956.21362 -762.96912
[19,] -773.24432 956.21362
[20,] -42.88716 -773.24432
[21,] 228.65216 -42.88716
[22,] 141.32910 228.65216
[23,] -604.61054 141.32910
[24,] -593.04641 -604.61054
[25,] 948.48769 -593.04641
[26,] -127.68326 948.48769
[27,] 817.39334 -127.68326
[28,] -970.19532 817.39334
[29,] 767.77226 -970.19532
[30,] -686.12514 767.77226
[31,] -356.37754 -686.12514
[32,] -281.39697 -356.37754
[33,] -261.78819 -281.39697
[34,] -24.88954 -261.78819
[35,] 929.91127 -24.88954
[36,] -935.70012 929.91127
[37,] 1243.78409 -935.70012
[38,] -255.58810 1243.78409
[39,] 502.25186 -255.58810
[40,] 381.38414 502.25186
[41,] 742.21023 381.38414
[42,] -2287.27956 742.21023
[43,] -82.17716 -2287.27956
[44,] 519.45799 -82.17716
[45,] -278.90015 519.45799
[46,] -803.04929 -278.90015
[47,] 206.32164 -803.04929
[48,] 1374.26196 206.32164
[49,] 238.00323 1374.26196
[50,] -304.88175 238.00323
[51,] -677.74844 -304.88175
[52,] 1310.50389 -677.74844
[53,] 161.98222 1310.50389
[54,] -621.66011 161.98222
[55,] -584.03293 -621.66011
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -897.65229 -755.37121
2 -338.77195 -897.65229
3 -1099.24071 -338.77195
4 50.27649 -1099.24071
5 -908.99559 50.27649
6 2638.85120 -908.99559
7 1795.83195 2638.85120
8 -195.17386 1795.83195
9 312.03618 -195.17386
10 686.60974 312.03618
11 -531.62237 686.60974
12 909.85578 -531.62237
13 -1532.62273 909.85578
14 1026.92505 -1532.62273
15 457.34394 1026.92505
16 -771.96919 457.34394
17 -762.96912 -771.96919
18 956.21362 -762.96912
19 -773.24432 956.21362
20 -42.88716 -773.24432
21 228.65216 -42.88716
22 141.32910 228.65216
23 -604.61054 141.32910
24 -593.04641 -604.61054
25 948.48769 -593.04641
26 -127.68326 948.48769
27 817.39334 -127.68326
28 -970.19532 817.39334
29 767.77226 -970.19532
30 -686.12514 767.77226
31 -356.37754 -686.12514
32 -281.39697 -356.37754
33 -261.78819 -281.39697
34 -24.88954 -261.78819
35 929.91127 -24.88954
36 -935.70012 929.91127
37 1243.78409 -935.70012
38 -255.58810 1243.78409
39 502.25186 -255.58810
40 381.38414 502.25186
41 742.21023 381.38414
42 -2287.27956 742.21023
43 -82.17716 -2287.27956
44 519.45799 -82.17716
45 -278.90015 519.45799
46 -803.04929 -278.90015
47 206.32164 -803.04929
48 1374.26196 206.32164
49 238.00323 1374.26196
50 -304.88175 238.00323
51 -677.74844 -304.88175
52 1310.50389 -677.74844
53 161.98222 1310.50389
54 -621.66011 161.98222
55 -584.03293 -621.66011
> 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/7w94a1258481536.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/8o6cn1258481536.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/9l0na1258481536.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/1076zu1258481536.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/11ky9j1258481536.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/12an151258481536.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/1368dh1258481536.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/142wt11258481536.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/15buuv1258481536.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/16pcrg1258481536.tab")
+ }
>
> system("convert tmp/1wmz11258481536.ps tmp/1wmz11258481536.png")
> system("convert tmp/24lpa1258481536.ps tmp/24lpa1258481536.png")
> system("convert tmp/3fm0p1258481536.ps tmp/3fm0p1258481536.png")
> system("convert tmp/47spz1258481536.ps tmp/47spz1258481536.png")
> system("convert tmp/5ycsv1258481536.ps tmp/5ycsv1258481536.png")
> system("convert tmp/6smbz1258481536.ps tmp/6smbz1258481536.png")
> system("convert tmp/7w94a1258481536.ps tmp/7w94a1258481536.png")
> system("convert tmp/8o6cn1258481536.ps tmp/8o6cn1258481536.png")
> system("convert tmp/9l0na1258481536.ps tmp/9l0na1258481536.png")
> system("convert tmp/1076zu1258481536.ps tmp/1076zu1258481536.png")
>
>
> proc.time()
user system elapsed
2.397 1.568 3.284