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(1
+ ,6
+ ,100
+ ,6
+ ,9
+ ,2
+ ,9
+ ,99
+ ,2
+ ,8
+ ,3
+ ,7
+ ,108
+ ,4
+ ,3
+ ,4
+ ,8
+ ,103
+ ,0
+ ,4
+ ,5
+ ,1
+ ,99
+ ,8
+ ,7
+ ,6
+ ,9
+ ,115
+ ,0
+ ,7
+ ,7
+ ,9
+ ,90
+ ,8
+ ,1
+ ,8
+ ,7
+ ,95
+ ,9
+ ,9
+ ,9
+ ,2
+ ,114
+ ,4
+ ,4
+ ,10
+ ,9
+ ,108
+ ,2
+ ,9
+ ,11
+ ,8
+ ,112
+ ,6
+ ,3
+ ,12
+ ,3
+ ,109
+ ,1
+ ,3
+ ,13
+ ,0
+ ,105
+ ,0
+ ,3
+ ,14
+ ,7
+ ,105
+ ,0
+ ,2
+ ,15
+ ,5
+ ,118
+ ,5
+ ,8
+ ,16
+ ,7
+ ,103
+ ,7
+ ,6
+ ,17
+ ,9
+ ,112
+ ,5
+ ,2
+ ,18
+ ,6
+ ,116
+ ,6
+ ,6
+ ,19
+ ,4
+ ,96
+ ,6
+ ,6
+ ,20
+ ,5
+ ,101
+ ,9
+ ,0
+ ,21
+ ,8
+ ,116
+ ,5
+ ,4
+ ,22
+ ,5
+ ,119
+ ,3
+ ,9
+ ,23
+ ,9
+ ,115
+ ,4
+ ,5
+ ,24
+ ,0
+ ,108
+ ,5
+ ,2
+ ,25
+ ,0
+ ,111
+ ,5
+ ,8
+ ,26
+ ,3
+ ,108
+ ,8
+ ,3
+ ,27
+ ,8
+ ,121
+ ,8
+ ,9
+ ,28
+ ,1
+ ,109
+ ,6
+ ,8
+ ,29
+ ,3
+ ,112
+ ,2
+ ,8
+ ,30
+ ,2
+ ,119
+ ,6
+ ,8
+ ,31
+ ,5
+ ,104
+ ,1
+ ,5
+ ,32
+ ,2
+ ,105
+ ,3
+ ,4
+ ,33
+ ,5
+ ,115
+ ,0
+ ,4
+ ,34
+ ,4
+ ,124
+ ,1
+ ,1
+ ,35
+ ,3
+ ,116
+ ,8
+ ,6
+ ,36
+ ,0
+ ,107
+ ,5
+ ,2
+ ,37
+ ,7
+ ,115
+ ,6
+ ,1
+ ,38
+ ,8
+ ,116
+ ,2
+ ,3
+ ,39
+ ,8
+ ,116
+ ,3
+ ,8
+ ,40
+ ,3
+ ,119
+ ,0
+ ,9
+ ,41
+ ,1
+ ,111
+ ,9
+ ,1
+ ,42
+ ,9
+ ,118
+ ,6
+ ,7
+ ,43
+ ,0
+ ,106
+ ,9
+ ,2
+ ,44
+ ,8
+ ,103
+ ,2
+ ,5
+ ,45
+ ,8
+ ,118
+ ,6
+ ,0
+ ,46
+ ,7
+ ,118
+ ,7
+ ,5
+ ,47
+ ,4
+ ,102
+ ,8
+ ,0
+ ,48
+ ,3
+ ,100
+ ,6
+ ,1
+ ,49
+ ,0
+ ,94
+ ,9
+ ,6
+ ,50
+ ,2
+ ,94
+ ,5
+ ,3
+ ,51
+ ,1
+ ,102
+ ,9
+ ,9
+ ,52
+ ,1
+ ,95
+ ,3
+ ,3
+ ,53
+ ,8
+ ,92
+ ,5
+ ,5
+ ,54
+ ,7
+ ,102
+ ,7
+ ,8
+ ,55
+ ,6
+ ,91
+ ,5
+ ,7
+ ,56
+ ,1
+ ,89
+ ,5
+ ,4
+ ,57
+ ,5
+ ,104
+ ,2
+ ,8
+ ,58
+ ,1
+ ,105
+ ,2
+ ,1
+ ,59
+ ,1
+ ,99
+ ,0
+ ,2
+ ,60
+ ,7
+ ,95
+ ,5
+ ,0
+ ,61
+ ,3
+ ,90
+ ,5
+ ,8
+ ,62
+ ,8
+ ,96
+ ,1
+ ,7
+ ,63
+ ,5
+ ,113
+ ,0
+ ,5
+ ,64
+ ,7
+ ,101
+ ,9
+ ,0
+ ,65
+ ,5
+ ,101
+ ,4
+ ,9
+ ,66
+ ,7
+ ,113
+ ,6
+ ,8
+ ,67
+ ,2
+ ,96
+ ,6
+ ,2
+ ,68
+ ,4
+ ,97
+ ,8
+ ,2
+ ,69
+ ,0
+ ,114
+ ,9
+ ,9
+ ,70
+ ,0
+ ,112
+ ,5
+ ,5
+ ,71
+ ,5
+ ,108
+ ,4
+ ,9
+ ,72
+ ,3
+ ,107
+ ,0
+ ,0
+ ,73
+ ,1
+ ,103
+ ,5
+ ,9
+ ,74
+ ,1
+ ,107
+ ,5
+ ,0
+ ,75
+ ,3
+ ,122
+ ,3
+ ,9)
+ ,dim=c(5
+ ,75)
+ ,dimnames=list(c('maand_i...n'
+ ,'steenkool'
+ ,'aardolie'
+ ,'uranium'
+ ,'metaal')
+ ,1:75))
> y <- array(NA,dim=c(5,75),dimnames=list(c('maand_i...n','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 = '3'
> 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
aardolie maand_i...n steenkool uranium metaal
1 100 1 6 6 9
2 99 2 9 2 8
3 108 3 7 4 3
4 103 4 8 0 4
5 99 5 1 8 7
6 115 6 9 0 7
7 90 7 9 8 1
8 95 8 7 9 9
9 114 9 2 4 4
10 108 10 9 2 9
11 112 11 8 6 3
12 109 12 3 1 3
13 105 13 0 0 3
14 105 14 7 0 2
15 118 15 5 5 8
16 103 16 7 7 6
17 112 17 9 5 2
18 116 18 6 6 6
19 96 19 4 6 6
20 101 20 5 9 0
21 116 21 8 5 4
22 119 22 5 3 9
23 115 23 9 4 5
24 108 24 0 5 2
25 111 25 0 5 8
26 108 26 3 8 3
27 121 27 8 8 9
28 109 28 1 6 8
29 112 29 3 2 8
30 119 30 2 6 8
31 104 31 5 1 5
32 105 32 2 3 4
33 115 33 5 0 4
34 124 34 4 1 1
35 116 35 3 8 6
36 107 36 0 5 2
37 115 37 7 6 1
38 116 38 8 2 3
39 116 39 8 3 8
40 119 40 3 0 9
41 111 41 1 9 1
42 118 42 9 6 7
43 106 43 0 9 2
44 103 44 8 2 5
45 118 45 8 6 0
46 118 46 7 7 5
47 102 47 4 8 0
48 100 48 3 6 1
49 94 49 0 9 6
50 94 50 2 5 3
51 102 51 1 9 9
52 95 52 1 3 3
53 92 53 8 5 5
54 102 54 7 7 8
55 91 55 6 5 7
56 89 56 1 5 4
57 104 57 5 2 8
58 105 58 1 2 1
59 99 59 1 0 2
60 95 60 7 5 0
61 90 61 3 5 8
62 96 62 8 1 7
63 113 63 5 0 5
64 101 64 7 9 0
65 101 65 5 4 9
66 113 66 7 6 8
67 96 67 2 6 2
68 97 68 4 8 2
69 114 69 0 9 9
70 112 70 0 5 5
71 108 71 5 4 9
72 107 72 3 0 0
73 103 73 1 5 9
74 107 74 1 5 0
75 122 75 3 3 9
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) maand_i...n steenkool uranium metaal
108.27802 -0.04899 0.15156 -0.48329 0.37445
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-16.3237 -5.8622 0.2151 6.8363 16.8901
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 108.27802 3.87182 27.966 <2e-16 ***
maand_i...n -0.04899 0.04928 -0.994 0.324
steenkool 0.15156 0.35894 0.422 0.674
uranium -0.48329 0.36604 -1.320 0.191
metaal 0.37445 0.33363 1.122 0.266
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.768 on 70 degrees of freedom
Multiple R-squared: 0.06935, Adjusted R-squared: 0.01617
F-statistic: 1.304 on 4 and 70 DF, p-value: 0.277
> 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.45515995 0.91031990 0.5448401
[2,] 0.30509541 0.61019082 0.6949046
[3,] 0.18796299 0.37592598 0.8120370
[4,] 0.22530912 0.45061824 0.7746909
[5,] 0.27059433 0.54118867 0.7294057
[6,] 0.33433865 0.66867730 0.6656613
[7,] 0.29338537 0.58677073 0.7066146
[8,] 0.31251593 0.62503186 0.6874841
[9,] 0.26328054 0.52656109 0.7367195
[10,] 0.21122522 0.42245044 0.7887748
[11,] 0.17150121 0.34300241 0.8284988
[12,] 0.35280165 0.70560331 0.6471983
[13,] 0.29650728 0.59301457 0.7034927
[14,] 0.25225063 0.50450125 0.7477494
[15,] 0.19592554 0.39185109 0.8040745
[16,] 0.14354142 0.28708284 0.8564586
[17,] 0.10620354 0.21240709 0.8937965
[18,] 0.07703663 0.15407326 0.9229634
[19,] 0.05247422 0.10494844 0.9475258
[20,] 0.04596787 0.09193574 0.9540321
[21,] 0.03671908 0.07343816 0.9632809
[22,] 0.03241846 0.06483693 0.9675815
[23,] 0.02595350 0.05190700 0.9740465
[24,] 0.05728014 0.11456027 0.9427199
[25,] 0.05728202 0.11456404 0.9427180
[26,] 0.03958255 0.07916510 0.9604174
[27,] 0.05537566 0.11075132 0.9446243
[28,] 0.04430000 0.08860001 0.9557000
[29,] 0.03426624 0.06853247 0.9657338
[30,] 0.02679420 0.05358841 0.9732058
[31,] 0.02173910 0.04347819 0.9782609
[32,] 0.01829544 0.03659089 0.9817046
[33,] 0.02094074 0.04188148 0.9790593
[34,] 0.01901422 0.03802843 0.9809858
[35,] 0.02396292 0.04792585 0.9760371
[36,] 0.02296810 0.04593619 0.9770319
[37,] 0.04503104 0.09006209 0.9549690
[38,] 0.11125205 0.22250410 0.8887479
[39,] 0.36440284 0.72880569 0.6355972
[40,] 0.47893318 0.95786636 0.5210668
[41,] 0.59200117 0.81599767 0.4079988
[42,] 0.65677157 0.68645685 0.3432284
[43,] 0.70648340 0.58703320 0.2935166
[44,] 0.72023331 0.55953339 0.2797667
[45,] 0.73173973 0.53652053 0.2682603
[46,] 0.77019954 0.45960093 0.2298005
[47,] 0.78189166 0.43621668 0.2181083
[48,] 0.78982222 0.42035555 0.2101778
[49,] 0.80391583 0.39216834 0.1960842
[50,] 0.75510262 0.48979477 0.2448974
[51,] 0.75982209 0.48035581 0.2401779
[52,] 0.70893991 0.58212018 0.2910601
[53,] 0.63577205 0.72845591 0.3642280
[54,] 0.68641643 0.62716714 0.3135836
[55,] 0.72299749 0.55400502 0.2770025
[56,] 0.73156888 0.53686225 0.2684311
[57,] 0.64131302 0.71737397 0.3586870
[58,] 0.56823513 0.86352974 0.4317649
[59,] 0.58296663 0.83406675 0.4170334
[60,] 0.45456445 0.90912891 0.5454355
> postscript(file="/var/wessaorg/rcomp/tmp/18eju1353426955.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/28gwk1353426955.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/3kus91353426955.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/4rwbl1353426955.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/5df4o1353426955.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
-9.608691552 -12.573069598 -0.382152647 -7.792318842 -7.939498575
6 7 8 9 10
3.030750171 -15.807272426 -12.967483399 6.295090706 -3.555630241
11 12 13 14 15
4.824754654 0.215082225 -3.764552930 -4.402004166 9.119831214
16 17 18 19 20
-4.418820684 4.858278292 8.347419779 -11.300483674 -2.706496177
21 22 23 24 25
8.456879242 9.121710618 6.545559934 2.565177963 3.367468205
26 27 28 29 30
3.283895930 13.328410010 1.846158156 2.658886035 11.792575090
31 32 33 34 35
-4.906191277 -2.061516321 6.082943261 16.890119431 10.601422756
36 37 38 39 40
2.153010959 9.998846657 8.214231155 6.874257638 8.856709771
41 42 43 44 45
8.553983187 10.693970766 3.429061278 -5.240750961 13.613629396
46 47 48 49 50
12.425211111 -0.715603640 -3.856085501 -9.774819451 -10.838743652
51 52 53 54 55
-2.951750439 -10.555790123 -14.350015410 -4.306248146 -15.697831393
56 57 58 59 60
-15.767721228 -4.272614106 0.003738054 -7.288299040 -8.983311063
61 62 63 64 65
-16.323698506 -12.591187016 5.178076442 -0.854218992 -6.288600880
66 67 68 69 70
6.798297915 -7.148244000 -5.435794512 10.081554287 7.695189858
71 72 73 74 75
1.005315618 1.794308188 -2.807204353 4.611825492 15.021083479
> postscript(file="/var/wessaorg/rcomp/tmp/6gxj51353426955.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 -9.608691552 NA
1 -12.573069598 -9.608691552
2 -0.382152647 -12.573069598
3 -7.792318842 -0.382152647
4 -7.939498575 -7.792318842
5 3.030750171 -7.939498575
6 -15.807272426 3.030750171
7 -12.967483399 -15.807272426
8 6.295090706 -12.967483399
9 -3.555630241 6.295090706
10 4.824754654 -3.555630241
11 0.215082225 4.824754654
12 -3.764552930 0.215082225
13 -4.402004166 -3.764552930
14 9.119831214 -4.402004166
15 -4.418820684 9.119831214
16 4.858278292 -4.418820684
17 8.347419779 4.858278292
18 -11.300483674 8.347419779
19 -2.706496177 -11.300483674
20 8.456879242 -2.706496177
21 9.121710618 8.456879242
22 6.545559934 9.121710618
23 2.565177963 6.545559934
24 3.367468205 2.565177963
25 3.283895930 3.367468205
26 13.328410010 3.283895930
27 1.846158156 13.328410010
28 2.658886035 1.846158156
29 11.792575090 2.658886035
30 -4.906191277 11.792575090
31 -2.061516321 -4.906191277
32 6.082943261 -2.061516321
33 16.890119431 6.082943261
34 10.601422756 16.890119431
35 2.153010959 10.601422756
36 9.998846657 2.153010959
37 8.214231155 9.998846657
38 6.874257638 8.214231155
39 8.856709771 6.874257638
40 8.553983187 8.856709771
41 10.693970766 8.553983187
42 3.429061278 10.693970766
43 -5.240750961 3.429061278
44 13.613629396 -5.240750961
45 12.425211111 13.613629396
46 -0.715603640 12.425211111
47 -3.856085501 -0.715603640
48 -9.774819451 -3.856085501
49 -10.838743652 -9.774819451
50 -2.951750439 -10.838743652
51 -10.555790123 -2.951750439
52 -14.350015410 -10.555790123
53 -4.306248146 -14.350015410
54 -15.697831393 -4.306248146
55 -15.767721228 -15.697831393
56 -4.272614106 -15.767721228
57 0.003738054 -4.272614106
58 -7.288299040 0.003738054
59 -8.983311063 -7.288299040
60 -16.323698506 -8.983311063
61 -12.591187016 -16.323698506
62 5.178076442 -12.591187016
63 -0.854218992 5.178076442
64 -6.288600880 -0.854218992
65 6.798297915 -6.288600880
66 -7.148244000 6.798297915
67 -5.435794512 -7.148244000
68 10.081554287 -5.435794512
69 7.695189858 10.081554287
70 1.005315618 7.695189858
71 1.794308188 1.005315618
72 -2.807204353 1.794308188
73 4.611825492 -2.807204353
74 15.021083479 4.611825492
75 NA 15.021083479
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -12.573069598 -9.608691552
[2,] -0.382152647 -12.573069598
[3,] -7.792318842 -0.382152647
[4,] -7.939498575 -7.792318842
[5,] 3.030750171 -7.939498575
[6,] -15.807272426 3.030750171
[7,] -12.967483399 -15.807272426
[8,] 6.295090706 -12.967483399
[9,] -3.555630241 6.295090706
[10,] 4.824754654 -3.555630241
[11,] 0.215082225 4.824754654
[12,] -3.764552930 0.215082225
[13,] -4.402004166 -3.764552930
[14,] 9.119831214 -4.402004166
[15,] -4.418820684 9.119831214
[16,] 4.858278292 -4.418820684
[17,] 8.347419779 4.858278292
[18,] -11.300483674 8.347419779
[19,] -2.706496177 -11.300483674
[20,] 8.456879242 -2.706496177
[21,] 9.121710618 8.456879242
[22,] 6.545559934 9.121710618
[23,] 2.565177963 6.545559934
[24,] 3.367468205 2.565177963
[25,] 3.283895930 3.367468205
[26,] 13.328410010 3.283895930
[27,] 1.846158156 13.328410010
[28,] 2.658886035 1.846158156
[29,] 11.792575090 2.658886035
[30,] -4.906191277 11.792575090
[31,] -2.061516321 -4.906191277
[32,] 6.082943261 -2.061516321
[33,] 16.890119431 6.082943261
[34,] 10.601422756 16.890119431
[35,] 2.153010959 10.601422756
[36,] 9.998846657 2.153010959
[37,] 8.214231155 9.998846657
[38,] 6.874257638 8.214231155
[39,] 8.856709771 6.874257638
[40,] 8.553983187 8.856709771
[41,] 10.693970766 8.553983187
[42,] 3.429061278 10.693970766
[43,] -5.240750961 3.429061278
[44,] 13.613629396 -5.240750961
[45,] 12.425211111 13.613629396
[46,] -0.715603640 12.425211111
[47,] -3.856085501 -0.715603640
[48,] -9.774819451 -3.856085501
[49,] -10.838743652 -9.774819451
[50,] -2.951750439 -10.838743652
[51,] -10.555790123 -2.951750439
[52,] -14.350015410 -10.555790123
[53,] -4.306248146 -14.350015410
[54,] -15.697831393 -4.306248146
[55,] -15.767721228 -15.697831393
[56,] -4.272614106 -15.767721228
[57,] 0.003738054 -4.272614106
[58,] -7.288299040 0.003738054
[59,] -8.983311063 -7.288299040
[60,] -16.323698506 -8.983311063
[61,] -12.591187016 -16.323698506
[62,] 5.178076442 -12.591187016
[63,] -0.854218992 5.178076442
[64,] -6.288600880 -0.854218992
[65,] 6.798297915 -6.288600880
[66,] -7.148244000 6.798297915
[67,] -5.435794512 -7.148244000
[68,] 10.081554287 -5.435794512
[69,] 7.695189858 10.081554287
[70,] 1.005315618 7.695189858
[71,] 1.794308188 1.005315618
[72,] -2.807204353 1.794308188
[73,] 4.611825492 -2.807204353
[74,] 15.021083479 4.611825492
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -12.573069598 -9.608691552
2 -0.382152647 -12.573069598
3 -7.792318842 -0.382152647
4 -7.939498575 -7.792318842
5 3.030750171 -7.939498575
6 -15.807272426 3.030750171
7 -12.967483399 -15.807272426
8 6.295090706 -12.967483399
9 -3.555630241 6.295090706
10 4.824754654 -3.555630241
11 0.215082225 4.824754654
12 -3.764552930 0.215082225
13 -4.402004166 -3.764552930
14 9.119831214 -4.402004166
15 -4.418820684 9.119831214
16 4.858278292 -4.418820684
17 8.347419779 4.858278292
18 -11.300483674 8.347419779
19 -2.706496177 -11.300483674
20 8.456879242 -2.706496177
21 9.121710618 8.456879242
22 6.545559934 9.121710618
23 2.565177963 6.545559934
24 3.367468205 2.565177963
25 3.283895930 3.367468205
26 13.328410010 3.283895930
27 1.846158156 13.328410010
28 2.658886035 1.846158156
29 11.792575090 2.658886035
30 -4.906191277 11.792575090
31 -2.061516321 -4.906191277
32 6.082943261 -2.061516321
33 16.890119431 6.082943261
34 10.601422756 16.890119431
35 2.153010959 10.601422756
36 9.998846657 2.153010959
37 8.214231155 9.998846657
38 6.874257638 8.214231155
39 8.856709771 6.874257638
40 8.553983187 8.856709771
41 10.693970766 8.553983187
42 3.429061278 10.693970766
43 -5.240750961 3.429061278
44 13.613629396 -5.240750961
45 12.425211111 13.613629396
46 -0.715603640 12.425211111
47 -3.856085501 -0.715603640
48 -9.774819451 -3.856085501
49 -10.838743652 -9.774819451
50 -2.951750439 -10.838743652
51 -10.555790123 -2.951750439
52 -14.350015410 -10.555790123
53 -4.306248146 -14.350015410
54 -15.697831393 -4.306248146
55 -15.767721228 -15.697831393
56 -4.272614106 -15.767721228
57 0.003738054 -4.272614106
58 -7.288299040 0.003738054
59 -8.983311063 -7.288299040
60 -16.323698506 -8.983311063
61 -12.591187016 -16.323698506
62 5.178076442 -12.591187016
63 -0.854218992 5.178076442
64 -6.288600880 -0.854218992
65 6.798297915 -6.288600880
66 -7.148244000 6.798297915
67 -5.435794512 -7.148244000
68 10.081554287 -5.435794512
69 7.695189858 10.081554287
70 1.005315618 7.695189858
71 1.794308188 1.005315618
72 -2.807204353 1.794308188
73 4.611825492 -2.807204353
74 15.021083479 4.611825492
> 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/70nhl1353426955.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/8kjw11353426955.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/9ix4u1353426955.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/10v0qb1353426955.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/1153fb1353426955.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/12aj7k1353426955.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/13hg0o1353426955.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/14kqnq1353426955.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/151myw1353426955.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/16seig1353426955.tab")
+ }
>
> try(system("convert tmp/18eju1353426955.ps tmp/18eju1353426955.png",intern=TRUE))
character(0)
> try(system("convert tmp/28gwk1353426955.ps tmp/28gwk1353426955.png",intern=TRUE))
character(0)
> try(system("convert tmp/3kus91353426955.ps tmp/3kus91353426955.png",intern=TRUE))
character(0)
> try(system("convert tmp/4rwbl1353426955.ps tmp/4rwbl1353426955.png",intern=TRUE))
character(0)
> try(system("convert tmp/5df4o1353426955.ps tmp/5df4o1353426955.png",intern=TRUE))
character(0)
> try(system("convert tmp/6gxj51353426955.ps tmp/6gxj51353426955.png",intern=TRUE))
character(0)
> try(system("convert tmp/70nhl1353426955.ps tmp/70nhl1353426955.png",intern=TRUE))
character(0)
> try(system("convert tmp/8kjw11353426955.ps tmp/8kjw11353426955.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ix4u1353426955.ps tmp/9ix4u1353426955.png",intern=TRUE))
character(0)
> try(system("convert tmp/10v0qb1353426955.ps tmp/10v0qb1353426955.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
7.385 1.079 8.473