R version 2.11.1 (2010-05-31)
Copyright (C) 2010 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(8
+ ,78
+ ,284
+ ,9.100000381
+ ,109
+ ,9.300000191
+ ,68
+ ,433
+ ,8.699999809
+ ,144
+ ,7.5
+ ,70
+ ,739
+ ,7.199999809
+ ,113
+ ,8.899999619
+ ,96
+ ,1792
+ ,8.899999619
+ ,97
+ ,10.19999981
+ ,74
+ ,477
+ ,8.300000191
+ ,206
+ ,8.300000191
+ ,111
+ ,362
+ ,10.89999962
+ ,124
+ ,8.800000191
+ ,77
+ ,671
+ ,10
+ ,152
+ ,8.800000191
+ ,168
+ ,636
+ ,9.100000381
+ ,162
+ ,10.69999981
+ ,82
+ ,329
+ ,8.699999809
+ ,150
+ ,11.69999981
+ ,89
+ ,634
+ ,7.599999905
+ ,134
+ ,8.5
+ ,149
+ ,631
+ ,10.80000019
+ ,292
+ ,8.300000191
+ ,60
+ ,257
+ ,9.5
+ ,108
+ ,8.199999809
+ ,96
+ ,284
+ ,8.800000191
+ ,111
+ ,7.900000095
+ ,83
+ ,603
+ ,9.5
+ ,182
+ ,10.30000019
+ ,130
+ ,686
+ ,8.699999809
+ ,129
+ ,7.400000095
+ ,145
+ ,345
+ ,11.19999981
+ ,158
+ ,9.600000381
+ ,112
+ ,1357
+ ,9.699999809
+ ,186
+ ,9.300000191
+ ,131
+ ,544
+ ,9.600000381
+ ,177
+ ,10.60000038
+ ,80
+ ,205
+ ,9.100000381
+ ,127
+ ,9.699999809
+ ,130
+ ,1264
+ ,9.199999809
+ ,179
+ ,11.60000038
+ ,140
+ ,688
+ ,8.300000191
+ ,80
+ ,8.100000381
+ ,154
+ ,354
+ ,8.399999619
+ ,103
+ ,9.800000191
+ ,118
+ ,1632
+ ,9.399999619
+ ,101
+ ,7.400000095
+ ,94
+ ,348
+ ,9.800000191
+ ,117
+ ,9.399999619
+ ,119
+ ,370
+ ,10.39999962
+ ,88
+ ,11.19999981
+ ,153
+ ,648
+ ,9.899999619
+ ,78
+ ,9.100000381
+ ,116
+ ,366
+ ,9.199999809
+ ,102
+ ,10.5
+ ,97
+ ,540
+ ,10.30000019
+ ,95
+ ,11.89999962
+ ,176
+ ,680
+ ,8.899999619
+ ,80
+ ,8.399999619
+ ,75
+ ,345
+ ,9.600000381
+ ,92
+ ,5
+ ,134
+ ,525
+ ,10.30000019
+ ,126
+ ,9.800000191
+ ,161
+ ,870
+ ,10.39999962
+ ,108
+ ,9.800000191
+ ,111
+ ,669
+ ,9.699999809
+ ,77
+ ,10.80000019
+ ,114
+ ,452
+ ,9.600000381
+ ,60
+ ,10.10000038
+ ,142
+ ,430
+ ,10.69999981
+ ,71
+ ,10.89999962
+ ,238
+ ,822
+ ,10.30000019
+ ,86
+ ,9.199999809
+ ,78
+ ,190
+ ,10.69999981
+ ,93
+ ,8.300000191
+ ,196
+ ,867
+ ,9.600000381
+ ,106
+ ,7.300000191
+ ,125
+ ,969
+ ,10.5
+ ,162
+ ,9.399999619
+ ,82
+ ,499
+ ,7.699999809
+ ,95
+ ,9.399999619
+ ,125
+ ,925
+ ,10.19999981
+ ,91
+ ,9.800000191
+ ,129
+ ,353
+ ,9.899999619
+ ,52
+ ,3.599999905
+ ,84
+ ,288
+ ,8.399999619
+ ,110
+ ,8.399999619
+ ,183
+ ,718
+ ,10.39999962
+ ,69
+ ,10.80000019
+ ,119
+ ,540
+ ,9.199999809
+ ,57
+ ,10.10000038
+ ,180
+ ,668
+ ,13
+ ,106
+ ,9
+ ,82
+ ,347
+ ,8.800000191
+ ,40
+ ,10
+ ,71
+ ,345
+ ,9.199999809
+ ,50
+ ,11.30000019
+ ,118
+ ,463
+ ,7.800000191
+ ,35
+ ,11.30000019
+ ,121
+ ,728
+ ,8.199999809
+ ,86
+ ,12.80000019
+ ,68
+ ,383
+ ,7.400000095
+ ,57
+ ,10
+ ,112
+ ,316
+ ,10.39999962
+ ,57
+ ,6.699999809
+ ,109
+ ,388
+ ,8.899999619
+ ,94)
+ ,dim=c(5
+ ,53)
+ ,dimnames=list(c('X1'
+ ,'X2'
+ ,'X3'
+ ,'X4'
+ ,'X5')
+ ,1:53))
> y <- array(NA,dim=c(5,53),dimnames=list(c('X1','X2','X3','X4','X5'),1:53))
> 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
> 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
X1 X2 X3 X4 X5
1 8.0 78 284 9.1 109
2 9.3 68 433 8.7 144
3 7.5 70 739 7.2 113
4 8.9 96 1792 8.9 97
5 10.2 74 477 8.3 206
6 8.3 111 362 10.9 124
7 8.8 77 671 10.0 152
8 8.8 168 636 9.1 162
9 10.7 82 329 8.7 150
10 11.7 89 634 7.6 134
11 8.5 149 631 10.8 292
12 8.3 60 257 9.5 108
13 8.2 96 284 8.8 111
14 7.9 83 603 9.5 182
15 10.3 130 686 8.7 129
16 7.4 145 345 11.2 158
17 9.6 112 1357 9.7 186
18 9.3 131 544 9.6 177
19 10.6 80 205 9.1 127
20 9.7 130 1264 9.2 179
21 11.6 140 688 8.3 80
22 8.1 154 354 8.4 103
23 9.8 118 1632 9.4 101
24 7.4 94 348 9.8 117
25 9.4 119 370 10.4 88
26 11.2 153 648 9.9 78
27 9.1 116 366 9.2 102
28 10.5 97 540 10.3 95
29 11.9 176 680 8.9 80
30 8.4 75 345 9.6 92
31 5.0 134 525 10.3 126
32 9.8 161 870 10.4 108
33 9.8 111 669 9.7 77
34 10.8 114 452 9.6 60
35 10.1 142 430 10.7 71
36 10.9 238 822 10.3 86
37 9.2 78 190 10.7 93
38 8.3 196 867 9.6 106
39 7.3 125 969 10.5 162
40 9.4 82 499 7.7 95
41 9.4 125 925 10.2 91
42 9.8 129 353 9.9 52
43 3.6 84 288 8.4 110
44 8.4 183 718 10.4 69
45 10.8 119 540 9.2 57
46 10.1 180 668 13.0 106
47 9.0 82 347 8.8 40
48 10.0 71 345 9.2 50
49 11.3 118 463 7.8 35
50 11.3 121 728 8.2 86
51 12.8 68 383 7.4 57
52 10.0 112 316 10.4 57
53 6.7 109 388 8.9 94
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X2 X3 X4 X5
12.2662552 0.0073916 0.0005837 -0.3302302 -0.0094629
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.6404 -0.7904 0.3053 0.9164 2.7906
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.2662552 2.0201467 6.072 1.95e-07 ***
X2 0.0073916 0.0069336 1.066 0.2917
X3 0.0005837 0.0007219 0.809 0.4228
X4 -0.3302302 0.2345518 -1.408 0.1656
X5 -0.0094629 0.0048868 -1.936 0.0587 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.601 on 48 degrees of freedom
Multiple R-squared: 0.1437, Adjusted R-squared: 0.07235
F-statistic: 2.014 on 4 and 48 DF, p-value: 0.1075
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.04735538 0.09471077 0.9526446
[2,] 0.18431886 0.36863772 0.8156811
[3,] 0.45654232 0.91308464 0.5434577
[4,] 0.44680065 0.89360131 0.5531993
[5,] 0.33208288 0.66416577 0.6679171
[6,] 0.24800068 0.49600136 0.7519993
[7,] 0.20672543 0.41345085 0.7932746
[8,] 0.17879923 0.35759846 0.8212008
[9,] 0.12315396 0.24630792 0.8768460
[10,] 0.08844399 0.17688797 0.9115560
[11,] 0.07078761 0.14157521 0.9292124
[12,] 0.13920769 0.27841538 0.8607923
[13,] 0.13791671 0.27583342 0.8620833
[14,] 0.18165772 0.36331544 0.8183423
[15,] 0.19592614 0.39185227 0.8040739
[16,] 0.18234504 0.36469009 0.8176550
[17,] 0.15320961 0.30641921 0.8467904
[18,] 0.13469996 0.26939992 0.8653000
[19,] 0.16569258 0.33138517 0.8343074
[20,] 0.14393095 0.28786191 0.8560690
[21,] 0.16867180 0.33734359 0.8313282
[22,] 0.22046957 0.44093914 0.7795304
[23,] 0.16535669 0.33071337 0.8346433
[24,] 0.39131485 0.78262971 0.6086851
[25,] 0.31706157 0.63412315 0.6829384
[26,] 0.25087979 0.50175957 0.7491202
[27,] 0.21140037 0.42280074 0.7885996
[28,] 0.16518466 0.33036932 0.8348153
[29,] 0.15909670 0.31819340 0.8409033
[30,] 0.15886440 0.31772879 0.8411356
[31,] 0.13129618 0.26259236 0.8687038
[32,] 0.09541635 0.19083270 0.9045837
[33,] 0.07491765 0.14983529 0.9250824
[34,] 0.12493212 0.24986424 0.8750679
[35,] 0.09577038 0.19154075 0.9042296
[36,] 0.31268277 0.62536555 0.6873172
[37,] 0.33328742 0.66657484 0.6667126
[38,] 0.21176195 0.42352390 0.7882381
> postscript(file="/var/www/rcomp/tmp/1v4ym1290178260.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/rcomp/tmp/2v4ym1290178260.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/rcomp/tmp/36vfp1290178260.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/rcomp/tmp/46vfp1290178260.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/rcomp/tmp/56vfp1290178260.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 = 53
Frequency = 1
1 2 3 4 5 6 7
-0.9720267 0.5140247 -2.2680705 -1.2649202 1.7985981 -0.2251221 0.3135783
8 9 10 11 12 13 14
-0.5412068 1.9280255 2.1835915 1.0937186 -0.4005880 -0.9852193 -0.4723074
15 16 17 18 19 20 21
0.7661213 -0.9457068 0.6771120 0.5930432 1.8296360 0.4669926 1.3952646
22 23 24 25 26 27 28
-1.7625878 -0.2311740 -1.4207860 0.3052957 1.4319640 -0.2339897 1.5018973
29 30 31 32 33 34 35
1.6319733 -0.5812128 -3.9694873 0.2922483 0.2546453 1.1652449 0.7384664
36 37 38 39 40 41 42
0.6099054 0.6598045 -1.7478168 -1.4554225 -0.3218953 -0.1006735 0.1355243
43 44 45 46 47 48 49
-5.6404098 -1.5506955 0.9164389 1.4093913 -0.7903754 0.5188206 0.7982710
50 51 52 53
1.2361106 2.7906403 0.6952086 -2.7698629
> postscript(file="/var/www/rcomp/tmp/6g4fa1290178260.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 = 53
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.9720267 NA
1 0.5140247 -0.9720267
2 -2.2680705 0.5140247
3 -1.2649202 -2.2680705
4 1.7985981 -1.2649202
5 -0.2251221 1.7985981
6 0.3135783 -0.2251221
7 -0.5412068 0.3135783
8 1.9280255 -0.5412068
9 2.1835915 1.9280255
10 1.0937186 2.1835915
11 -0.4005880 1.0937186
12 -0.9852193 -0.4005880
13 -0.4723074 -0.9852193
14 0.7661213 -0.4723074
15 -0.9457068 0.7661213
16 0.6771120 -0.9457068
17 0.5930432 0.6771120
18 1.8296360 0.5930432
19 0.4669926 1.8296360
20 1.3952646 0.4669926
21 -1.7625878 1.3952646
22 -0.2311740 -1.7625878
23 -1.4207860 -0.2311740
24 0.3052957 -1.4207860
25 1.4319640 0.3052957
26 -0.2339897 1.4319640
27 1.5018973 -0.2339897
28 1.6319733 1.5018973
29 -0.5812128 1.6319733
30 -3.9694873 -0.5812128
31 0.2922483 -3.9694873
32 0.2546453 0.2922483
33 1.1652449 0.2546453
34 0.7384664 1.1652449
35 0.6099054 0.7384664
36 0.6598045 0.6099054
37 -1.7478168 0.6598045
38 -1.4554225 -1.7478168
39 -0.3218953 -1.4554225
40 -0.1006735 -0.3218953
41 0.1355243 -0.1006735
42 -5.6404098 0.1355243
43 -1.5506955 -5.6404098
44 0.9164389 -1.5506955
45 1.4093913 0.9164389
46 -0.7903754 1.4093913
47 0.5188206 -0.7903754
48 0.7982710 0.5188206
49 1.2361106 0.7982710
50 2.7906403 1.2361106
51 0.6952086 2.7906403
52 -2.7698629 0.6952086
53 NA -2.7698629
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.5140247 -0.9720267
[2,] -2.2680705 0.5140247
[3,] -1.2649202 -2.2680705
[4,] 1.7985981 -1.2649202
[5,] -0.2251221 1.7985981
[6,] 0.3135783 -0.2251221
[7,] -0.5412068 0.3135783
[8,] 1.9280255 -0.5412068
[9,] 2.1835915 1.9280255
[10,] 1.0937186 2.1835915
[11,] -0.4005880 1.0937186
[12,] -0.9852193 -0.4005880
[13,] -0.4723074 -0.9852193
[14,] 0.7661213 -0.4723074
[15,] -0.9457068 0.7661213
[16,] 0.6771120 -0.9457068
[17,] 0.5930432 0.6771120
[18,] 1.8296360 0.5930432
[19,] 0.4669926 1.8296360
[20,] 1.3952646 0.4669926
[21,] -1.7625878 1.3952646
[22,] -0.2311740 -1.7625878
[23,] -1.4207860 -0.2311740
[24,] 0.3052957 -1.4207860
[25,] 1.4319640 0.3052957
[26,] -0.2339897 1.4319640
[27,] 1.5018973 -0.2339897
[28,] 1.6319733 1.5018973
[29,] -0.5812128 1.6319733
[30,] -3.9694873 -0.5812128
[31,] 0.2922483 -3.9694873
[32,] 0.2546453 0.2922483
[33,] 1.1652449 0.2546453
[34,] 0.7384664 1.1652449
[35,] 0.6099054 0.7384664
[36,] 0.6598045 0.6099054
[37,] -1.7478168 0.6598045
[38,] -1.4554225 -1.7478168
[39,] -0.3218953 -1.4554225
[40,] -0.1006735 -0.3218953
[41,] 0.1355243 -0.1006735
[42,] -5.6404098 0.1355243
[43,] -1.5506955 -5.6404098
[44,] 0.9164389 -1.5506955
[45,] 1.4093913 0.9164389
[46,] -0.7903754 1.4093913
[47,] 0.5188206 -0.7903754
[48,] 0.7982710 0.5188206
[49,] 1.2361106 0.7982710
[50,] 2.7906403 1.2361106
[51,] 0.6952086 2.7906403
[52,] -2.7698629 0.6952086
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.5140247 -0.9720267
2 -2.2680705 0.5140247
3 -1.2649202 -2.2680705
4 1.7985981 -1.2649202
5 -0.2251221 1.7985981
6 0.3135783 -0.2251221
7 -0.5412068 0.3135783
8 1.9280255 -0.5412068
9 2.1835915 1.9280255
10 1.0937186 2.1835915
11 -0.4005880 1.0937186
12 -0.9852193 -0.4005880
13 -0.4723074 -0.9852193
14 0.7661213 -0.4723074
15 -0.9457068 0.7661213
16 0.6771120 -0.9457068
17 0.5930432 0.6771120
18 1.8296360 0.5930432
19 0.4669926 1.8296360
20 1.3952646 0.4669926
21 -1.7625878 1.3952646
22 -0.2311740 -1.7625878
23 -1.4207860 -0.2311740
24 0.3052957 -1.4207860
25 1.4319640 0.3052957
26 -0.2339897 1.4319640
27 1.5018973 -0.2339897
28 1.6319733 1.5018973
29 -0.5812128 1.6319733
30 -3.9694873 -0.5812128
31 0.2922483 -3.9694873
32 0.2546453 0.2922483
33 1.1652449 0.2546453
34 0.7384664 1.1652449
35 0.6099054 0.7384664
36 0.6598045 0.6099054
37 -1.7478168 0.6598045
38 -1.4554225 -1.7478168
39 -0.3218953 -1.4554225
40 -0.1006735 -0.3218953
41 0.1355243 -0.1006735
42 -5.6404098 0.1355243
43 -1.5506955 -5.6404098
44 0.9164389 -1.5506955
45 1.4093913 0.9164389
46 -0.7903754 1.4093913
47 0.5188206 -0.7903754
48 0.7982710 0.5188206
49 1.2361106 0.7982710
50 2.7906403 1.2361106
51 0.6952086 2.7906403
52 -2.7698629 0.6952086
> 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/rcomp/tmp/7rwwv1290178260.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/rcomp/tmp/8rwwv1290178260.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/rcomp/tmp/9rwwv1290178260.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/rcomp/tmp/10k5vg1290178260.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/1155c41290178260.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/rcomp/tmp/1286ss1290178260.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/rcomp/tmp/134gqi1290178260.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/rcomp/tmp/148yoo1290178260.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/rcomp/tmp/15tznu1290178260.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/rcomp/tmp/16fzli1290178260.tab")
+ }
>
> try(system("convert tmp/1v4ym1290178260.ps tmp/1v4ym1290178260.png",intern=TRUE))
character(0)
> try(system("convert tmp/2v4ym1290178260.ps tmp/2v4ym1290178260.png",intern=TRUE))
character(0)
> try(system("convert tmp/36vfp1290178260.ps tmp/36vfp1290178260.png",intern=TRUE))
character(0)
> try(system("convert tmp/46vfp1290178260.ps tmp/46vfp1290178260.png",intern=TRUE))
character(0)
> try(system("convert tmp/56vfp1290178260.ps tmp/56vfp1290178260.png",intern=TRUE))
character(0)
> try(system("convert tmp/6g4fa1290178260.ps tmp/6g4fa1290178260.png",intern=TRUE))
character(0)
> try(system("convert tmp/7rwwv1290178260.ps tmp/7rwwv1290178260.png",intern=TRUE))
character(0)
> try(system("convert tmp/8rwwv1290178260.ps tmp/8rwwv1290178260.png",intern=TRUE))
character(0)
> try(system("convert tmp/9rwwv1290178260.ps tmp/9rwwv1290178260.png",intern=TRUE))
character(0)
> try(system("convert tmp/10k5vg1290178260.ps tmp/10k5vg1290178260.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.58 2.04 5.58