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 = 'Do not include Seasonal 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 t
1 4.0 7.2 102.9 271244 1
2 4.1 7.4 97.4 269907 2
3 4.0 8.8 111.4 271296 3
4 3.8 9.3 87.4 270157 4
5 4.7 9.3 96.8 271322 5
6 4.3 8.7 114.1 267179 6
7 3.9 8.2 110.3 264101 7
8 4.0 8.3 103.9 265518 8
9 4.3 8.5 101.6 269419 9
10 4.8 8.6 94.6 268714 10
11 4.4 8.5 95.9 272482 11
12 4.3 8.2 104.7 268351 12
13 4.7 8.1 102.8 268175 13
14 4.7 7.9 98.1 270674 14
15 4.9 8.6 113.9 272764 15
16 5.0 8.7 80.9 272599 16
17 4.2 8.7 95.7 270333 17
18 4.3 8.5 113.2 270846 18
19 4.8 8.4 105.9 270491 19
20 4.8 8.5 108.8 269160 20
21 4.8 8.7 102.3 274027 21
22 4.2 8.7 99.0 273784 22
23 4.6 8.6 100.7 276663 23
24 4.8 8.5 115.5 274525 24
25 4.5 8.3 100.7 271344 25
26 4.4 8.0 109.9 271115 26
27 4.3 8.2 114.6 270798 27
28 3.9 8.1 85.4 273911 28
29 3.7 8.1 100.5 273985 29
30 4.0 8.0 114.8 271917 30
31 4.1 7.9 116.5 273338 31
32 3.7 7.9 112.9 270601 32
33 3.8 8.0 102.0 273547 33
34 3.8 8.0 106.0 275363 34
35 3.8 7.9 105.3 281229 35
36 3.3 8.0 118.8 277793 36
37 3.3 7.7 106.1 279913 37
38 3.3 7.2 109.3 282500 38
39 3.2 7.5 117.2 280041 39
40 3.4 7.3 92.5 282166 40
41 4.2 7.0 104.2 290304 41
42 4.9 7.0 112.5 283519 42
43 5.1 7.0 122.4 287816 43
44 5.5 7.2 113.3 285226 44
45 5.6 7.3 100.0 287595 45
46 6.4 7.1 110.7 289741 46
47 6.1 6.8 112.8 289148 47
48 7.1 6.4 109.8 288301 48
49 7.8 6.1 117.3 290155 49
50 7.9 6.5 109.1 289648 50
51 7.4 7.7 115.9 288225 51
52 7.5 7.9 96.0 289351 52
53 6.8 7.5 99.8 294735 53
54 5.2 6.9 116.8 305333 54
55 4.7 6.6 115.7 309030 55
56 4.1 6.9 99.4 310215 56
57 3.9 7.7 94.3 321935 57
58 2.6 8.0 91.0 325734 58
59 2.7 8.0 93.2 320846 59
60 1.8 7.7 103.1 323023 60
61 1.0 7.3 94.1 319753 61
62 0.3 7.4 91.8 321753 62
63 1.3 8.1 102.7 320757 63
64 1.0 8.3 82.6 324479 64
65 1.1 8.2 89.1 324641 65
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Werkl.graad Industr.prod. BrutoSchuld t
3.105e+01 -8.557e-01 4.182e-03 -7.585e-05 3.361e-02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.48193 -0.79897 0.05863 0.69848 3.00682
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.105e+01 6.889e+00 4.507 3.1e-05 ***
Werkl.graad -8.557e-01 3.043e-01 -2.812 0.006647 **
Industr.prod. 4.182e-03 1.955e-02 0.214 0.831299
BrutoSchuld -7.585e-05 2.023e-05 -3.748 0.000402 ***
t 3.361e-02 2.078e-02 1.617 0.111032
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.215 on 60 degrees of freedom
Multiple R-squared: 0.4138, Adjusted R-squared: 0.3747
F-statistic: 10.59 on 4 and 60 DF, p-value: 1.478e-06
> 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,] 1.720470e-02 3.440939e-02 0.98279530
[2,] 6.276862e-03 1.255372e-02 0.99372314
[3,] 2.102922e-03 4.205844e-03 0.99789708
[4,] 1.107748e-03 2.215496e-03 0.99889225
[5,] 2.790702e-04 5.581404e-04 0.99972093
[6,] 7.302160e-05 1.460432e-04 0.99992698
[7,] 1.473948e-05 2.947896e-05 0.99998526
[8,] 3.093253e-06 6.186506e-06 0.99999691
[9,] 7.073745e-07 1.414749e-06 0.99999929
[10,] 1.582663e-06 3.165326e-06 0.99999842
[11,] 9.085087e-07 1.817017e-06 0.99999909
[12,] 2.275128e-07 4.550256e-07 0.99999977
[13,] 5.607007e-08 1.121401e-07 0.99999994
[14,] 2.008247e-08 4.016494e-08 0.99999998
[15,] 6.388913e-08 1.277783e-07 0.99999994
[16,] 4.391286e-08 8.782571e-08 0.99999996
[17,] 2.640518e-08 5.281037e-08 0.99999997
[18,] 1.487164e-08 2.974328e-08 0.99999999
[19,] 7.117190e-09 1.423438e-08 0.99999999
[20,] 3.938714e-09 7.877429e-09 1.00000000
[21,] 1.609311e-08 3.218622e-08 0.99999998
[22,] 5.697236e-08 1.139447e-07 0.99999994
[23,] 3.047884e-08 6.095767e-08 0.99999997
[24,] 1.194607e-08 2.389214e-08 0.99999999
[25,] 8.621920e-09 1.724384e-08 0.99999999
[26,] 4.308678e-09 8.617356e-09 1.00000000
[27,] 2.081744e-09 4.163488e-09 1.00000000
[28,] 1.239117e-09 2.478235e-09 1.00000000
[29,] 1.865505e-09 3.731010e-09 1.00000000
[30,] 1.278780e-09 2.557560e-09 1.00000000
[31,] 6.756047e-10 1.351209e-09 1.00000000
[32,] 2.595417e-09 5.190834e-09 1.00000000
[33,] 4.790137e-09 9.580274e-09 1.00000000
[34,] 1.461468e-08 2.922935e-08 0.99999999
[35,] 7.898309e-07 1.579662e-06 0.99999921
[36,] 1.044726e-05 2.089451e-05 0.99998955
[37,] 4.488990e-04 8.977980e-04 0.99955110
[38,] 8.910460e-03 1.782092e-02 0.99108954
[39,] 8.843259e-02 1.768652e-01 0.91156741
[40,] 6.795518e-01 6.408965e-01 0.32044824
[41,] 9.201776e-01 1.596449e-01 0.07982243
[42,] 9.377078e-01 1.245845e-01 0.06229225
[43,] 9.520915e-01 9.581701e-02 0.04790851
[44,] 9.663499e-01 6.730014e-02 0.03365007
[45,] 9.615084e-01 7.698320e-02 0.03849160
[46,] 9.471806e-01 1.056388e-01 0.05281938
[47,] 9.385491e-01 1.229017e-01 0.06145087
[48,] 9.294663e-01 1.410674e-01 0.07053371
[49,] 8.913552e-01 2.172897e-01 0.10864483
[50,] 9.604980e-01 7.900393e-02 0.03950197
> postscript(file="/var/www/html/rcomp/tmp/132bn1258566447.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/273m21258566447.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/3endx1258566447.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/4uia71258566447.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/52ujw1258566447.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.78024571 -0.62112828 0.48999863 0.69820808 1.61364631 0.28004142
7 8 9 10 11 12
-0.79897114 -0.51277134 0.23025414 0.75801401 0.51919440 -0.22124956
13 14 15 16 17 18
0.05417003 0.05862644 0.91642581 1.19388433 0.12650365 -0.01252183
19 20 21 22 23 24
0.37190560 0.31077952 0.84463992 0.20639993 0.69847857 0.55523947
25 26 27 28 29 30
-0.12887883 -0.57503706 -0.58121456 -0.74215329 -1.03330449 -1.06914307
31 32 33 34 35 36
-0.98765073 -1.61380079 -1.19280897 -1.10540910 -0.77673505 -1.54185382
37 38 39 40 41 42
-1.61825204 -1.89686236 -1.99332268 -1.73358666 -0.65558098 -0.53853332
43 44 45 46 47 48
-0.08763050 0.29150587 0.67877104 1.39204514 0.74797286 1.32039815
49 50 51 52 53 54
1.83934117 2.24383863 2.60066002 3.00681659 2.32341053 0.90913486
55 56 57 58 59 60
0.40383374 0.18497588 1.54616875 0.77120704 0.45765050 -0.60894521
61 62 63 64 65
-1.99520546 -2.48193389 -1.03770782 -0.83381352 -0.86788896
> postscript(file="/var/www/html/rcomp/tmp/637dy1258566447.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.78024571 NA
1 -0.62112828 -0.78024571
2 0.48999863 -0.62112828
3 0.69820808 0.48999863
4 1.61364631 0.69820808
5 0.28004142 1.61364631
6 -0.79897114 0.28004142
7 -0.51277134 -0.79897114
8 0.23025414 -0.51277134
9 0.75801401 0.23025414
10 0.51919440 0.75801401
11 -0.22124956 0.51919440
12 0.05417003 -0.22124956
13 0.05862644 0.05417003
14 0.91642581 0.05862644
15 1.19388433 0.91642581
16 0.12650365 1.19388433
17 -0.01252183 0.12650365
18 0.37190560 -0.01252183
19 0.31077952 0.37190560
20 0.84463992 0.31077952
21 0.20639993 0.84463992
22 0.69847857 0.20639993
23 0.55523947 0.69847857
24 -0.12887883 0.55523947
25 -0.57503706 -0.12887883
26 -0.58121456 -0.57503706
27 -0.74215329 -0.58121456
28 -1.03330449 -0.74215329
29 -1.06914307 -1.03330449
30 -0.98765073 -1.06914307
31 -1.61380079 -0.98765073
32 -1.19280897 -1.61380079
33 -1.10540910 -1.19280897
34 -0.77673505 -1.10540910
35 -1.54185382 -0.77673505
36 -1.61825204 -1.54185382
37 -1.89686236 -1.61825204
38 -1.99332268 -1.89686236
39 -1.73358666 -1.99332268
40 -0.65558098 -1.73358666
41 -0.53853332 -0.65558098
42 -0.08763050 -0.53853332
43 0.29150587 -0.08763050
44 0.67877104 0.29150587
45 1.39204514 0.67877104
46 0.74797286 1.39204514
47 1.32039815 0.74797286
48 1.83934117 1.32039815
49 2.24383863 1.83934117
50 2.60066002 2.24383863
51 3.00681659 2.60066002
52 2.32341053 3.00681659
53 0.90913486 2.32341053
54 0.40383374 0.90913486
55 0.18497588 0.40383374
56 1.54616875 0.18497588
57 0.77120704 1.54616875
58 0.45765050 0.77120704
59 -0.60894521 0.45765050
60 -1.99520546 -0.60894521
61 -2.48193389 -1.99520546
62 -1.03770782 -2.48193389
63 -0.83381352 -1.03770782
64 -0.86788896 -0.83381352
65 NA -0.86788896
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.62112828 -0.78024571
[2,] 0.48999863 -0.62112828
[3,] 0.69820808 0.48999863
[4,] 1.61364631 0.69820808
[5,] 0.28004142 1.61364631
[6,] -0.79897114 0.28004142
[7,] -0.51277134 -0.79897114
[8,] 0.23025414 -0.51277134
[9,] 0.75801401 0.23025414
[10,] 0.51919440 0.75801401
[11,] -0.22124956 0.51919440
[12,] 0.05417003 -0.22124956
[13,] 0.05862644 0.05417003
[14,] 0.91642581 0.05862644
[15,] 1.19388433 0.91642581
[16,] 0.12650365 1.19388433
[17,] -0.01252183 0.12650365
[18,] 0.37190560 -0.01252183
[19,] 0.31077952 0.37190560
[20,] 0.84463992 0.31077952
[21,] 0.20639993 0.84463992
[22,] 0.69847857 0.20639993
[23,] 0.55523947 0.69847857
[24,] -0.12887883 0.55523947
[25,] -0.57503706 -0.12887883
[26,] -0.58121456 -0.57503706
[27,] -0.74215329 -0.58121456
[28,] -1.03330449 -0.74215329
[29,] -1.06914307 -1.03330449
[30,] -0.98765073 -1.06914307
[31,] -1.61380079 -0.98765073
[32,] -1.19280897 -1.61380079
[33,] -1.10540910 -1.19280897
[34,] -0.77673505 -1.10540910
[35,] -1.54185382 -0.77673505
[36,] -1.61825204 -1.54185382
[37,] -1.89686236 -1.61825204
[38,] -1.99332268 -1.89686236
[39,] -1.73358666 -1.99332268
[40,] -0.65558098 -1.73358666
[41,] -0.53853332 -0.65558098
[42,] -0.08763050 -0.53853332
[43,] 0.29150587 -0.08763050
[44,] 0.67877104 0.29150587
[45,] 1.39204514 0.67877104
[46,] 0.74797286 1.39204514
[47,] 1.32039815 0.74797286
[48,] 1.83934117 1.32039815
[49,] 2.24383863 1.83934117
[50,] 2.60066002 2.24383863
[51,] 3.00681659 2.60066002
[52,] 2.32341053 3.00681659
[53,] 0.90913486 2.32341053
[54,] 0.40383374 0.90913486
[55,] 0.18497588 0.40383374
[56,] 1.54616875 0.18497588
[57,] 0.77120704 1.54616875
[58,] 0.45765050 0.77120704
[59,] -0.60894521 0.45765050
[60,] -1.99520546 -0.60894521
[61,] -2.48193389 -1.99520546
[62,] -1.03770782 -2.48193389
[63,] -0.83381352 -1.03770782
[64,] -0.86788896 -0.83381352
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.62112828 -0.78024571
2 0.48999863 -0.62112828
3 0.69820808 0.48999863
4 1.61364631 0.69820808
5 0.28004142 1.61364631
6 -0.79897114 0.28004142
7 -0.51277134 -0.79897114
8 0.23025414 -0.51277134
9 0.75801401 0.23025414
10 0.51919440 0.75801401
11 -0.22124956 0.51919440
12 0.05417003 -0.22124956
13 0.05862644 0.05417003
14 0.91642581 0.05862644
15 1.19388433 0.91642581
16 0.12650365 1.19388433
17 -0.01252183 0.12650365
18 0.37190560 -0.01252183
19 0.31077952 0.37190560
20 0.84463992 0.31077952
21 0.20639993 0.84463992
22 0.69847857 0.20639993
23 0.55523947 0.69847857
24 -0.12887883 0.55523947
25 -0.57503706 -0.12887883
26 -0.58121456 -0.57503706
27 -0.74215329 -0.58121456
28 -1.03330449 -0.74215329
29 -1.06914307 -1.03330449
30 -0.98765073 -1.06914307
31 -1.61380079 -0.98765073
32 -1.19280897 -1.61380079
33 -1.10540910 -1.19280897
34 -0.77673505 -1.10540910
35 -1.54185382 -0.77673505
36 -1.61825204 -1.54185382
37 -1.89686236 -1.61825204
38 -1.99332268 -1.89686236
39 -1.73358666 -1.99332268
40 -0.65558098 -1.73358666
41 -0.53853332 -0.65558098
42 -0.08763050 -0.53853332
43 0.29150587 -0.08763050
44 0.67877104 0.29150587
45 1.39204514 0.67877104
46 0.74797286 1.39204514
47 1.32039815 0.74797286
48 1.83934117 1.32039815
49 2.24383863 1.83934117
50 2.60066002 2.24383863
51 3.00681659 2.60066002
52 2.32341053 3.00681659
53 0.90913486 2.32341053
54 0.40383374 0.90913486
55 0.18497588 0.40383374
56 1.54616875 0.18497588
57 0.77120704 1.54616875
58 0.45765050 0.77120704
59 -0.60894521 0.45765050
60 -1.99520546 -0.60894521
61 -2.48193389 -1.99520546
62 -1.03770782 -2.48193389
63 -0.83381352 -1.03770782
64 -0.86788896 -0.83381352
> 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/7al5c1258566447.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/8fslg1258566447.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/9zfai1258566447.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/10xevx1258566447.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/11h7zg1258566447.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/12qd751258566447.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/13e3qw1258566447.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/14twx51258566447.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/15xxvg1258566447.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/16sdcy1258566447.tab")
+ }
>
> system("convert tmp/132bn1258566447.ps tmp/132bn1258566447.png")
> system("convert tmp/273m21258566447.ps tmp/273m21258566447.png")
> system("convert tmp/3endx1258566447.ps tmp/3endx1258566447.png")
> system("convert tmp/4uia71258566447.ps tmp/4uia71258566447.png")
> system("convert tmp/52ujw1258566447.ps tmp/52ujw1258566447.png")
> system("convert tmp/637dy1258566447.ps tmp/637dy1258566447.png")
> system("convert tmp/7al5c1258566447.ps tmp/7al5c1258566447.png")
> system("convert tmp/8fslg1258566447.ps tmp/8fslg1258566447.png")
> system("convert tmp/9zfai1258566447.ps tmp/9zfai1258566447.png")
> system("convert tmp/10xevx1258566447.ps tmp/10xevx1258566447.png")
>
>
> proc.time()
user system elapsed
2.500 1.586 2.967