R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
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(6,100,6,9,9,99,2,8,7,108,4,3,8,103,0,4,1,99,8,7,9,115,0,7,9,90,8,1,7,95,9,9,2,114,4,4,9,108,2,9,8,112,6,3,3,109,1,3,0,105,0,3,7,105,0,2,5,118,5,8,7,103,7,6,9,112,5,2,6,116,6,6,4,96,6,6,5,101,9,0,8,116,5,4,5,119,3,9,9,115,4,5,0,108,5,2,0,111,5,8,3,108,8,3,8,121,8,9,1,109,6,8,3,112,2,8,2,119,6,8,5,104,1,5,2,105,3,4,5,115,0,4,4,124,1,1,3,116,8,6,0,107,5,2,7,115,6,1,8,116,2,3,8,116,3,8,3,119,0,9,1,111,9,1,9,118,6,7,0,106,9,2,8,103,2,5,8,118,6,0,7,118,7,5,4,102,8,0,3,100,6,1,0,94,9,6,2,94,5,3,1,102,9,9,1,95,3,3,8,92,5,5,7,102,7,8,6,91,5,7,1,89,5,4,5,104,2,8,1,105,2,1,1,99,0,2,7,95,5,0,3,90,5,8,8,96,1,7,5,113,0,5,7,101,9,0,5,101,4,9,7,113,6,8,2,96,6,2,4,97,8,2,0,114,9,9,0,112,5,5,5,108,4,9,3,107,0,0,1,103,5,9,1,107,5,0,3,122,3,9),dim=c(4,75),dimnames=list(c('steenkool','aardolie','uranium','metaal
'),1:75))
> y <- array(NA,dim=c(4,75),dimnames=list(c('steenkool','aardolie','uranium','metaal
'),1:75))
> 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 = '4'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
metaal\r steenkool aardolie uranium
1 9 6 100 6
2 8 9 99 2
3 3 7 108 4
4 4 8 103 0
5 7 1 99 8
6 7 9 115 0
7 1 9 90 8
8 9 7 95 9
9 4 2 114 4
10 9 9 108 2
11 3 8 112 6
12 3 3 109 1
13 3 0 105 0
14 2 7 105 0
15 8 5 118 5
16 6 7 103 7
17 2 9 112 5
18 6 6 116 6
19 6 4 96 6
20 0 5 101 9
21 4 8 116 5
22 9 5 119 3
23 5 9 115 4
24 2 0 108 5
25 8 0 111 5
26 3 3 108 8
27 9 8 121 8
28 8 1 109 6
29 8 3 112 2
30 8 2 119 6
31 5 5 104 1
32 4 2 105 3
33 4 5 115 0
34 1 4 124 1
35 6 3 116 8
36 2 0 107 5
37 1 7 115 6
38 3 8 116 2
39 8 8 116 3
40 9 3 119 0
41 1 1 111 9
42 7 9 118 6
43 2 0 106 9
44 5 8 103 2
45 0 8 118 6
46 5 7 118 7
47 0 4 102 8
48 1 3 100 6
49 6 0 94 9
50 3 2 94 5
51 9 1 102 9
52 3 1 95 3
53 5 8 92 5
54 8 7 102 7
55 7 6 91 5
56 4 1 89 5
57 8 5 104 2
58 1 1 105 2
59 2 1 99 0
60 0 7 95 5
61 8 3 90 5
62 7 8 96 1
63 5 5 113 0
64 0 7 101 9
65 9 5 101 4
66 8 7 113 6
67 2 2 96 6
68 2 4 97 8
69 9 0 114 9
70 5 0 112 5
71 9 5 108 4
72 0 3 107 0
73 9 1 103 5
74 0 1 107 5
75 9 3 122 3
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) steenkool aardolie uranium
-0.57169 0.06576 0.04788 0.01009
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.6644 -2.4722 0.0941 2.8810 4.5317
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.57169 4.55984 -0.125 0.901
steenkool 0.06576 0.12073 0.545 0.588
aardolie 0.04788 0.04146 1.155 0.252
uranium 0.01009 0.13066 0.077 0.939
Residual standard error: 3.092 on 71 degrees of freedom
Multiple R-squared: 0.02485, Adjusted R-squared: -0.01635
F-statistic: 0.6032 on 3 and 71 DF, p-value: 0.6151
> 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.72799538 0.54400923 0.2720046
[2,] 0.77898103 0.44203795 0.2210190
[3,] 0.74139544 0.51720912 0.2586046
[4,] 0.72925487 0.54149026 0.2707451
[5,] 0.68193016 0.63613968 0.3180698
[6,] 0.65276064 0.69447872 0.3472394
[7,] 0.57533975 0.84932050 0.4246603
[8,] 0.55673214 0.88653571 0.4432679
[9,] 0.50471273 0.99057454 0.4952873
[10,] 0.41492378 0.82984755 0.5850762
[11,] 0.46458643 0.92917287 0.5354136
[12,] 0.37889746 0.75779493 0.6211025
[13,] 0.30760640 0.61521280 0.6923936
[14,] 0.45927676 0.91855352 0.5407232
[15,] 0.39204511 0.78409023 0.6079549
[16,] 0.42194271 0.84388542 0.5780573
[17,] 0.34816888 0.69633775 0.6518311
[18,] 0.32772814 0.65545629 0.6722719
[19,] 0.32991913 0.65983825 0.6700809
[20,] 0.28982552 0.57965105 0.7101745
[21,] 0.28576405 0.57152811 0.7142359
[22,] 0.28087995 0.56175990 0.7191200
[23,] 0.26853255 0.53706509 0.7314675
[24,] 0.24186591 0.48373182 0.7581341
[25,] 0.18931589 0.37863179 0.8106841
[26,] 0.14782335 0.29564669 0.8521767
[27,] 0.11728527 0.23457054 0.8827147
[28,] 0.17016589 0.34033178 0.8298341
[29,] 0.13171224 0.26342447 0.8682878
[30,] 0.12406529 0.24813058 0.8759347
[31,] 0.16607211 0.33214423 0.8339279
[32,] 0.15083020 0.30166039 0.8491698
[33,] 0.13555681 0.27111362 0.8644432
[34,] 0.14886145 0.29772290 0.8511385
[35,] 0.17090729 0.34181457 0.8290927
[36,] 0.13643919 0.27287838 0.8635608
[37,] 0.12396809 0.24793618 0.8760319
[38,] 0.09184358 0.18368717 0.9081564
[39,] 0.17547088 0.35094176 0.8245291
[40,] 0.14250394 0.28500789 0.8574961
[41,] 0.21177068 0.42354136 0.7882293
[42,] 0.23077553 0.46155106 0.7692245
[43,] 0.19782540 0.39565079 0.8021746
[44,] 0.15447910 0.30895820 0.8455209
[45,] 0.18723066 0.37446132 0.8127693
[46,] 0.14302644 0.28605288 0.8569736
[47,] 0.10488387 0.20976774 0.8951161
[48,] 0.09177622 0.18355243 0.9082238
[49,] 0.08447581 0.16895162 0.9155242
[50,] 0.06152141 0.12304282 0.9384786
[51,] 0.05847382 0.11694764 0.9415262
[52,] 0.05950302 0.11900605 0.9404970
[53,] 0.04556208 0.09112417 0.9544379
[54,] 0.05905910 0.11811821 0.9409409
[55,] 0.08835565 0.17671131 0.9116443
[56,] 0.09171217 0.18342435 0.9082878
[57,] 0.05933980 0.11867960 0.9406602
[58,] 0.14115636 0.28231272 0.8588436
[59,] 0.22116834 0.44233669 0.7788317
[60,] 0.15713462 0.31426923 0.8428654
[61,] 0.09266161 0.18532322 0.9073384
[62,] 0.13163524 0.26327048 0.8683648
> postscript(file="/var/wessaorg/rcomp/tmp/1pvph1353424948.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/24hb91353424948.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/3004l1353424948.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/4q75f1353424948.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/5ca6a1353424948.ps",horizontal=F,onefile=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 = 75
Frequency = 1
1 2 3 4 5 6
4.32891924 3.21987255 -2.09968276 -0.88570904 2.68540350 1.47400837
7 8 9 10 11 12
-3.40974719 4.47229226 -1.05816755 3.78897741 -2.37711908 -1.85427964
13 14 15 16 17 18
-1.45541768 -2.91570779 2.54297050 1.10944598 -3.43278899 0.56288343
19 20 21 22 23 24
1.65193965 -4.68345971 -1.55854221 3.51526489 -0.56633489 -2.64947847
25 26 27 28 29 30
3.20688981 -1.87700311 3.17181415 3.22680274 2.99200283 2.68227463
31 32 33 34 35 36
0.25359509 -0.61718659 -1.26296871 -4.63819394 0.73997899 -2.60160124
37 38 39 40 41 42
-4.45499506 -2.52828477 2.46162942 3.67703379 -3.89920918 1.26986176
43 44 45 46 47 48
-2.59406726 0.09411933 -5.66438251 -0.60871259 -4.65549541 -3.47381357
49 50 51 52 53 54
1.98045960 -1.11070860 4.53168597 -1.07265848 0.59051150 3.15732322
55 56 57 58 59 60
2.76990020 0.19443332 3.24350928 -3.54134504 -2.23390998 -4.48736448
61 62 63 64 65 66
4.01504463 2.43934581 -0.16721424 -4.81497117 4.36696936 2.64075941
67 68 69 70 71 72
-2.21654889 -2.41610922 4.02291484 0.15901257 4.03182869 -4.74843935
73 74 75
4.52415199 -4.66735697 3.50314463
> postscript(file="/var/wessaorg/rcomp/tmp/6bcw51353424948.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 75
Frequency = 1
lag(myerror, k = 1) myerror
0 4.32891924 NA
1 3.21987255 4.32891924
2 -2.09968276 3.21987255
3 -0.88570904 -2.09968276
4 2.68540350 -0.88570904
5 1.47400837 2.68540350
6 -3.40974719 1.47400837
7 4.47229226 -3.40974719
8 -1.05816755 4.47229226
9 3.78897741 -1.05816755
10 -2.37711908 3.78897741
11 -1.85427964 -2.37711908
12 -1.45541768 -1.85427964
13 -2.91570779 -1.45541768
14 2.54297050 -2.91570779
15 1.10944598 2.54297050
16 -3.43278899 1.10944598
17 0.56288343 -3.43278899
18 1.65193965 0.56288343
19 -4.68345971 1.65193965
20 -1.55854221 -4.68345971
21 3.51526489 -1.55854221
22 -0.56633489 3.51526489
23 -2.64947847 -0.56633489
24 3.20688981 -2.64947847
25 -1.87700311 3.20688981
26 3.17181415 -1.87700311
27 3.22680274 3.17181415
28 2.99200283 3.22680274
29 2.68227463 2.99200283
30 0.25359509 2.68227463
31 -0.61718659 0.25359509
32 -1.26296871 -0.61718659
33 -4.63819394 -1.26296871
34 0.73997899 -4.63819394
35 -2.60160124 0.73997899
36 -4.45499506 -2.60160124
37 -2.52828477 -4.45499506
38 2.46162942 -2.52828477
39 3.67703379 2.46162942
40 -3.89920918 3.67703379
41 1.26986176 -3.89920918
42 -2.59406726 1.26986176
43 0.09411933 -2.59406726
44 -5.66438251 0.09411933
45 -0.60871259 -5.66438251
46 -4.65549541 -0.60871259
47 -3.47381357 -4.65549541
48 1.98045960 -3.47381357
49 -1.11070860 1.98045960
50 4.53168597 -1.11070860
51 -1.07265848 4.53168597
52 0.59051150 -1.07265848
53 3.15732322 0.59051150
54 2.76990020 3.15732322
55 0.19443332 2.76990020
56 3.24350928 0.19443332
57 -3.54134504 3.24350928
58 -2.23390998 -3.54134504
59 -4.48736448 -2.23390998
60 4.01504463 -4.48736448
61 2.43934581 4.01504463
62 -0.16721424 2.43934581
63 -4.81497117 -0.16721424
64 4.36696936 -4.81497117
65 2.64075941 4.36696936
66 -2.21654889 2.64075941
67 -2.41610922 -2.21654889
68 4.02291484 -2.41610922
69 0.15901257 4.02291484
70 4.03182869 0.15901257
71 -4.74843935 4.03182869
72 4.52415199 -4.74843935
73 -4.66735697 4.52415199
74 3.50314463 -4.66735697
75 NA 3.50314463
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.21987255 4.32891924
[2,] -2.09968276 3.21987255
[3,] -0.88570904 -2.09968276
[4,] 2.68540350 -0.88570904
[5,] 1.47400837 2.68540350
[6,] -3.40974719 1.47400837
[7,] 4.47229226 -3.40974719
[8,] -1.05816755 4.47229226
[9,] 3.78897741 -1.05816755
[10,] -2.37711908 3.78897741
[11,] -1.85427964 -2.37711908
[12,] -1.45541768 -1.85427964
[13,] -2.91570779 -1.45541768
[14,] 2.54297050 -2.91570779
[15,] 1.10944598 2.54297050
[16,] -3.43278899 1.10944598
[17,] 0.56288343 -3.43278899
[18,] 1.65193965 0.56288343
[19,] -4.68345971 1.65193965
[20,] -1.55854221 -4.68345971
[21,] 3.51526489 -1.55854221
[22,] -0.56633489 3.51526489
[23,] -2.64947847 -0.56633489
[24,] 3.20688981 -2.64947847
[25,] -1.87700311 3.20688981
[26,] 3.17181415 -1.87700311
[27,] 3.22680274 3.17181415
[28,] 2.99200283 3.22680274
[29,] 2.68227463 2.99200283
[30,] 0.25359509 2.68227463
[31,] -0.61718659 0.25359509
[32,] -1.26296871 -0.61718659
[33,] -4.63819394 -1.26296871
[34,] 0.73997899 -4.63819394
[35,] -2.60160124 0.73997899
[36,] -4.45499506 -2.60160124
[37,] -2.52828477 -4.45499506
[38,] 2.46162942 -2.52828477
[39,] 3.67703379 2.46162942
[40,] -3.89920918 3.67703379
[41,] 1.26986176 -3.89920918
[42,] -2.59406726 1.26986176
[43,] 0.09411933 -2.59406726
[44,] -5.66438251 0.09411933
[45,] -0.60871259 -5.66438251
[46,] -4.65549541 -0.60871259
[47,] -3.47381357 -4.65549541
[48,] 1.98045960 -3.47381357
[49,] -1.11070860 1.98045960
[50,] 4.53168597 -1.11070860
[51,] -1.07265848 4.53168597
[52,] 0.59051150 -1.07265848
[53,] 3.15732322 0.59051150
[54,] 2.76990020 3.15732322
[55,] 0.19443332 2.76990020
[56,] 3.24350928 0.19443332
[57,] -3.54134504 3.24350928
[58,] -2.23390998 -3.54134504
[59,] -4.48736448 -2.23390998
[60,] 4.01504463 -4.48736448
[61,] 2.43934581 4.01504463
[62,] -0.16721424 2.43934581
[63,] -4.81497117 -0.16721424
[64,] 4.36696936 -4.81497117
[65,] 2.64075941 4.36696936
[66,] -2.21654889 2.64075941
[67,] -2.41610922 -2.21654889
[68,] 4.02291484 -2.41610922
[69,] 0.15901257 4.02291484
[70,] 4.03182869 0.15901257
[71,] -4.74843935 4.03182869
[72,] 4.52415199 -4.74843935
[73,] -4.66735697 4.52415199
[74,] 3.50314463 -4.66735697
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.21987255 4.32891924
2 -2.09968276 3.21987255
3 -0.88570904 -2.09968276
4 2.68540350 -0.88570904
5 1.47400837 2.68540350
6 -3.40974719 1.47400837
7 4.47229226 -3.40974719
8 -1.05816755 4.47229226
9 3.78897741 -1.05816755
10 -2.37711908 3.78897741
11 -1.85427964 -2.37711908
12 -1.45541768 -1.85427964
13 -2.91570779 -1.45541768
14 2.54297050 -2.91570779
15 1.10944598 2.54297050
16 -3.43278899 1.10944598
17 0.56288343 -3.43278899
18 1.65193965 0.56288343
19 -4.68345971 1.65193965
20 -1.55854221 -4.68345971
21 3.51526489 -1.55854221
22 -0.56633489 3.51526489
23 -2.64947847 -0.56633489
24 3.20688981 -2.64947847
25 -1.87700311 3.20688981
26 3.17181415 -1.87700311
27 3.22680274 3.17181415
28 2.99200283 3.22680274
29 2.68227463 2.99200283
30 0.25359509 2.68227463
31 -0.61718659 0.25359509
32 -1.26296871 -0.61718659
33 -4.63819394 -1.26296871
34 0.73997899 -4.63819394
35 -2.60160124 0.73997899
36 -4.45499506 -2.60160124
37 -2.52828477 -4.45499506
38 2.46162942 -2.52828477
39 3.67703379 2.46162942
40 -3.89920918 3.67703379
41 1.26986176 -3.89920918
42 -2.59406726 1.26986176
43 0.09411933 -2.59406726
44 -5.66438251 0.09411933
45 -0.60871259 -5.66438251
46 -4.65549541 -0.60871259
47 -3.47381357 -4.65549541
48 1.98045960 -3.47381357
49 -1.11070860 1.98045960
50 4.53168597 -1.11070860
51 -1.07265848 4.53168597
52 0.59051150 -1.07265848
53 3.15732322 0.59051150
54 2.76990020 3.15732322
55 0.19443332 2.76990020
56 3.24350928 0.19443332
57 -3.54134504 3.24350928
58 -2.23390998 -3.54134504
59 -4.48736448 -2.23390998
60 4.01504463 -4.48736448
61 2.43934581 4.01504463
62 -0.16721424 2.43934581
63 -4.81497117 -0.16721424
64 4.36696936 -4.81497117
65 2.64075941 4.36696936
66 -2.21654889 2.64075941
67 -2.41610922 -2.21654889
68 4.02291484 -2.41610922
69 0.15901257 4.02291484
70 4.03182869 0.15901257
71 -4.74843935 4.03182869
72 4.52415199 -4.74843935
73 -4.66735697 4.52415199
74 3.50314463 -4.66735697
> 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/wessaorg/rcomp/tmp/7gqz81353424948.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/85hin1353424948.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/9mfm71353424948.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/105of11353424948.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11o7v21353424948.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/wessaorg/rcomp/tmp/12xnd31353424948.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/wessaorg/rcomp/tmp/132mr21353424948.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/wessaorg/rcomp/tmp/142vsb1353424948.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/wessaorg/rcomp/tmp/15cwtn1353424948.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/wessaorg/rcomp/tmp/166jae1353424948.tab")
+ }
>
> try(system("convert tmp/1pvph1353424948.ps tmp/1pvph1353424948.png",intern=TRUE))
character(0)
> try(system("convert tmp/24hb91353424948.ps tmp/24hb91353424948.png",intern=TRUE))
character(0)
> try(system("convert tmp/3004l1353424948.ps tmp/3004l1353424948.png",intern=TRUE))
character(0)
> try(system("convert tmp/4q75f1353424948.ps tmp/4q75f1353424948.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ca6a1353424948.ps tmp/5ca6a1353424948.png",intern=TRUE))
character(0)
> try(system("convert tmp/6bcw51353424948.ps tmp/6bcw51353424948.png",intern=TRUE))
character(0)
> try(system("convert tmp/7gqz81353424948.ps tmp/7gqz81353424948.png",intern=TRUE))
character(0)
> try(system("convert tmp/85hin1353424948.ps tmp/85hin1353424948.png",intern=TRUE))
character(0)
> try(system("convert tmp/9mfm71353424948.ps tmp/9mfm71353424948.png",intern=TRUE))
character(0)
> try(system("convert tmp/105of11353424948.ps tmp/105of11353424948.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.730 0.989 7.740