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(4
+ ,7.2
+ ,102.9
+ ,271244
+ ,4.1
+ ,7.4
+ ,97.4
+ ,269907
+ ,4
+ ,8.8
+ ,111.4
+ ,271296
+ ,3.8
+ ,9.3
+ ,87.4
+ ,270157
+ ,4.7
+ ,9.3
+ ,96.8
+ ,271322
+ ,4.3
+ ,8.7
+ ,114.1
+ ,267179
+ ,3.9
+ ,8.2
+ ,110.3
+ ,264101
+ ,4
+ ,8.3
+ ,103.9
+ ,265518
+ ,4.3
+ ,8.5
+ ,101.6
+ ,269419
+ ,4.8
+ ,8.6
+ ,94.6
+ ,268714
+ ,4.4
+ ,8.5
+ ,95.9
+ ,272482
+ ,4.3
+ ,8.2
+ ,104.7
+ ,268351
+ ,4.7
+ ,8.1
+ ,102.8
+ ,268175
+ ,4.7
+ ,7.9
+ ,98.1
+ ,270674
+ ,4.9
+ ,8.6
+ ,113.9
+ ,272764
+ ,5
+ ,8.7
+ ,80.9
+ ,272599
+ ,4.2
+ ,8.7
+ ,95.7
+ ,270333
+ ,4.3
+ ,8.5
+ ,113.2
+ ,270846
+ ,4.8
+ ,8.4
+ ,105.9
+ ,270491
+ ,4.8
+ ,8.5
+ ,108.8
+ ,269160
+ ,4.8
+ ,8.7
+ ,102.3
+ ,274027
+ ,4.2
+ ,8.7
+ ,99
+ ,273784
+ ,4.6
+ ,8.6
+ ,100.7
+ ,276663
+ ,4.8
+ ,8.5
+ ,115.5
+ ,274525
+ ,4.5
+ ,8.3
+ ,100.7
+ ,271344
+ ,4.4
+ ,8
+ ,109.9
+ ,271115
+ ,4.3
+ ,8.2
+ ,114.6
+ ,270798
+ ,3.9
+ ,8.1
+ ,85.4
+ ,273911
+ ,3.7
+ ,8.1
+ ,100.5
+ ,273985
+ ,4
+ ,8
+ ,114.8
+ ,271917
+ ,4.1
+ ,7.9
+ ,116.5
+ ,273338
+ ,3.7
+ ,7.9
+ ,112.9
+ ,270601
+ ,3.8
+ ,8
+ ,102
+ ,273547
+ ,3.8
+ ,8
+ ,106
+ ,275363
+ ,3.8
+ ,7.9
+ ,105.3
+ ,281229
+ ,3.3
+ ,8
+ ,118.8
+ ,277793
+ ,3.3
+ ,7.7
+ ,106.1
+ ,279913
+ ,3.3
+ ,7.2
+ ,109.3
+ ,282500
+ ,3.2
+ ,7.5
+ ,117.2
+ ,280041
+ ,3.4
+ ,7.3
+ ,92.5
+ ,282166
+ ,4.2
+ ,7
+ ,104.2
+ ,290304
+ ,4.9
+ ,7
+ ,112.5
+ ,283519
+ ,5.1
+ ,7
+ ,122.4
+ ,287816
+ ,5.5
+ ,7.2
+ ,113.3
+ ,285226
+ ,5.6
+ ,7.3
+ ,100
+ ,287595
+ ,6.4
+ ,7.1
+ ,110.7
+ ,289741
+ ,6.1
+ ,6.8
+ ,112.8
+ ,289148
+ ,7.1
+ ,6.4
+ ,109.8
+ ,288301
+ ,7.8
+ ,6.1
+ ,117.3
+ ,290155
+ ,7.9
+ ,6.5
+ ,109.1
+ ,289648
+ ,7.4
+ ,7.7
+ ,115.9
+ ,288225
+ ,7.5
+ ,7.9
+ ,96
+ ,289351
+ ,6.8
+ ,7.5
+ ,99.8
+ ,294735
+ ,5.2
+ ,6.9
+ ,116.8
+ ,305333
+ ,4.7
+ ,6.6
+ ,115.7
+ ,309030
+ ,4.1
+ ,6.9
+ ,99.4
+ ,310215
+ ,3.9
+ ,7.7
+ ,94.3
+ ,321935
+ ,2.6
+ ,8
+ ,91
+ ,325734
+ ,2.7
+ ,8
+ ,93.2
+ ,320846
+ ,1.8
+ ,7.7
+ ,103.1
+ ,323023
+ ,1
+ ,7.3
+ ,94.1
+ ,319753
+ ,0.3
+ ,7.4
+ ,91.8
+ ,321753
+ ,1.3
+ ,8.1
+ ,102.7
+ ,320757
+ ,1
+ ,8.3
+ ,82.6
+ ,324479
+ ,1.1
+ ,8.2
+ ,89.1
+ ,324641)
+ ,dim=c(4
+ ,65)
+ ,dimnames=list(c('Cons.index'
+ ,'Werkl.graad'
+ ,'Industr.prod.'
+ ,'BrutoSchuld')
+ ,1:65))
> y <- array(NA,dim=c(4,65),dimnames=list(c('Cons.index','Werkl.graad','Industr.prod.','BrutoSchuld'),1:65))
> 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
Cons.index Werkl.graad Industr.prod. BrutoSchuld M1 M2 M3 M4 M5 M6 M7 M8 M9
1 4.0 7.2 102.9 271244 1 0 0 0 0 0 0 0 0
2 4.1 7.4 97.4 269907 0 1 0 0 0 0 0 0 0
3 4.0 8.8 111.4 271296 0 0 1 0 0 0 0 0 0
4 3.8 9.3 87.4 270157 0 0 0 1 0 0 0 0 0
5 4.7 9.3 96.8 271322 0 0 0 0 1 0 0 0 0
6 4.3 8.7 114.1 267179 0 0 0 0 0 1 0 0 0
7 3.9 8.2 110.3 264101 0 0 0 0 0 0 1 0 0
8 4.0 8.3 103.9 265518 0 0 0 0 0 0 0 1 0
9 4.3 8.5 101.6 269419 0 0 0 0 0 0 0 0 1
10 4.8 8.6 94.6 268714 0 0 0 0 0 0 0 0 0
11 4.4 8.5 95.9 272482 0 0 0 0 0 0 0 0 0
12 4.3 8.2 104.7 268351 0 0 0 0 0 0 0 0 0
13 4.7 8.1 102.8 268175 1 0 0 0 0 0 0 0 0
14 4.7 7.9 98.1 270674 0 1 0 0 0 0 0 0 0
15 4.9 8.6 113.9 272764 0 0 1 0 0 0 0 0 0
16 5.0 8.7 80.9 272599 0 0 0 1 0 0 0 0 0
17 4.2 8.7 95.7 270333 0 0 0 0 1 0 0 0 0
18 4.3 8.5 113.2 270846 0 0 0 0 0 1 0 0 0
19 4.8 8.4 105.9 270491 0 0 0 0 0 0 1 0 0
20 4.8 8.5 108.8 269160 0 0 0 0 0 0 0 1 0
21 4.8 8.7 102.3 274027 0 0 0 0 0 0 0 0 1
22 4.2 8.7 99.0 273784 0 0 0 0 0 0 0 0 0
23 4.6 8.6 100.7 276663 0 0 0 0 0 0 0 0 0
24 4.8 8.5 115.5 274525 0 0 0 0 0 0 0 0 0
25 4.5 8.3 100.7 271344 1 0 0 0 0 0 0 0 0
26 4.4 8.0 109.9 271115 0 1 0 0 0 0 0 0 0
27 4.3 8.2 114.6 270798 0 0 1 0 0 0 0 0 0
28 3.9 8.1 85.4 273911 0 0 0 1 0 0 0 0 0
29 3.7 8.1 100.5 273985 0 0 0 0 1 0 0 0 0
30 4.0 8.0 114.8 271917 0 0 0 0 0 1 0 0 0
31 4.1 7.9 116.5 273338 0 0 0 0 0 0 1 0 0
32 3.7 7.9 112.9 270601 0 0 0 0 0 0 0 1 0
33 3.8 8.0 102.0 273547 0 0 0 0 0 0 0 0 1
34 3.8 8.0 106.0 275363 0 0 0 0 0 0 0 0 0
35 3.8 7.9 105.3 281229 0 0 0 0 0 0 0 0 0
36 3.3 8.0 118.8 277793 0 0 0 0 0 0 0 0 0
37 3.3 7.7 106.1 279913 1 0 0 0 0 0 0 0 0
38 3.3 7.2 109.3 282500 0 1 0 0 0 0 0 0 0
39 3.2 7.5 117.2 280041 0 0 1 0 0 0 0 0 0
40 3.4 7.3 92.5 282166 0 0 0 1 0 0 0 0 0
41 4.2 7.0 104.2 290304 0 0 0 0 1 0 0 0 0
42 4.9 7.0 112.5 283519 0 0 0 0 0 1 0 0 0
43 5.1 7.0 122.4 287816 0 0 0 0 0 0 1 0 0
44 5.5 7.2 113.3 285226 0 0 0 0 0 0 0 1 0
45 5.6 7.3 100.0 287595 0 0 0 0 0 0 0 0 1
46 6.4 7.1 110.7 289741 0 0 0 0 0 0 0 0 0
47 6.1 6.8 112.8 289148 0 0 0 0 0 0 0 0 0
48 7.1 6.4 109.8 288301 0 0 0 0 0 0 0 0 0
49 7.8 6.1 117.3 290155 1 0 0 0 0 0 0 0 0
50 7.9 6.5 109.1 289648 0 1 0 0 0 0 0 0 0
51 7.4 7.7 115.9 288225 0 0 1 0 0 0 0 0 0
52 7.5 7.9 96.0 289351 0 0 0 1 0 0 0 0 0
53 6.8 7.5 99.8 294735 0 0 0 0 1 0 0 0 0
54 5.2 6.9 116.8 305333 0 0 0 0 0 1 0 0 0
55 4.7 6.6 115.7 309030 0 0 0 0 0 0 1 0 0
56 4.1 6.9 99.4 310215 0 0 0 0 0 0 0 1 0
57 3.9 7.7 94.3 321935 0 0 0 0 0 0 0 0 1
58 2.6 8.0 91.0 325734 0 0 0 0 0 0 0 0 0
59 2.7 8.0 93.2 320846 0 0 0 0 0 0 0 0 0
60 1.8 7.7 103.1 323023 0 0 0 0 0 0 0 0 0
61 1.0 7.3 94.1 319753 1 0 0 0 0 0 0 0 0
62 0.3 7.4 91.8 321753 0 1 0 0 0 0 0 0 0
63 1.3 8.1 102.7 320757 0 0 1 0 0 0 0 0 0
64 1.0 8.3 82.6 324479 0 0 0 1 0 0 0 0 0
65 1.1 8.2 89.1 324641 0 0 0 0 1 0 0 0 0
M10 M11 t
1 0 0 1
2 0 0 2
3 0 0 3
4 0 0 4
5 0 0 5
6 0 0 6
7 0 0 7
8 0 0 8
9 0 0 9
10 1 0 10
11 0 1 11
12 0 0 12
13 0 0 13
14 0 0 14
15 0 0 15
16 0 0 16
17 0 0 17
18 0 0 18
19 0 0 19
20 0 0 20
21 0 0 21
22 1 0 22
23 0 1 23
24 0 0 24
25 0 0 25
26 0 0 26
27 0 0 27
28 0 0 28
29 0 0 29
30 0 0 30
31 0 0 31
32 0 0 32
33 0 0 33
34 1 0 34
35 0 1 35
36 0 0 36
37 0 0 37
38 0 0 38
39 0 0 39
40 0 0 40
41 0 0 41
42 0 0 42
43 0 0 43
44 0 0 44
45 0 0 45
46 1 0 46
47 0 1 47
48 0 0 48
49 0 0 49
50 0 0 50
51 0 0 51
52 0 0 52
53 0 0 53
54 0 0 54
55 0 0 55
56 0 0 56
57 0 0 57
58 1 0 58
59 0 1 59
60 0 0 60
61 0 0 61
62 0 0 62
63 0 0 63
64 0 0 64
65 0 0 65
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Werkl.graad Industr.prod. BrutoSchuld M1
19.8323313 -1.1152724 0.0592692 -0.0000473 -0.1383295
M2 M3 M4 M5 M6
-0.1749844 0.1185306 1.7227645 1.0825725 -0.1838970
M7 M8 M9 M10 M11
-0.3656494 0.0349714 1.1005706 1.0741931 0.8862475
t
0.0023985
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.2458 -0.7513 0.0713 0.5752 2.6778
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.983e+01 1.155e+01 1.717 0.09230 .
Werkl.graad -1.115e+00 3.818e-01 -2.921 0.00526 **
Industr.prod. 5.927e-02 4.285e-02 1.383 0.17288
BrutoSchuld -4.730e-05 2.633e-05 -1.796 0.07862 .
M1 -1.383e-01 8.362e-01 -0.165 0.86928
M2 -1.750e-01 8.631e-01 -0.203 0.84018
M3 1.185e-01 7.802e-01 0.152 0.87987
M4 1.723e+00 1.191e+00 1.446 0.15441
M5 1.083e+00 8.943e-01 1.211 0.23187
M6 -1.839e-01 8.046e-01 -0.229 0.82016
M7 -3.656e-01 8.045e-01 -0.455 0.65147
M8 3.497e-02 8.146e-01 0.043 0.96593
M9 1.101e+00 8.849e-01 1.244 0.21952
M10 1.074e+00 8.750e-01 1.228 0.22545
M11 8.862e-01 8.520e-01 1.040 0.30336
t 2.398e-03 2.539e-02 0.094 0.92512
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.252 on 49 degrees of freedom
Multiple R-squared: 0.4915, Adjusted R-squared: 0.3358
F-statistic: 3.157 on 15 and 49 DF, p-value: 0.001195
> 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.430276e-02 1.086055e-01 0.9456972
[2,] 1.792860e-02 3.585721e-02 0.9820714
[3,] 5.457346e-03 1.091469e-02 0.9945427
[4,] 4.635666e-03 9.271333e-03 0.9953643
[5,] 1.689128e-03 3.378257e-03 0.9983109
[6,] 9.145533e-04 1.829107e-03 0.9990854
[7,] 1.045998e-03 2.091996e-03 0.9989540
[8,] 4.830121e-04 9.660242e-04 0.9995170
[9,] 2.446219e-04 4.892439e-04 0.9997554
[10,] 2.698341e-04 5.396682e-04 0.9997302
[11,] 2.032850e-04 4.065700e-04 0.9997967
[12,] 9.895553e-05 1.979111e-04 0.9999010
[13,] 5.687242e-05 1.137448e-04 0.9999431
[14,] 2.396101e-05 4.792202e-05 0.9999760
[15,] 1.639731e-05 3.279462e-05 0.9999836
[16,] 6.992559e-06 1.398512e-05 0.9999930
[17,] 3.338219e-06 6.676437e-06 0.9999967
[18,] 6.361221e-06 1.272244e-05 0.9999936
[19,] 1.720726e-05 3.441451e-05 0.9999828
[20,] 6.433230e-06 1.286646e-05 0.9999936
[21,] 3.123672e-06 6.247343e-06 0.9999969
[22,] 2.061492e-06 4.122983e-06 0.9999979
[23,] 4.551514e-06 9.103028e-06 0.9999954
[24,] 3.874918e-06 7.749836e-06 0.9999961
[25,] 6.442635e-06 1.288527e-05 0.9999936
[26,] 6.900617e-05 1.380123e-04 0.9999310
[27,] 3.575987e-03 7.151975e-03 0.9964240
[28,] 1.586332e-02 3.172664e-02 0.9841367
> postscript(file="/var/www/html/rcomp/tmp/1k6d41258650132.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/2zbbs1258650132.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/3bsrx1258650132.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/41u371258650132.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/537mi1258650132.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 = 65
Frequency = 1
1 2 3 4 5 6
-0.93573600 -0.31568284 0.08571392 0.20530582 1.24107165 0.21466295
7 8 9 10 11 12
-0.48398187 -0.22912858 -0.45324015 0.56380504 0.33899663 0.07130363
13 14 15 16 17 18
0.59999429 0.80796117 0.65513909 1.20811401 0.06154400 0.18961345
19 20 21 22 23 24
1.17331439 0.64698693 0.41749669 0.02557056 0.53500586 0.52901845
25 26 27 28 29 30
0.86862203 -0.08781152 -0.55422958 -0.79448683 -1.24815850 -0.74097841
31 32 33 34 35 36
-0.60669785 -1.32580468 -1.39689851 -1.52410180 -1.13113896 -1.59841570
37 38 39 40 41 42
-0.94407398 -1.53475337 -2.18061964 -2.24584680 -1.45116761 -0.29995323
43 44 45 46 47 48
-0.30412179 0.33275949 0.37662040 0.44486744 -0.15668073 1.41880497
49 50 51 52 53 54
1.56332697 2.60571948 2.67779602 2.62693410 2.14805209 0.63665523
55 56 57 58 59 60
0.22148712 0.57518685 1.05602157 0.48985876 0.41381720 -0.42071135
61 62 63 64 65
-1.15213332 -1.47543293 -0.68379981 -1.00002030 -0.75134163
> postscript(file="/var/www/html/rcomp/tmp/6emag1258650132.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 = 65
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.93573600 NA
1 -0.31568284 -0.93573600
2 0.08571392 -0.31568284
3 0.20530582 0.08571392
4 1.24107165 0.20530582
5 0.21466295 1.24107165
6 -0.48398187 0.21466295
7 -0.22912858 -0.48398187
8 -0.45324015 -0.22912858
9 0.56380504 -0.45324015
10 0.33899663 0.56380504
11 0.07130363 0.33899663
12 0.59999429 0.07130363
13 0.80796117 0.59999429
14 0.65513909 0.80796117
15 1.20811401 0.65513909
16 0.06154400 1.20811401
17 0.18961345 0.06154400
18 1.17331439 0.18961345
19 0.64698693 1.17331439
20 0.41749669 0.64698693
21 0.02557056 0.41749669
22 0.53500586 0.02557056
23 0.52901845 0.53500586
24 0.86862203 0.52901845
25 -0.08781152 0.86862203
26 -0.55422958 -0.08781152
27 -0.79448683 -0.55422958
28 -1.24815850 -0.79448683
29 -0.74097841 -1.24815850
30 -0.60669785 -0.74097841
31 -1.32580468 -0.60669785
32 -1.39689851 -1.32580468
33 -1.52410180 -1.39689851
34 -1.13113896 -1.52410180
35 -1.59841570 -1.13113896
36 -0.94407398 -1.59841570
37 -1.53475337 -0.94407398
38 -2.18061964 -1.53475337
39 -2.24584680 -2.18061964
40 -1.45116761 -2.24584680
41 -0.29995323 -1.45116761
42 -0.30412179 -0.29995323
43 0.33275949 -0.30412179
44 0.37662040 0.33275949
45 0.44486744 0.37662040
46 -0.15668073 0.44486744
47 1.41880497 -0.15668073
48 1.56332697 1.41880497
49 2.60571948 1.56332697
50 2.67779602 2.60571948
51 2.62693410 2.67779602
52 2.14805209 2.62693410
53 0.63665523 2.14805209
54 0.22148712 0.63665523
55 0.57518685 0.22148712
56 1.05602157 0.57518685
57 0.48985876 1.05602157
58 0.41381720 0.48985876
59 -0.42071135 0.41381720
60 -1.15213332 -0.42071135
61 -1.47543293 -1.15213332
62 -0.68379981 -1.47543293
63 -1.00002030 -0.68379981
64 -0.75134163 -1.00002030
65 NA -0.75134163
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.31568284 -0.93573600
[2,] 0.08571392 -0.31568284
[3,] 0.20530582 0.08571392
[4,] 1.24107165 0.20530582
[5,] 0.21466295 1.24107165
[6,] -0.48398187 0.21466295
[7,] -0.22912858 -0.48398187
[8,] -0.45324015 -0.22912858
[9,] 0.56380504 -0.45324015
[10,] 0.33899663 0.56380504
[11,] 0.07130363 0.33899663
[12,] 0.59999429 0.07130363
[13,] 0.80796117 0.59999429
[14,] 0.65513909 0.80796117
[15,] 1.20811401 0.65513909
[16,] 0.06154400 1.20811401
[17,] 0.18961345 0.06154400
[18,] 1.17331439 0.18961345
[19,] 0.64698693 1.17331439
[20,] 0.41749669 0.64698693
[21,] 0.02557056 0.41749669
[22,] 0.53500586 0.02557056
[23,] 0.52901845 0.53500586
[24,] 0.86862203 0.52901845
[25,] -0.08781152 0.86862203
[26,] -0.55422958 -0.08781152
[27,] -0.79448683 -0.55422958
[28,] -1.24815850 -0.79448683
[29,] -0.74097841 -1.24815850
[30,] -0.60669785 -0.74097841
[31,] -1.32580468 -0.60669785
[32,] -1.39689851 -1.32580468
[33,] -1.52410180 -1.39689851
[34,] -1.13113896 -1.52410180
[35,] -1.59841570 -1.13113896
[36,] -0.94407398 -1.59841570
[37,] -1.53475337 -0.94407398
[38,] -2.18061964 -1.53475337
[39,] -2.24584680 -2.18061964
[40,] -1.45116761 -2.24584680
[41,] -0.29995323 -1.45116761
[42,] -0.30412179 -0.29995323
[43,] 0.33275949 -0.30412179
[44,] 0.37662040 0.33275949
[45,] 0.44486744 0.37662040
[46,] -0.15668073 0.44486744
[47,] 1.41880497 -0.15668073
[48,] 1.56332697 1.41880497
[49,] 2.60571948 1.56332697
[50,] 2.67779602 2.60571948
[51,] 2.62693410 2.67779602
[52,] 2.14805209 2.62693410
[53,] 0.63665523 2.14805209
[54,] 0.22148712 0.63665523
[55,] 0.57518685 0.22148712
[56,] 1.05602157 0.57518685
[57,] 0.48985876 1.05602157
[58,] 0.41381720 0.48985876
[59,] -0.42071135 0.41381720
[60,] -1.15213332 -0.42071135
[61,] -1.47543293 -1.15213332
[62,] -0.68379981 -1.47543293
[63,] -1.00002030 -0.68379981
[64,] -0.75134163 -1.00002030
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.31568284 -0.93573600
2 0.08571392 -0.31568284
3 0.20530582 0.08571392
4 1.24107165 0.20530582
5 0.21466295 1.24107165
6 -0.48398187 0.21466295
7 -0.22912858 -0.48398187
8 -0.45324015 -0.22912858
9 0.56380504 -0.45324015
10 0.33899663 0.56380504
11 0.07130363 0.33899663
12 0.59999429 0.07130363
13 0.80796117 0.59999429
14 0.65513909 0.80796117
15 1.20811401 0.65513909
16 0.06154400 1.20811401
17 0.18961345 0.06154400
18 1.17331439 0.18961345
19 0.64698693 1.17331439
20 0.41749669 0.64698693
21 0.02557056 0.41749669
22 0.53500586 0.02557056
23 0.52901845 0.53500586
24 0.86862203 0.52901845
25 -0.08781152 0.86862203
26 -0.55422958 -0.08781152
27 -0.79448683 -0.55422958
28 -1.24815850 -0.79448683
29 -0.74097841 -1.24815850
30 -0.60669785 -0.74097841
31 -1.32580468 -0.60669785
32 -1.39689851 -1.32580468
33 -1.52410180 -1.39689851
34 -1.13113896 -1.52410180
35 -1.59841570 -1.13113896
36 -0.94407398 -1.59841570
37 -1.53475337 -0.94407398
38 -2.18061964 -1.53475337
39 -2.24584680 -2.18061964
40 -1.45116761 -2.24584680
41 -0.29995323 -1.45116761
42 -0.30412179 -0.29995323
43 0.33275949 -0.30412179
44 0.37662040 0.33275949
45 0.44486744 0.37662040
46 -0.15668073 0.44486744
47 1.41880497 -0.15668073
48 1.56332697 1.41880497
49 2.60571948 1.56332697
50 2.67779602 2.60571948
51 2.62693410 2.67779602
52 2.14805209 2.62693410
53 0.63665523 2.14805209
54 0.22148712 0.63665523
55 0.57518685 0.22148712
56 1.05602157 0.57518685
57 0.48985876 1.05602157
58 0.41381720 0.48985876
59 -0.42071135 0.41381720
60 -1.15213332 -0.42071135
61 -1.47543293 -1.15213332
62 -0.68379981 -1.47543293
63 -1.00002030 -0.68379981
64 -0.75134163 -1.00002030
> 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/78bm11258650132.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/8dle21258650132.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/9d9hx1258650132.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/10e5qq1258650132.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/11bl691258650132.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/12k7pa1258650132.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/136n4z1258650132.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/14l3221258650132.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/15c2d01258650132.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/1617ka1258650133.tab")
+ }
>
> system("convert tmp/1k6d41258650132.ps tmp/1k6d41258650132.png")
> system("convert tmp/2zbbs1258650132.ps tmp/2zbbs1258650132.png")
> system("convert tmp/3bsrx1258650132.ps tmp/3bsrx1258650132.png")
> system("convert tmp/41u371258650132.ps tmp/41u371258650132.png")
> system("convert tmp/537mi1258650132.ps tmp/537mi1258650132.png")
> system("convert tmp/6emag1258650132.ps tmp/6emag1258650132.png")
> system("convert tmp/78bm11258650132.ps tmp/78bm11258650132.png")
> system("convert tmp/8dle21258650132.ps tmp/8dle21258650132.png")
> system("convert tmp/9d9hx1258650132.ps tmp/9d9hx1258650132.png")
> system("convert tmp/10e5qq1258650132.ps tmp/10e5qq1258650132.png")
>
>
> proc.time()
user system elapsed
2.415 1.550 2.942