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 = 'No 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
1 4.0 7.2 102.9 271244
2 4.1 7.4 97.4 269907
3 4.0 8.8 111.4 271296
4 3.8 9.3 87.4 270157
5 4.7 9.3 96.8 271322
6 4.3 8.7 114.1 267179
7 3.9 8.2 110.3 264101
8 4.0 8.3 103.9 265518
9 4.3 8.5 101.6 269419
10 4.8 8.6 94.6 268714
11 4.4 8.5 95.9 272482
12 4.3 8.2 104.7 268351
13 4.7 8.1 102.8 268175
14 4.7 7.9 98.1 270674
15 4.9 8.6 113.9 272764
16 5.0 8.7 80.9 272599
17 4.2 8.7 95.7 270333
18 4.3 8.5 113.2 270846
19 4.8 8.4 105.9 270491
20 4.8 8.5 108.8 269160
21 4.8 8.7 102.3 274027
22 4.2 8.7 99.0 273784
23 4.6 8.6 100.7 276663
24 4.8 8.5 115.5 274525
25 4.5 8.3 100.7 271344
26 4.4 8.0 109.9 271115
27 4.3 8.2 114.6 270798
28 3.9 8.1 85.4 273911
29 3.7 8.1 100.5 273985
30 4.0 8.0 114.8 271917
31 4.1 7.9 116.5 273338
32 3.7 7.9 112.9 270601
33 3.8 8.0 102.0 273547
34 3.8 8.0 106.0 275363
35 3.8 7.9 105.3 281229
36 3.3 8.0 118.8 277793
37 3.3 7.7 106.1 279913
38 3.3 7.2 109.3 282500
39 3.2 7.5 117.2 280041
40 3.4 7.3 92.5 282166
41 4.2 7.0 104.2 290304
42 4.9 7.0 112.5 283519
43 5.1 7.0 122.4 287816
44 5.5 7.2 113.3 285226
45 5.6 7.3 100.0 287595
46 6.4 7.1 110.7 289741
47 6.1 6.8 112.8 289148
48 7.1 6.4 109.8 288301
49 7.8 6.1 117.3 290155
50 7.9 6.5 109.1 289648
51 7.4 7.7 115.9 288225
52 7.5 7.9 96.0 289351
53 6.8 7.5 99.8 294735
54 5.2 6.9 116.8 305333
55 4.7 6.6 115.7 309030
56 4.1 6.9 99.4 310215
57 3.9 7.7 94.3 321935
58 2.6 8.0 91.0 325734
59 2.7 8.0 93.2 320846
60 1.8 7.7 103.1 323023
61 1.0 7.3 94.1 319753
62 0.3 7.4 91.8 321753
63 1.3 8.1 102.7 320757
64 1.0 8.3 82.6 324479
65 1.1 8.2 89.1 324641
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Werkl.graad Industr.prod. BrutoSchuld
2.506e+01 -1.068e+00 1.142e-02 -4.764e-05
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.57270 -0.80356 -0.03798 0.65090 3.57008
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.506e+01 5.884e+00 4.258 7.23e-05 ***
Werkl.graad -1.068e+00 2.780e-01 -3.843 0.000292 ***
Industr.prod. 1.142e-02 1.928e-02 0.592 0.555959
BrutoSchuld -4.764e-05 1.039e-05 -4.584 2.31e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.231 on 61 degrees of freedom
Multiple R-squared: 0.3882, Adjusted R-squared: 0.3581
F-statistic: 12.9 on 3 and 61 DF, p-value: 1.248e-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,] 3.820036e-02 7.640071e-02 0.9617996
[2,] 9.584512e-03 1.916902e-02 0.9904155
[3,] 2.474112e-03 4.948224e-03 0.9975259
[4,] 3.139327e-03 6.278654e-03 0.9968607
[5,] 8.634089e-04 1.726818e-03 0.9991366
[6,] 2.352979e-04 4.705958e-04 0.9997647
[7,] 1.620151e-04 3.240303e-04 0.9998380
[8,] 7.557450e-05 1.511490e-04 0.9999244
[9,] 3.822256e-05 7.644511e-05 0.9999618
[10,] 2.301333e-05 4.602666e-05 0.9999770
[11,] 7.524891e-06 1.504978e-05 0.9999925
[12,] 2.004641e-06 4.009281e-06 0.9999980
[13,] 9.118979e-07 1.823796e-06 0.9999991
[14,] 4.781220e-07 9.562441e-07 0.9999995
[15,] 1.490136e-07 2.980271e-07 0.9999999
[16,] 8.399317e-08 1.679863e-07 0.9999999
[17,] 2.794349e-08 5.588698e-08 1.0000000
[18,] 9.286366e-09 1.857273e-08 1.0000000
[19,] 2.382375e-09 4.764751e-09 1.0000000
[20,] 5.490997e-10 1.098199e-09 1.0000000
[21,] 1.331092e-10 2.662184e-10 1.0000000
[22,] 1.183425e-10 2.366850e-10 1.0000000
[23,] 2.142421e-10 4.284842e-10 1.0000000
[24,] 8.496132e-11 1.699226e-10 1.0000000
[25,] 2.538465e-11 5.076929e-11 1.0000000
[26,] 1.889492e-11 3.778983e-11 1.0000000
[27,] 9.909244e-12 1.981849e-11 1.0000000
[28,] 4.859553e-12 9.719106e-12 1.0000000
[29,] 1.777231e-12 3.554461e-12 1.0000000
[30,] 3.659502e-12 7.319005e-12 1.0000000
[31,] 4.494615e-12 8.989230e-12 1.0000000
[32,] 7.092156e-12 1.418431e-11 1.0000000
[33,] 8.620384e-11 1.724077e-10 1.0000000
[34,] 5.149808e-10 1.029962e-09 1.0000000
[35,] 9.022131e-09 1.804426e-08 1.0000000
[36,] 9.401217e-07 1.880243e-06 0.9999991
[37,] 2.143307e-05 4.286614e-05 0.9999786
[38,] 5.723303e-04 1.144661e-03 0.9994277
[39,] 3.405736e-03 6.811473e-03 0.9965943
[40,] 1.395721e-02 2.791443e-02 0.9860428
[41,] 4.088631e-02 8.177261e-02 0.9591137
[42,] 1.009893e-01 2.019786e-01 0.8990107
[43,] 1.676137e-01 3.352273e-01 0.8323863
[44,] 2.129940e-01 4.259881e-01 0.7870060
[45,] 2.346637e-01 4.693274e-01 0.7653363
[46,] 2.095707e-01 4.191413e-01 0.7904293
[47,] 1.446089e-01 2.892179e-01 0.8553911
[48,] 1.222898e-01 2.445796e-01 0.8777102
[49,] 1.097411e-01 2.194821e-01 0.8902589
[50,] 1.390661e-01 2.781323e-01 0.8609339
[51,] 4.996712e-01 9.993424e-01 0.5003288
[52,] 6.046652e-01 7.906696e-01 0.3953348
> postscript(file="/var/www/html/rcomp/tmp/12qca1258648001.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/21b2p1258648001.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/3olbr1258648001.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/4v0l21258648001.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/5ucew1258648001.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
-1.619133021 -1.306351982 -0.004198960 0.549726320 1.397918058 -0.037980122
7 8 9 10 11 12
-1.075439835 -0.728038799 -0.002270159 0.650896776 0.308705824 -0.409061619
13 14 15 16 17 18
-0.102599368 -0.143591240 0.723506269 1.299194975 0.222305331 -0.066711537
19 20 21 22 23 24
0.392866640 0.403201825 0.922931064 0.349026117 0.759920756 0.582285548
25 26 27 28 29 30
0.086016757 -0.450441688 -0.405508370 -0.430733448 -0.799579777 -0.868172813
31 32 33 34 35 36
-0.826731258 -1.316015443 -0.844410222 -0.803564793 -0.622985064 -1.333925792
37 38 39 40 41 42
-1.408492020 -1.856001849 -1.842790625 -1.473291841 -0.739719353 -0.457675848
43 44 45 46 47 48
-0.165996258 0.428192212 0.899708318 1.466105009 0.793355906 1.359882491
49 50 51 52 53 54
1.742055592 2.338881674 2.975587080 3.570075899 2.655797330 0.725523962
55 56 57 58 59 60
0.093661734 0.056708858 1.327962910 0.567130969 0.409173173 -0.820664119
61 62 63 64 65
-2.101067283 -2.572697366 -0.996669038 -0.676234461 -0.749560006
> postscript(file="/var/www/html/rcomp/tmp/6ruh01258648001.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 -1.619133021 NA
1 -1.306351982 -1.619133021
2 -0.004198960 -1.306351982
3 0.549726320 -0.004198960
4 1.397918058 0.549726320
5 -0.037980122 1.397918058
6 -1.075439835 -0.037980122
7 -0.728038799 -1.075439835
8 -0.002270159 -0.728038799
9 0.650896776 -0.002270159
10 0.308705824 0.650896776
11 -0.409061619 0.308705824
12 -0.102599368 -0.409061619
13 -0.143591240 -0.102599368
14 0.723506269 -0.143591240
15 1.299194975 0.723506269
16 0.222305331 1.299194975
17 -0.066711537 0.222305331
18 0.392866640 -0.066711537
19 0.403201825 0.392866640
20 0.922931064 0.403201825
21 0.349026117 0.922931064
22 0.759920756 0.349026117
23 0.582285548 0.759920756
24 0.086016757 0.582285548
25 -0.450441688 0.086016757
26 -0.405508370 -0.450441688
27 -0.430733448 -0.405508370
28 -0.799579777 -0.430733448
29 -0.868172813 -0.799579777
30 -0.826731258 -0.868172813
31 -1.316015443 -0.826731258
32 -0.844410222 -1.316015443
33 -0.803564793 -0.844410222
34 -0.622985064 -0.803564793
35 -1.333925792 -0.622985064
36 -1.408492020 -1.333925792
37 -1.856001849 -1.408492020
38 -1.842790625 -1.856001849
39 -1.473291841 -1.842790625
40 -0.739719353 -1.473291841
41 -0.457675848 -0.739719353
42 -0.165996258 -0.457675848
43 0.428192212 -0.165996258
44 0.899708318 0.428192212
45 1.466105009 0.899708318
46 0.793355906 1.466105009
47 1.359882491 0.793355906
48 1.742055592 1.359882491
49 2.338881674 1.742055592
50 2.975587080 2.338881674
51 3.570075899 2.975587080
52 2.655797330 3.570075899
53 0.725523962 2.655797330
54 0.093661734 0.725523962
55 0.056708858 0.093661734
56 1.327962910 0.056708858
57 0.567130969 1.327962910
58 0.409173173 0.567130969
59 -0.820664119 0.409173173
60 -2.101067283 -0.820664119
61 -2.572697366 -2.101067283
62 -0.996669038 -2.572697366
63 -0.676234461 -0.996669038
64 -0.749560006 -0.676234461
65 NA -0.749560006
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.306351982 -1.619133021
[2,] -0.004198960 -1.306351982
[3,] 0.549726320 -0.004198960
[4,] 1.397918058 0.549726320
[5,] -0.037980122 1.397918058
[6,] -1.075439835 -0.037980122
[7,] -0.728038799 -1.075439835
[8,] -0.002270159 -0.728038799
[9,] 0.650896776 -0.002270159
[10,] 0.308705824 0.650896776
[11,] -0.409061619 0.308705824
[12,] -0.102599368 -0.409061619
[13,] -0.143591240 -0.102599368
[14,] 0.723506269 -0.143591240
[15,] 1.299194975 0.723506269
[16,] 0.222305331 1.299194975
[17,] -0.066711537 0.222305331
[18,] 0.392866640 -0.066711537
[19,] 0.403201825 0.392866640
[20,] 0.922931064 0.403201825
[21,] 0.349026117 0.922931064
[22,] 0.759920756 0.349026117
[23,] 0.582285548 0.759920756
[24,] 0.086016757 0.582285548
[25,] -0.450441688 0.086016757
[26,] -0.405508370 -0.450441688
[27,] -0.430733448 -0.405508370
[28,] -0.799579777 -0.430733448
[29,] -0.868172813 -0.799579777
[30,] -0.826731258 -0.868172813
[31,] -1.316015443 -0.826731258
[32,] -0.844410222 -1.316015443
[33,] -0.803564793 -0.844410222
[34,] -0.622985064 -0.803564793
[35,] -1.333925792 -0.622985064
[36,] -1.408492020 -1.333925792
[37,] -1.856001849 -1.408492020
[38,] -1.842790625 -1.856001849
[39,] -1.473291841 -1.842790625
[40,] -0.739719353 -1.473291841
[41,] -0.457675848 -0.739719353
[42,] -0.165996258 -0.457675848
[43,] 0.428192212 -0.165996258
[44,] 0.899708318 0.428192212
[45,] 1.466105009 0.899708318
[46,] 0.793355906 1.466105009
[47,] 1.359882491 0.793355906
[48,] 1.742055592 1.359882491
[49,] 2.338881674 1.742055592
[50,] 2.975587080 2.338881674
[51,] 3.570075899 2.975587080
[52,] 2.655797330 3.570075899
[53,] 0.725523962 2.655797330
[54,] 0.093661734 0.725523962
[55,] 0.056708858 0.093661734
[56,] 1.327962910 0.056708858
[57,] 0.567130969 1.327962910
[58,] 0.409173173 0.567130969
[59,] -0.820664119 0.409173173
[60,] -2.101067283 -0.820664119
[61,] -2.572697366 -2.101067283
[62,] -0.996669038 -2.572697366
[63,] -0.676234461 -0.996669038
[64,] -0.749560006 -0.676234461
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.306351982 -1.619133021
2 -0.004198960 -1.306351982
3 0.549726320 -0.004198960
4 1.397918058 0.549726320
5 -0.037980122 1.397918058
6 -1.075439835 -0.037980122
7 -0.728038799 -1.075439835
8 -0.002270159 -0.728038799
9 0.650896776 -0.002270159
10 0.308705824 0.650896776
11 -0.409061619 0.308705824
12 -0.102599368 -0.409061619
13 -0.143591240 -0.102599368
14 0.723506269 -0.143591240
15 1.299194975 0.723506269
16 0.222305331 1.299194975
17 -0.066711537 0.222305331
18 0.392866640 -0.066711537
19 0.403201825 0.392866640
20 0.922931064 0.403201825
21 0.349026117 0.922931064
22 0.759920756 0.349026117
23 0.582285548 0.759920756
24 0.086016757 0.582285548
25 -0.450441688 0.086016757
26 -0.405508370 -0.450441688
27 -0.430733448 -0.405508370
28 -0.799579777 -0.430733448
29 -0.868172813 -0.799579777
30 -0.826731258 -0.868172813
31 -1.316015443 -0.826731258
32 -0.844410222 -1.316015443
33 -0.803564793 -0.844410222
34 -0.622985064 -0.803564793
35 -1.333925792 -0.622985064
36 -1.408492020 -1.333925792
37 -1.856001849 -1.408492020
38 -1.842790625 -1.856001849
39 -1.473291841 -1.842790625
40 -0.739719353 -1.473291841
41 -0.457675848 -0.739719353
42 -0.165996258 -0.457675848
43 0.428192212 -0.165996258
44 0.899708318 0.428192212
45 1.466105009 0.899708318
46 0.793355906 1.466105009
47 1.359882491 0.793355906
48 1.742055592 1.359882491
49 2.338881674 1.742055592
50 2.975587080 2.338881674
51 3.570075899 2.975587080
52 2.655797330 3.570075899
53 0.725523962 2.655797330
54 0.093661734 0.725523962
55 0.056708858 0.093661734
56 1.327962910 0.056708858
57 0.567130969 1.327962910
58 0.409173173 0.567130969
59 -0.820664119 0.409173173
60 -2.101067283 -0.820664119
61 -2.572697366 -2.101067283
62 -0.996669038 -2.572697366
63 -0.676234461 -0.996669038
64 -0.749560006 -0.676234461
> 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/789qw1258648001.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/8lpee1258648001.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/9upqg1258648001.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/10js881258648001.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/11qtzz1258648001.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/12g6zd1258648001.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/132z0z1258648001.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/1435no1258648001.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/152xi11258648001.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/16zk2y1258648001.tab")
+ }
>
> system("convert tmp/12qca1258648001.ps tmp/12qca1258648001.png")
> system("convert tmp/21b2p1258648001.ps tmp/21b2p1258648001.png")
> system("convert tmp/3olbr1258648001.ps tmp/3olbr1258648001.png")
> system("convert tmp/4v0l21258648001.ps tmp/4v0l21258648001.png")
> system("convert tmp/5ucew1258648001.ps tmp/5ucew1258648001.png")
> system("convert tmp/6ruh01258648001.ps tmp/6ruh01258648001.png")
> system("convert tmp/789qw1258648001.ps tmp/789qw1258648001.png")
> system("convert tmp/8lpee1258648001.ps tmp/8lpee1258648001.png")
> system("convert tmp/9upqg1258648001.ps tmp/9upqg1258648001.png")
> system("convert tmp/10js881258648001.ps tmp/10js881258648001.png")
>
>
> proc.time()
user system elapsed
2.531 1.574 3.025