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(130
+ ,0
+ ,135
+ ,139
+ ,149
+ ,161
+ ,127
+ ,0
+ ,130
+ ,135
+ ,139
+ ,149
+ ,122
+ ,0
+ ,127
+ ,130
+ ,135
+ ,139
+ ,117
+ ,0
+ ,122
+ ,127
+ ,130
+ ,135
+ ,112
+ ,0
+ ,117
+ ,122
+ ,127
+ ,130
+ ,113
+ ,0
+ ,112
+ ,117
+ ,122
+ ,127
+ ,149
+ ,0
+ ,113
+ ,112
+ ,117
+ ,122
+ ,157
+ ,0
+ ,149
+ ,113
+ ,112
+ ,117
+ ,157
+ ,0
+ ,157
+ ,149
+ ,113
+ ,112
+ ,147
+ ,0
+ ,157
+ ,157
+ ,149
+ ,113
+ ,137
+ ,0
+ ,147
+ ,157
+ ,157
+ ,149
+ ,132
+ ,0
+ ,137
+ ,147
+ ,157
+ ,157
+ ,125
+ ,0
+ ,132
+ ,137
+ ,147
+ ,157
+ ,123
+ ,0
+ ,125
+ ,132
+ ,137
+ ,147
+ ,117
+ ,0
+ ,123
+ ,125
+ ,132
+ ,137
+ ,114
+ ,0
+ ,117
+ ,123
+ ,125
+ ,132
+ ,111
+ ,0
+ ,114
+ ,117
+ ,123
+ ,125
+ ,112
+ ,0
+ ,111
+ ,114
+ ,117
+ ,123
+ ,144
+ ,0
+ ,112
+ ,111
+ ,114
+ ,117
+ ,150
+ ,0
+ ,144
+ ,112
+ ,111
+ ,114
+ ,149
+ ,0
+ ,150
+ ,144
+ ,112
+ ,111
+ ,134
+ ,0
+ ,149
+ ,150
+ ,144
+ ,112
+ ,123
+ ,0
+ ,134
+ ,149
+ ,150
+ ,144
+ ,116
+ ,0
+ ,123
+ ,134
+ ,149
+ ,150
+ ,117
+ ,0
+ ,116
+ ,123
+ ,134
+ ,149
+ ,111
+ ,0
+ ,117
+ ,116
+ ,123
+ ,134
+ ,105
+ ,0
+ ,111
+ ,117
+ ,116
+ ,123
+ ,102
+ ,0
+ ,105
+ ,111
+ ,117
+ ,116
+ ,95
+ ,0
+ ,102
+ ,105
+ ,111
+ ,117
+ ,93
+ ,0
+ ,95
+ ,102
+ ,105
+ ,111
+ ,124
+ ,0
+ ,93
+ ,95
+ ,102
+ ,105
+ ,130
+ ,0
+ ,124
+ ,93
+ ,95
+ ,102
+ ,124
+ ,0
+ ,130
+ ,124
+ ,93
+ ,95
+ ,115
+ ,0
+ ,124
+ ,130
+ ,124
+ ,93
+ ,106
+ ,0
+ ,115
+ ,124
+ ,130
+ ,124
+ ,105
+ ,0
+ ,106
+ ,115
+ ,124
+ ,130
+ ,105
+ ,1
+ ,105
+ ,106
+ ,115
+ ,124
+ ,101
+ ,1
+ ,105
+ ,105
+ ,106
+ ,115
+ ,95
+ ,1
+ ,101
+ ,105
+ ,105
+ ,106
+ ,93
+ ,1
+ ,95
+ ,101
+ ,105
+ ,105
+ ,84
+ ,1
+ ,93
+ ,95
+ ,101
+ ,105
+ ,87
+ ,1
+ ,84
+ ,93
+ ,95
+ ,101
+ ,116
+ ,1
+ ,87
+ ,84
+ ,93
+ ,95
+ ,120
+ ,1
+ ,116
+ ,87
+ ,84
+ ,93
+ ,117
+ ,1
+ ,120
+ ,116
+ ,87
+ ,84
+ ,109
+ ,1
+ ,117
+ ,120
+ ,116
+ ,87
+ ,105
+ ,1
+ ,109
+ ,117
+ ,120
+ ,116
+ ,107
+ ,1
+ ,105
+ ,109
+ ,117
+ ,120
+ ,109
+ ,1
+ ,107
+ ,105
+ ,109
+ ,117
+ ,109
+ ,1
+ ,109
+ ,107
+ ,105
+ ,109
+ ,108
+ ,1
+ ,109
+ ,109
+ ,107
+ ,105
+ ,107
+ ,1
+ ,108
+ ,109
+ ,109
+ ,107
+ ,99
+ ,1
+ ,107
+ ,108
+ ,109
+ ,109
+ ,103
+ ,1
+ ,99
+ ,107
+ ,108
+ ,109
+ ,131
+ ,1
+ ,103
+ ,99
+ ,107
+ ,108
+ ,137
+ ,1
+ ,131
+ ,103
+ ,99
+ ,107)
+ ,dim=c(6
+ ,56)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3'
+ ,'Y4')
+ ,1:56))
> y <- array(NA,dim=c(6,56),dimnames=list(c('Y','X','Y1','Y2','Y3','Y4'),1:56))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly 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
Y X Y1 Y2 Y3 Y4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 130 0 135 139 149 161 1 0 0 0 0 0 0 0 0 0 0 1
2 127 0 130 135 139 149 0 1 0 0 0 0 0 0 0 0 0 2
3 122 0 127 130 135 139 0 0 1 0 0 0 0 0 0 0 0 3
4 117 0 122 127 130 135 0 0 0 1 0 0 0 0 0 0 0 4
5 112 0 117 122 127 130 0 0 0 0 1 0 0 0 0 0 0 5
6 113 0 112 117 122 127 0 0 0 0 0 1 0 0 0 0 0 6
7 149 0 113 112 117 122 0 0 0 0 0 0 1 0 0 0 0 7
8 157 0 149 113 112 117 0 0 0 0 0 0 0 1 0 0 0 8
9 157 0 157 149 113 112 0 0 0 0 0 0 0 0 1 0 0 9
10 147 0 157 157 149 113 0 0 0 0 0 0 0 0 0 1 0 10
11 137 0 147 157 157 149 0 0 0 0 0 0 0 0 0 0 1 11
12 132 0 137 147 157 157 0 0 0 0 0 0 0 0 0 0 0 12
13 125 0 132 137 147 157 1 0 0 0 0 0 0 0 0 0 0 13
14 123 0 125 132 137 147 0 1 0 0 0 0 0 0 0 0 0 14
15 117 0 123 125 132 137 0 0 1 0 0 0 0 0 0 0 0 15
16 114 0 117 123 125 132 0 0 0 1 0 0 0 0 0 0 0 16
17 111 0 114 117 123 125 0 0 0 0 1 0 0 0 0 0 0 17
18 112 0 111 114 117 123 0 0 0 0 0 1 0 0 0 0 0 18
19 144 0 112 111 114 117 0 0 0 0 0 0 1 0 0 0 0 19
20 150 0 144 112 111 114 0 0 0 0 0 0 0 1 0 0 0 20
21 149 0 150 144 112 111 0 0 0 0 0 0 0 0 1 0 0 21
22 134 0 149 150 144 112 0 0 0 0 0 0 0 0 0 1 0 22
23 123 0 134 149 150 144 0 0 0 0 0 0 0 0 0 0 1 23
24 116 0 123 134 149 150 0 0 0 0 0 0 0 0 0 0 0 24
25 117 0 116 123 134 149 1 0 0 0 0 0 0 0 0 0 0 25
26 111 0 117 116 123 134 0 1 0 0 0 0 0 0 0 0 0 26
27 105 0 111 117 116 123 0 0 1 0 0 0 0 0 0 0 0 27
28 102 0 105 111 117 116 0 0 0 1 0 0 0 0 0 0 0 28
29 95 0 102 105 111 117 0 0 0 0 1 0 0 0 0 0 0 29
30 93 0 95 102 105 111 0 0 0 0 0 1 0 0 0 0 0 30
31 124 0 93 95 102 105 0 0 0 0 0 0 1 0 0 0 0 31
32 130 0 124 93 95 102 0 0 0 0 0 0 0 1 0 0 0 32
33 124 0 130 124 93 95 0 0 0 0 0 0 0 0 1 0 0 33
34 115 0 124 130 124 93 0 0 0 0 0 0 0 0 0 1 0 34
35 106 0 115 124 130 124 0 0 0 0 0 0 0 0 0 0 1 35
36 105 0 106 115 124 130 0 0 0 0 0 0 0 0 0 0 0 36
37 105 1 105 106 115 124 1 0 0 0 0 0 0 0 0 0 0 37
38 101 1 105 105 106 115 0 1 0 0 0 0 0 0 0 0 0 38
39 95 1 101 105 105 106 0 0 1 0 0 0 0 0 0 0 0 39
40 93 1 95 101 105 105 0 0 0 1 0 0 0 0 0 0 0 40
41 84 1 93 95 101 105 0 0 0 0 1 0 0 0 0 0 0 41
42 87 1 84 93 95 101 0 0 0 0 0 1 0 0 0 0 0 42
43 116 1 87 84 93 95 0 0 0 0 0 0 1 0 0 0 0 43
44 120 1 116 87 84 93 0 0 0 0 0 0 0 1 0 0 0 44
45 117 1 120 116 87 84 0 0 0 0 0 0 0 0 1 0 0 45
46 109 1 117 120 116 87 0 0 0 0 0 0 0 0 0 1 0 46
47 105 1 109 117 120 116 0 0 0 0 0 0 0 0 0 0 1 47
48 107 1 105 109 117 120 0 0 0 0 0 0 0 0 0 0 0 48
49 109 1 107 105 109 117 1 0 0 0 0 0 0 0 0 0 0 49
50 109 1 109 107 105 109 0 1 0 0 0 0 0 0 0 0 0 50
51 108 1 109 109 107 105 0 0 1 0 0 0 0 0 0 0 0 51
52 107 1 108 109 109 107 0 0 0 1 0 0 0 0 0 0 0 52
53 99 1 107 108 109 109 0 0 0 0 1 0 0 0 0 0 0 53
54 103 1 99 107 108 109 0 0 0 0 0 1 0 0 0 0 0 54
55 131 1 103 99 107 108 0 0 0 0 0 0 1 0 0 0 0 55
56 137 1 131 103 99 107 0 0 0 0 0 0 0 1 0 0 0 56
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 Y3 Y4
18.77594 0.43898 1.06333 0.24961 -0.34966 -0.07581
M1 M2 M3 M4 M5 M6
-0.40873 -4.56723 -7.36961 -5.09917 -8.43471 -1.36185
M7 M8 M9 M10 M11 t
28.67850 -1.22253 -18.43288 -16.45172 -8.56438 -0.07445
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3.7222 -1.7476 0.2134 1.2111 4.5959
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.77594 9.86485 1.903 0.06459 .
X 0.43898 1.45633 0.301 0.76473
Y1 1.06333 0.16220 6.556 9.88e-08 ***
Y2 0.24961 0.23670 1.055 0.29829
Y3 -0.34966 0.24421 -1.432 0.16038
Y4 -0.07581 0.18460 -0.411 0.68361
M1 -0.40873 2.59320 -0.158 0.87559
M2 -4.56723 2.70862 -1.686 0.09996 .
M3 -7.36961 2.48626 -2.964 0.00522 **
M4 -5.09917 2.40403 -2.121 0.04050 *
M5 -8.43471 2.38102 -3.542 0.00107 **
M6 -1.36185 2.59008 -0.526 0.60209
M7 28.67850 2.71199 10.575 7.06e-13 ***
M8 -1.22253 6.71886 -0.182 0.85659
M9 -18.43288 7.22139 -2.553 0.01484 *
M10 -16.45172 6.61876 -2.486 0.01745 *
M11 -8.56438 2.68245 -3.193 0.00283 **
t -0.07445 0.05423 -1.373 0.17784
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.569 on 38 degrees of freedom
Multiple R-squared: 0.9849, Adjusted R-squared: 0.9781
F-statistic: 145.6 on 17 and 38 DF, p-value: < 2.2e-16
> 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.2499587 0.4999175 0.75004127
[2,] 0.7745526 0.4508948 0.22544740
[3,] 0.6642796 0.6714407 0.33572036
[4,] 0.5619266 0.8761468 0.43807338
[5,] 0.7661158 0.4677685 0.23388423
[6,] 0.7724170 0.4551660 0.22758300
[7,] 0.7554187 0.4891625 0.24458126
[8,] 0.6699025 0.6601950 0.33009750
[9,] 0.7712082 0.4575835 0.22879177
[10,] 0.7834255 0.4331490 0.21657451
[11,] 0.9094508 0.1810984 0.09054919
[12,] 0.9214875 0.1570249 0.07851247
[13,] 0.8930900 0.2138200 0.10691002
[14,] 0.9096069 0.1807863 0.09039313
[15,] 0.8550021 0.2899959 0.14499793
> postscript(file="/var/www/html/rcomp/tmp/1ot2m1258711373.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/2jz0c1258711373.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/34c5n1258711373.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/4ml1f1258711373.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/55p1n1258711373.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 = 56
Frequency = 1
1 2 3 4 5 6
-2.23319845 0.90850320 1.06662432 -2.11543982 1.43121726 0.02176663
7 8 9 10 11 12
4.11322355 1.43182693 1.19458868 -0.04535724 -1.69848499 -1.45251826
13 14 15 16 17 18
-3.65320737 3.01638667 1.26073708 0.11725708 3.98490669 -0.32427980
19 20 21 22 23 24
-0.10853022 0.31434050 2.35382817 -3.72221508 -1.81164454 -1.75557411
25 26 27 28 29 30
4.59588613 -1.47065457 -1.74500673 0.75563822 -0.16888167 -3.52798645
31 32 33 34 35 36
-0.12380061 0.71258396 -3.35054520 1.31297131 0.01580458 0.69923414
37 38 39 40 41 42
0.45143642 -2.89525522 -2.79706175 0.30955882 -3.05477657 0.61478867
43 44 45 46 47 48
-2.44879314 -3.35731846 -0.19787165 2.45460101 3.49432495 2.50885823
49 50 51 52 53 54
0.83908328 0.44101992 2.21470707 0.93298571 -2.19246571 3.21571094
55 56
-1.43209958 0.89856707
> postscript(file="/var/www/html/rcomp/tmp/6gqko1258711373.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 = 56
Frequency = 1
lag(myerror, k = 1) myerror
0 -2.23319845 NA
1 0.90850320 -2.23319845
2 1.06662432 0.90850320
3 -2.11543982 1.06662432
4 1.43121726 -2.11543982
5 0.02176663 1.43121726
6 4.11322355 0.02176663
7 1.43182693 4.11322355
8 1.19458868 1.43182693
9 -0.04535724 1.19458868
10 -1.69848499 -0.04535724
11 -1.45251826 -1.69848499
12 -3.65320737 -1.45251826
13 3.01638667 -3.65320737
14 1.26073708 3.01638667
15 0.11725708 1.26073708
16 3.98490669 0.11725708
17 -0.32427980 3.98490669
18 -0.10853022 -0.32427980
19 0.31434050 -0.10853022
20 2.35382817 0.31434050
21 -3.72221508 2.35382817
22 -1.81164454 -3.72221508
23 -1.75557411 -1.81164454
24 4.59588613 -1.75557411
25 -1.47065457 4.59588613
26 -1.74500673 -1.47065457
27 0.75563822 -1.74500673
28 -0.16888167 0.75563822
29 -3.52798645 -0.16888167
30 -0.12380061 -3.52798645
31 0.71258396 -0.12380061
32 -3.35054520 0.71258396
33 1.31297131 -3.35054520
34 0.01580458 1.31297131
35 0.69923414 0.01580458
36 0.45143642 0.69923414
37 -2.89525522 0.45143642
38 -2.79706175 -2.89525522
39 0.30955882 -2.79706175
40 -3.05477657 0.30955882
41 0.61478867 -3.05477657
42 -2.44879314 0.61478867
43 -3.35731846 -2.44879314
44 -0.19787165 -3.35731846
45 2.45460101 -0.19787165
46 3.49432495 2.45460101
47 2.50885823 3.49432495
48 0.83908328 2.50885823
49 0.44101992 0.83908328
50 2.21470707 0.44101992
51 0.93298571 2.21470707
52 -2.19246571 0.93298571
53 3.21571094 -2.19246571
54 -1.43209958 3.21571094
55 0.89856707 -1.43209958
56 NA 0.89856707
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.90850320 -2.23319845
[2,] 1.06662432 0.90850320
[3,] -2.11543982 1.06662432
[4,] 1.43121726 -2.11543982
[5,] 0.02176663 1.43121726
[6,] 4.11322355 0.02176663
[7,] 1.43182693 4.11322355
[8,] 1.19458868 1.43182693
[9,] -0.04535724 1.19458868
[10,] -1.69848499 -0.04535724
[11,] -1.45251826 -1.69848499
[12,] -3.65320737 -1.45251826
[13,] 3.01638667 -3.65320737
[14,] 1.26073708 3.01638667
[15,] 0.11725708 1.26073708
[16,] 3.98490669 0.11725708
[17,] -0.32427980 3.98490669
[18,] -0.10853022 -0.32427980
[19,] 0.31434050 -0.10853022
[20,] 2.35382817 0.31434050
[21,] -3.72221508 2.35382817
[22,] -1.81164454 -3.72221508
[23,] -1.75557411 -1.81164454
[24,] 4.59588613 -1.75557411
[25,] -1.47065457 4.59588613
[26,] -1.74500673 -1.47065457
[27,] 0.75563822 -1.74500673
[28,] -0.16888167 0.75563822
[29,] -3.52798645 -0.16888167
[30,] -0.12380061 -3.52798645
[31,] 0.71258396 -0.12380061
[32,] -3.35054520 0.71258396
[33,] 1.31297131 -3.35054520
[34,] 0.01580458 1.31297131
[35,] 0.69923414 0.01580458
[36,] 0.45143642 0.69923414
[37,] -2.89525522 0.45143642
[38,] -2.79706175 -2.89525522
[39,] 0.30955882 -2.79706175
[40,] -3.05477657 0.30955882
[41,] 0.61478867 -3.05477657
[42,] -2.44879314 0.61478867
[43,] -3.35731846 -2.44879314
[44,] -0.19787165 -3.35731846
[45,] 2.45460101 -0.19787165
[46,] 3.49432495 2.45460101
[47,] 2.50885823 3.49432495
[48,] 0.83908328 2.50885823
[49,] 0.44101992 0.83908328
[50,] 2.21470707 0.44101992
[51,] 0.93298571 2.21470707
[52,] -2.19246571 0.93298571
[53,] 3.21571094 -2.19246571
[54,] -1.43209958 3.21571094
[55,] 0.89856707 -1.43209958
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.90850320 -2.23319845
2 1.06662432 0.90850320
3 -2.11543982 1.06662432
4 1.43121726 -2.11543982
5 0.02176663 1.43121726
6 4.11322355 0.02176663
7 1.43182693 4.11322355
8 1.19458868 1.43182693
9 -0.04535724 1.19458868
10 -1.69848499 -0.04535724
11 -1.45251826 -1.69848499
12 -3.65320737 -1.45251826
13 3.01638667 -3.65320737
14 1.26073708 3.01638667
15 0.11725708 1.26073708
16 3.98490669 0.11725708
17 -0.32427980 3.98490669
18 -0.10853022 -0.32427980
19 0.31434050 -0.10853022
20 2.35382817 0.31434050
21 -3.72221508 2.35382817
22 -1.81164454 -3.72221508
23 -1.75557411 -1.81164454
24 4.59588613 -1.75557411
25 -1.47065457 4.59588613
26 -1.74500673 -1.47065457
27 0.75563822 -1.74500673
28 -0.16888167 0.75563822
29 -3.52798645 -0.16888167
30 -0.12380061 -3.52798645
31 0.71258396 -0.12380061
32 -3.35054520 0.71258396
33 1.31297131 -3.35054520
34 0.01580458 1.31297131
35 0.69923414 0.01580458
36 0.45143642 0.69923414
37 -2.89525522 0.45143642
38 -2.79706175 -2.89525522
39 0.30955882 -2.79706175
40 -3.05477657 0.30955882
41 0.61478867 -3.05477657
42 -2.44879314 0.61478867
43 -3.35731846 -2.44879314
44 -0.19787165 -3.35731846
45 2.45460101 -0.19787165
46 3.49432495 2.45460101
47 2.50885823 3.49432495
48 0.83908328 2.50885823
49 0.44101992 0.83908328
50 2.21470707 0.44101992
51 0.93298571 2.21470707
52 -2.19246571 0.93298571
53 3.21571094 -2.19246571
54 -1.43209958 3.21571094
55 0.89856707 -1.43209958
> 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/7ncqr1258711373.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/8ra1e1258711373.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/9cyml1258711373.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/10y2q01258711373.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/11lmar1258711373.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/12h62t1258711373.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/137l111258711373.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/14089z1258711373.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/158f8m1258711373.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/1654ia1258711373.tab")
+ }
>
> system("convert tmp/1ot2m1258711373.ps tmp/1ot2m1258711373.png")
> system("convert tmp/2jz0c1258711373.ps tmp/2jz0c1258711373.png")
> system("convert tmp/34c5n1258711373.ps tmp/34c5n1258711373.png")
> system("convert tmp/4ml1f1258711373.ps tmp/4ml1f1258711373.png")
> system("convert tmp/55p1n1258711373.ps tmp/55p1n1258711373.png")
> system("convert tmp/6gqko1258711373.ps tmp/6gqko1258711373.png")
> system("convert tmp/7ncqr1258711373.ps tmp/7ncqr1258711373.png")
> system("convert tmp/8ra1e1258711373.ps tmp/8ra1e1258711373.png")
> system("convert tmp/9cyml1258711373.ps tmp/9cyml1258711373.png")
> system("convert tmp/10y2q01258711373.ps tmp/10y2q01258711373.png")
>
>
> proc.time()
user system elapsed
2.239 1.550 3.637