R version 2.7.0 (2008-04-22)
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(105.4,109.1,107.1,111.4,110.7,114.1,117.1,121.8,118.7,127.6,126.5,129.9,127.5,128,134.6,123.5,131.8,124,135.9,127.4,142.7,127.6,141.7,128.4,153.4,131.4,145,135.1,137.7,134,148.3,144.5,152.2,147.3,169.4,150.9,168.6,148.7,161.1,141.4,174.1,138.9,179,139.8,190.6,145.6,190,147.9,181.6,148.5,174.8,151.1,180.5,157.5,196.8,167.5,193.8,172.3,197,173.5,216.3,187.5,221.4,205.5,217.9,195.1,229.7,204.5,227.4,204.5,204.2,201.7,196.6,207,198.8,206.6,207.5,210.6,190.7,211.1,201.6,215,210.5,223.9,223.5,238.2,223.8,238.9,231.2,229.6,244,232.2,234.7,222.1,250.2,221.6,265.7,227.3,287.6,221,283.3,213.6,295.4,243.4,312.3,253.8,333.8,265.3,347.7,268.2,383.2,268.5,407.1,266.9,413.6,268.4,362.7,250.8,321.9,231.2,239.4,192),dim=c(2,61),dimnames=list(c('alg_indexcijfer_grondstoffen','indexcijfer_industr_grondstoffen'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('alg_indexcijfer_grondstoffen','indexcijfer_industr_grondstoffen'),1:61))
> 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
alg_indexcijfer_grondstoffen indexcijfer_industr_grondstoffen M1 M2 M3 M4 M5
1 105.4 109.1 1 0 0 0 0
2 107.1 111.4 0 1 0 0 0
3 110.7 114.1 0 0 1 0 0
4 117.1 121.8 0 0 0 1 0
5 118.7 127.6 0 0 0 0 1
6 126.5 129.9 0 0 0 0 0
7 127.5 128.0 0 0 0 0 0
8 134.6 123.5 0 0 0 0 0
9 131.8 124.0 0 0 0 0 0
10 135.9 127.4 0 0 0 0 0
11 142.7 127.6 0 0 0 0 0
12 141.7 128.4 0 0 0 0 0
13 153.4 131.4 1 0 0 0 0
14 145.0 135.1 0 1 0 0 0
15 137.7 134.0 0 0 1 0 0
16 148.3 144.5 0 0 0 1 0
17 152.2 147.3 0 0 0 0 1
18 169.4 150.9 0 0 0 0 0
19 168.6 148.7 0 0 0 0 0
20 161.1 141.4 0 0 0 0 0
21 174.1 138.9 0 0 0 0 0
22 179.0 139.8 0 0 0 0 0
23 190.6 145.6 0 0 0 0 0
24 190.0 147.9 0 0 0 0 0
25 181.6 148.5 1 0 0 0 0
26 174.8 151.1 0 1 0 0 0
27 180.5 157.5 0 0 1 0 0
28 196.8 167.5 0 0 0 1 0
29 193.8 172.3 0 0 0 0 1
30 197.0 173.5 0 0 0 0 0
31 216.3 187.5 0 0 0 0 0
32 221.4 205.5 0 0 0 0 0
33 217.9 195.1 0 0 0 0 0
34 229.7 204.5 0 0 0 0 0
35 227.4 204.5 0 0 0 0 0
36 204.2 201.7 0 0 0 0 0
37 196.6 207.0 1 0 0 0 0
38 198.8 206.6 0 1 0 0 0
39 207.5 210.6 0 0 1 0 0
40 190.7 211.1 0 0 0 1 0
41 201.6 215.0 0 0 0 0 1
42 210.5 223.9 0 0 0 0 0
43 223.5 238.2 0 0 0 0 0
44 223.8 238.9 0 0 0 0 0
45 231.2 229.6 0 0 0 0 0
46 244.0 232.2 0 0 0 0 0
47 234.7 222.1 0 0 0 0 0
48 250.2 221.6 0 0 0 0 0
49 265.7 227.3 1 0 0 0 0
50 287.6 221.0 0 1 0 0 0
51 283.3 213.6 0 0 1 0 0
52 295.4 243.4 0 0 0 1 0
53 312.3 253.8 0 0 0 0 1
54 333.8 265.3 0 0 0 0 0
55 347.7 268.2 0 0 0 0 0
56 383.2 268.5 0 0 0 0 0
57 407.1 266.9 0 0 0 0 0
58 413.6 268.4 0 0 0 0 0
59 362.7 250.8 0 0 0 0 0
60 321.9 231.2 0 0 0 0 0
61 239.4 192.0 1 0 0 0 0
M6 M7 M8 M9 M10 M11 t
1 0 0 0 0 0 0 1
2 0 0 0 0 0 0 2
3 0 0 0 0 0 0 3
4 0 0 0 0 0 0 4
5 0 0 0 0 0 0 5
6 1 0 0 0 0 0 6
7 0 1 0 0 0 0 7
8 0 0 1 0 0 0 8
9 0 0 0 1 0 0 9
10 0 0 0 0 1 0 10
11 0 0 0 0 0 1 11
12 0 0 0 0 0 0 12
13 0 0 0 0 0 0 13
14 0 0 0 0 0 0 14
15 0 0 0 0 0 0 15
16 0 0 0 0 0 0 16
17 0 0 0 0 0 0 17
18 1 0 0 0 0 0 18
19 0 1 0 0 0 0 19
20 0 0 1 0 0 0 20
21 0 0 0 1 0 0 21
22 0 0 0 0 1 0 22
23 0 0 0 0 0 1 23
24 0 0 0 0 0 0 24
25 0 0 0 0 0 0 25
26 0 0 0 0 0 0 26
27 0 0 0 0 0 0 27
28 0 0 0 0 0 0 28
29 0 0 0 0 0 0 29
30 1 0 0 0 0 0 30
31 0 1 0 0 0 0 31
32 0 0 1 0 0 0 32
33 0 0 0 1 0 0 33
34 0 0 0 0 1 0 34
35 0 0 0 0 0 1 35
36 0 0 0 0 0 0 36
37 0 0 0 0 0 0 37
38 0 0 0 0 0 0 38
39 0 0 0 0 0 0 39
40 0 0 0 0 0 0 40
41 0 0 0 0 0 0 41
42 1 0 0 0 0 0 42
43 0 1 0 0 0 0 43
44 0 0 1 0 0 0 44
45 0 0 0 1 0 0 45
46 0 0 0 0 1 0 46
47 0 0 0 0 0 1 47
48 0 0 0 0 0 0 48
49 0 0 0 0 0 0 49
50 0 0 0 0 0 0 50
51 0 0 0 0 0 0 51
52 0 0 0 0 0 0 52
53 0 0 0 0 0 0 53
54 1 0 0 0 0 0 54
55 0 1 0 0 0 0 55
56 0 0 1 0 0 0 56
57 0 0 0 1 0 0 57
58 0 0 0 0 1 0 58
59 0 0 0 0 0 1 59
60 0 0 0 0 0 0 60
61 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) indexcijfer_industr_grondstoffen
28.6032 0.6138
M1 M2
-9.9156 -4.1081
M3 M4
-5.5796 -9.2283
M5 M6
-8.7558 -2.5987
M7 M8
1.1675 6.1968
M9 M10
14.4706 18.1185
M11 t
9.7759 2.1867
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-53.865 -16.903 3.850 13.377 75.546
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.6032 32.9578 0.868 0.3899
indexcijfer_industr_grondstoffen 0.6138 0.3131 1.960 0.0559 .
M1 -9.9156 19.1491 -0.518 0.6070
M2 -4.1081 20.1274 -0.204 0.8392
M3 -5.5796 20.0653 -0.278 0.7822
M4 -9.2283 20.3960 -0.452 0.6530
M5 -8.7558 20.5713 -0.426 0.6723
M6 -2.5987 20.7842 -0.125 0.9010
M7 1.1675 21.0278 0.056 0.9560
M8 6.1968 20.9015 0.296 0.7682
M9 14.4706 20.3271 0.712 0.4801
M10 18.1185 20.3789 0.889 0.3785
M11 9.7759 20.0336 0.488 0.6278
t 2.1867 0.8571 2.551 0.0140 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 31.51 on 47 degrees of freedom
Multiple R-squared: 0.8629, Adjusted R-squared: 0.8249
F-statistic: 22.75 on 13 and 47 DF, p-value: 5.939e-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,] 5.190769e-03 1.038154e-02 0.9948092
[2,] 1.142261e-03 2.284522e-03 0.9988577
[3,] 1.906249e-04 3.812499e-04 0.9998094
[4,] 2.702352e-05 5.404705e-05 0.9999730
[5,] 1.191244e-05 2.382488e-05 0.9999881
[6,] 2.086529e-06 4.173058e-06 0.9999979
[7,] 6.193332e-07 1.238666e-06 0.9999994
[8,] 2.116425e-07 4.232849e-07 0.9999998
[9,] 7.291144e-08 1.458229e-07 0.9999999
[10,] 2.304544e-08 4.609089e-08 1.0000000
[11,] 3.974753e-09 7.949505e-09 1.0000000
[12,] 1.984508e-09 3.969017e-09 1.0000000
[13,] 5.467775e-10 1.093555e-09 1.0000000
[14,] 4.540016e-10 9.080032e-10 1.0000000
[15,] 9.782929e-10 1.956586e-09 1.0000000
[16,] 6.639171e-10 1.327834e-09 1.0000000
[17,] 9.024995e-10 1.804999e-09 1.0000000
[18,] 1.790122e-09 3.580245e-09 1.0000000
[19,] 5.120765e-08 1.024153e-07 0.9999999
[20,] 4.705822e-05 9.411644e-05 0.9999529
[21,] 2.157026e-03 4.314052e-03 0.9978430
[22,] 2.273771e-03 4.547541e-03 0.9977262
[23,] 1.297528e-03 2.595056e-03 0.9987025
[24,] 1.077023e-02 2.154045e-02 0.9892298
[25,] 5.689429e-02 1.137886e-01 0.9431057
[26,] 3.919704e-01 7.839408e-01 0.6080296
[27,] 4.615792e-01 9.231584e-01 0.5384208
[28,] 4.289375e-01 8.578750e-01 0.5710625
> postscript(file="/var/www/html/rcomp/tmp/1tnvz1227789357.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/2fwwd1227789357.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/3u7u91227789357.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/4ehi51227789357.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/5kqn01227789357.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 = 61
Frequency = 1
1 2 3 4 5 6
17.5548529 9.8486958 11.0760453 14.2114397 9.5918390 7.6361545
7 8 9 10 11 12
3.8495262 6.4957868 -7.0716719 -10.8934561 1.9396713 8.0377520
13 14 15 16 17 18
25.6251006 6.9595555 -0.3804704 5.2361480 4.7580931 11.4044054
19 20 21 22 23 24
6.0019317 -4.2330317 -0.1589446 -1.6461074 12.5494679 18.1267758
25 26 27 28 29 30
17.0873610 0.6970493 1.7531591 13.3767017 4.7709496 -1.1095014
31 32 33 34 35 36
3.6436777 -9.5216550 -17.0981641 -16.9030398 -13.0471427 -26.9392071
37 38 39 40 41 42
-30.0637102 -35.6124761 -30.0831298 -45.7280256 -39.8813139 -44.7883991
43 44 45 46 47 48
-46.5193746 -53.8651266 -51.2168691 -45.8475745 -42.7918066 -19.3957227
49 50 51 52 53 54
0.3342347 18.1071755 17.6343957 12.9037363 20.7604321 26.8573406
55 56 57 58 59 60
33.0242390 61.1240264 75.5456497 75.2901778 41.3498101 20.1704021
61
-30.5378389
> postscript(file="/var/www/html/rcomp/tmp/60rdy1227789357.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 17.5548529 NA
1 9.8486958 17.5548529
2 11.0760453 9.8486958
3 14.2114397 11.0760453
4 9.5918390 14.2114397
5 7.6361545 9.5918390
6 3.8495262 7.6361545
7 6.4957868 3.8495262
8 -7.0716719 6.4957868
9 -10.8934561 -7.0716719
10 1.9396713 -10.8934561
11 8.0377520 1.9396713
12 25.6251006 8.0377520
13 6.9595555 25.6251006
14 -0.3804704 6.9595555
15 5.2361480 -0.3804704
16 4.7580931 5.2361480
17 11.4044054 4.7580931
18 6.0019317 11.4044054
19 -4.2330317 6.0019317
20 -0.1589446 -4.2330317
21 -1.6461074 -0.1589446
22 12.5494679 -1.6461074
23 18.1267758 12.5494679
24 17.0873610 18.1267758
25 0.6970493 17.0873610
26 1.7531591 0.6970493
27 13.3767017 1.7531591
28 4.7709496 13.3767017
29 -1.1095014 4.7709496
30 3.6436777 -1.1095014
31 -9.5216550 3.6436777
32 -17.0981641 -9.5216550
33 -16.9030398 -17.0981641
34 -13.0471427 -16.9030398
35 -26.9392071 -13.0471427
36 -30.0637102 -26.9392071
37 -35.6124761 -30.0637102
38 -30.0831298 -35.6124761
39 -45.7280256 -30.0831298
40 -39.8813139 -45.7280256
41 -44.7883991 -39.8813139
42 -46.5193746 -44.7883991
43 -53.8651266 -46.5193746
44 -51.2168691 -53.8651266
45 -45.8475745 -51.2168691
46 -42.7918066 -45.8475745
47 -19.3957227 -42.7918066
48 0.3342347 -19.3957227
49 18.1071755 0.3342347
50 17.6343957 18.1071755
51 12.9037363 17.6343957
52 20.7604321 12.9037363
53 26.8573406 20.7604321
54 33.0242390 26.8573406
55 61.1240264 33.0242390
56 75.5456497 61.1240264
57 75.2901778 75.5456497
58 41.3498101 75.2901778
59 20.1704021 41.3498101
60 -30.5378389 20.1704021
61 NA -30.5378389
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 9.8486958 17.5548529
[2,] 11.0760453 9.8486958
[3,] 14.2114397 11.0760453
[4,] 9.5918390 14.2114397
[5,] 7.6361545 9.5918390
[6,] 3.8495262 7.6361545
[7,] 6.4957868 3.8495262
[8,] -7.0716719 6.4957868
[9,] -10.8934561 -7.0716719
[10,] 1.9396713 -10.8934561
[11,] 8.0377520 1.9396713
[12,] 25.6251006 8.0377520
[13,] 6.9595555 25.6251006
[14,] -0.3804704 6.9595555
[15,] 5.2361480 -0.3804704
[16,] 4.7580931 5.2361480
[17,] 11.4044054 4.7580931
[18,] 6.0019317 11.4044054
[19,] -4.2330317 6.0019317
[20,] -0.1589446 -4.2330317
[21,] -1.6461074 -0.1589446
[22,] 12.5494679 -1.6461074
[23,] 18.1267758 12.5494679
[24,] 17.0873610 18.1267758
[25,] 0.6970493 17.0873610
[26,] 1.7531591 0.6970493
[27,] 13.3767017 1.7531591
[28,] 4.7709496 13.3767017
[29,] -1.1095014 4.7709496
[30,] 3.6436777 -1.1095014
[31,] -9.5216550 3.6436777
[32,] -17.0981641 -9.5216550
[33,] -16.9030398 -17.0981641
[34,] -13.0471427 -16.9030398
[35,] -26.9392071 -13.0471427
[36,] -30.0637102 -26.9392071
[37,] -35.6124761 -30.0637102
[38,] -30.0831298 -35.6124761
[39,] -45.7280256 -30.0831298
[40,] -39.8813139 -45.7280256
[41,] -44.7883991 -39.8813139
[42,] -46.5193746 -44.7883991
[43,] -53.8651266 -46.5193746
[44,] -51.2168691 -53.8651266
[45,] -45.8475745 -51.2168691
[46,] -42.7918066 -45.8475745
[47,] -19.3957227 -42.7918066
[48,] 0.3342347 -19.3957227
[49,] 18.1071755 0.3342347
[50,] 17.6343957 18.1071755
[51,] 12.9037363 17.6343957
[52,] 20.7604321 12.9037363
[53,] 26.8573406 20.7604321
[54,] 33.0242390 26.8573406
[55,] 61.1240264 33.0242390
[56,] 75.5456497 61.1240264
[57,] 75.2901778 75.5456497
[58,] 41.3498101 75.2901778
[59,] 20.1704021 41.3498101
[60,] -30.5378389 20.1704021
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 9.8486958 17.5548529
2 11.0760453 9.8486958
3 14.2114397 11.0760453
4 9.5918390 14.2114397
5 7.6361545 9.5918390
6 3.8495262 7.6361545
7 6.4957868 3.8495262
8 -7.0716719 6.4957868
9 -10.8934561 -7.0716719
10 1.9396713 -10.8934561
11 8.0377520 1.9396713
12 25.6251006 8.0377520
13 6.9595555 25.6251006
14 -0.3804704 6.9595555
15 5.2361480 -0.3804704
16 4.7580931 5.2361480
17 11.4044054 4.7580931
18 6.0019317 11.4044054
19 -4.2330317 6.0019317
20 -0.1589446 -4.2330317
21 -1.6461074 -0.1589446
22 12.5494679 -1.6461074
23 18.1267758 12.5494679
24 17.0873610 18.1267758
25 0.6970493 17.0873610
26 1.7531591 0.6970493
27 13.3767017 1.7531591
28 4.7709496 13.3767017
29 -1.1095014 4.7709496
30 3.6436777 -1.1095014
31 -9.5216550 3.6436777
32 -17.0981641 -9.5216550
33 -16.9030398 -17.0981641
34 -13.0471427 -16.9030398
35 -26.9392071 -13.0471427
36 -30.0637102 -26.9392071
37 -35.6124761 -30.0637102
38 -30.0831298 -35.6124761
39 -45.7280256 -30.0831298
40 -39.8813139 -45.7280256
41 -44.7883991 -39.8813139
42 -46.5193746 -44.7883991
43 -53.8651266 -46.5193746
44 -51.2168691 -53.8651266
45 -45.8475745 -51.2168691
46 -42.7918066 -45.8475745
47 -19.3957227 -42.7918066
48 0.3342347 -19.3957227
49 18.1071755 0.3342347
50 17.6343957 18.1071755
51 12.9037363 17.6343957
52 20.7604321 12.9037363
53 26.8573406 20.7604321
54 33.0242390 26.8573406
55 61.1240264 33.0242390
56 75.5456497 61.1240264
57 75.2901778 75.5456497
58 41.3498101 75.2901778
59 20.1704021 41.3498101
60 -30.5378389 20.1704021
> 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/73lwb1227789357.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/8pduw1227789357.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/99dd91227789357.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/10106s1227789357.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/117snx1227789357.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/12okoj1227789357.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/136sns1227789357.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/14cbzk1227789357.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/15mu9x1227789357.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/16kqnq1227789357.tab")
+ }
>
> system("convert tmp/1tnvz1227789357.ps tmp/1tnvz1227789357.png")
> system("convert tmp/2fwwd1227789357.ps tmp/2fwwd1227789357.png")
> system("convert tmp/3u7u91227789357.ps tmp/3u7u91227789357.png")
> system("convert tmp/4ehi51227789357.ps tmp/4ehi51227789357.png")
> system("convert tmp/5kqn01227789357.ps tmp/5kqn01227789357.png")
> system("convert tmp/60rdy1227789357.ps tmp/60rdy1227789357.png")
> system("convert tmp/73lwb1227789357.ps tmp/73lwb1227789357.png")
> system("convert tmp/8pduw1227789357.ps tmp/8pduw1227789357.png")
> system("convert tmp/99dd91227789357.ps tmp/99dd91227789357.png")
> system("convert tmp/10106s1227789357.ps tmp/10106s1227789357.png")
>
>
> proc.time()
user system elapsed
4.983 2.780 5.372