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(14.2,-0.8,13.5,-0.2,11.9,0.2,14.6,1,15.6,0,14.1,-0.2,14.9,1,14.2,0.4,14.6,1,17.2,1.7,15.4,3.1,14.3,3.3,17.5,3.1,14.5,3.5,14.4,6,16.6,5.7,16.7,4.7,16.6,4.2,16.9,3.6,15.7,4.4,16.4,2.5,18.4,-0.6,16.9,-1.9,16.5,-1.9,18.3,0.7,15.1,-0.9,15.7,-1.7,18.1,-3.1,16.8,-2.1,18.9,0.2,19,1.2,18.1,3.8,17.8,4,21.5,6.6,17.1,5.3,18.7,7.6,19,4.7,16.4,6.6,16.9,4.4,18.6,4.6,19.3,6,19.4,4.8,17.6,4,18.6,2.7,18.1,3,20.4,4.1,18.1,4,19.6,2.7,19.9,2.6,19.2,3.1,17.8,4.4,19.2,3,22,2,21.1,1.3,19.5,1.5,22.2,1.3,20.9,3.2,22.2,1.8,23.5,3.3,21.5,1,24.3,2.4,22.8,0.4,20.3,-0.1,23.7,1.3,23.3,-1.1,19.6,-4.4,18,-7.5,17.3,-12.2,16.8,-14.5,18.2,-16,16.5,-16.7,16,-16.3,18.4,-16.9),dim=c(2,73),dimnames=list(c('Y','X'),1:73))
> y <- array(NA,dim=c(2,73),dimnames=list(c('Y','X'),1:73))
> 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 14.2 -0.8 1 0 0 0 0 0 0 0 0 0 0 1
2 13.5 -0.2 0 1 0 0 0 0 0 0 0 0 0 2
3 11.9 0.2 0 0 1 0 0 0 0 0 0 0 0 3
4 14.6 1.0 0 0 0 1 0 0 0 0 0 0 0 4
5 15.6 0.0 0 0 0 0 1 0 0 0 0 0 0 5
6 14.1 -0.2 0 0 0 0 0 1 0 0 0 0 0 6
7 14.9 1.0 0 0 0 0 0 0 1 0 0 0 0 7
8 14.2 0.4 0 0 0 0 0 0 0 1 0 0 0 8
9 14.6 1.0 0 0 0 0 0 0 0 0 1 0 0 9
10 17.2 1.7 0 0 0 0 0 0 0 0 0 1 0 10
11 15.4 3.1 0 0 0 0 0 0 0 0 0 0 1 11
12 14.3 3.3 0 0 0 0 0 0 0 0 0 0 0 12
13 17.5 3.1 1 0 0 0 0 0 0 0 0 0 0 13
14 14.5 3.5 0 1 0 0 0 0 0 0 0 0 0 14
15 14.4 6.0 0 0 1 0 0 0 0 0 0 0 0 15
16 16.6 5.7 0 0 0 1 0 0 0 0 0 0 0 16
17 16.7 4.7 0 0 0 0 1 0 0 0 0 0 0 17
18 16.6 4.2 0 0 0 0 0 1 0 0 0 0 0 18
19 16.9 3.6 0 0 0 0 0 0 1 0 0 0 0 19
20 15.7 4.4 0 0 0 0 0 0 0 1 0 0 0 20
21 16.4 2.5 0 0 0 0 0 0 0 0 1 0 0 21
22 18.4 -0.6 0 0 0 0 0 0 0 0 0 1 0 22
23 16.9 -1.9 0 0 0 0 0 0 0 0 0 0 1 23
24 16.5 -1.9 0 0 0 0 0 0 0 0 0 0 0 24
25 18.3 0.7 1 0 0 0 0 0 0 0 0 0 0 25
26 15.1 -0.9 0 1 0 0 0 0 0 0 0 0 0 26
27 15.7 -1.7 0 0 1 0 0 0 0 0 0 0 0 27
28 18.1 -3.1 0 0 0 1 0 0 0 0 0 0 0 28
29 16.8 -2.1 0 0 0 0 1 0 0 0 0 0 0 29
30 18.9 0.2 0 0 0 0 0 1 0 0 0 0 0 30
31 19.0 1.2 0 0 0 0 0 0 1 0 0 0 0 31
32 18.1 3.8 0 0 0 0 0 0 0 1 0 0 0 32
33 17.8 4.0 0 0 0 0 0 0 0 0 1 0 0 33
34 21.5 6.6 0 0 0 0 0 0 0 0 0 1 0 34
35 17.1 5.3 0 0 0 0 0 0 0 0 0 0 1 35
36 18.7 7.6 0 0 0 0 0 0 0 0 0 0 0 36
37 19.0 4.7 1 0 0 0 0 0 0 0 0 0 0 37
38 16.4 6.6 0 1 0 0 0 0 0 0 0 0 0 38
39 16.9 4.4 0 0 1 0 0 0 0 0 0 0 0 39
40 18.6 4.6 0 0 0 1 0 0 0 0 0 0 0 40
41 19.3 6.0 0 0 0 0 1 0 0 0 0 0 0 41
42 19.4 4.8 0 0 0 0 0 1 0 0 0 0 0 42
43 17.6 4.0 0 0 0 0 0 0 1 0 0 0 0 43
44 18.6 2.7 0 0 0 0 0 0 0 1 0 0 0 44
45 18.1 3.0 0 0 0 0 0 0 0 0 1 0 0 45
46 20.4 4.1 0 0 0 0 0 0 0 0 0 1 0 46
47 18.1 4.0 0 0 0 0 0 0 0 0 0 0 1 47
48 19.6 2.7 0 0 0 0 0 0 0 0 0 0 0 48
49 19.9 2.6 1 0 0 0 0 0 0 0 0 0 0 49
50 19.2 3.1 0 1 0 0 0 0 0 0 0 0 0 50
51 17.8 4.4 0 0 1 0 0 0 0 0 0 0 0 51
52 19.2 3.0 0 0 0 1 0 0 0 0 0 0 0 52
53 22.0 2.0 0 0 0 0 1 0 0 0 0 0 0 53
54 21.1 1.3 0 0 0 0 0 1 0 0 0 0 0 54
55 19.5 1.5 0 0 0 0 0 0 1 0 0 0 0 55
56 22.2 1.3 0 0 0 0 0 0 0 1 0 0 0 56
57 20.9 3.2 0 0 0 0 0 0 0 0 1 0 0 57
58 22.2 1.8 0 0 0 0 0 0 0 0 0 1 0 58
59 23.5 3.3 0 0 0 0 0 0 0 0 0 0 1 59
60 21.5 1.0 0 0 0 0 0 0 0 0 0 0 0 60
61 24.3 2.4 1 0 0 0 0 0 0 0 0 0 0 61
62 22.8 0.4 0 1 0 0 0 0 0 0 0 0 0 62
63 20.3 -0.1 0 0 1 0 0 0 0 0 0 0 0 63
64 23.7 1.3 0 0 0 1 0 0 0 0 0 0 0 64
65 23.3 -1.1 0 0 0 0 1 0 0 0 0 0 0 65
66 19.6 -4.4 0 0 0 0 0 1 0 0 0 0 0 66
67 18.0 -7.5 0 0 0 0 0 0 1 0 0 0 0 67
68 17.3 -12.2 0 0 0 0 0 0 0 1 0 0 0 68
69 16.8 -14.5 0 0 0 0 0 0 0 0 1 0 0 69
70 18.2 -16.0 0 0 0 0 0 0 0 0 0 1 0 70
71 16.5 -16.7 0 0 0 0 0 0 0 0 0 0 1 71
72 16.0 -16.3 0 0 0 0 0 0 0 0 0 0 0 72
73 18.4 -16.9 1 0 0 0 0 0 0 0 0 0 0 73
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
13.0121 0.2561 1.6176 -0.3686 -1.2653 0.9477
M5 M6 M7 M8 M9 M10
1.4422 0.8123 0.1518 0.2134 -0.1023 2.0658
M11 t
0.2370 0.1169
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.3746 -0.8105 -0.2486 0.6619 2.8086
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.012145 0.587842 22.135 < 2e-16 ***
X 0.256103 0.029271 8.749 3.00e-12 ***
M1 1.617641 0.681200 2.375 0.02083 *
M2 -0.368594 0.710692 -0.519 0.60595
M3 -1.265335 0.710527 -1.781 0.08009 .
M4 0.947683 0.709946 1.335 0.18705
M5 1.442206 0.708683 2.035 0.04635 *
M6 0.812339 0.707550 1.148 0.25556
M7 0.151781 0.706947 0.215 0.83074
M8 0.213377 0.706391 0.302 0.76366
M9 -0.102263 0.706137 -0.145 0.88535
M10 2.065836 0.705959 2.926 0.00486 **
M11 0.236983 0.705858 0.336 0.73826
t 0.116862 0.007561 15.456 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.223 on 59 degrees of freedom
Multiple R-squared: 0.8282, Adjusted R-squared: 0.7903
F-statistic: 21.88 on 13 and 59 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.1839827165 0.3679654329 0.8160173
[2,] 0.0894082296 0.1788164592 0.9105918
[3,] 0.0369825528 0.0739651056 0.9630174
[4,] 0.0168593194 0.0337186389 0.9831407
[5,] 0.0060627413 0.0121254826 0.9939373
[6,] 0.0021387337 0.0042774674 0.9978613
[7,] 0.0008779529 0.0017559057 0.9991220
[8,] 0.0005001272 0.0010002545 0.9994999
[9,] 0.0001759493 0.0003518987 0.9998241
[10,] 0.0003277664 0.0006555327 0.9996722
[11,] 0.0001761703 0.0003523406 0.9998238
[12,] 0.0001392804 0.0002785608 0.9998607
[13,] 0.0003872968 0.0007745937 0.9996127
[14,] 0.0012170251 0.0024340502 0.9987830
[15,] 0.0096615311 0.0193230621 0.9903385
[16,] 0.0079170508 0.0158341016 0.9920829
[17,] 0.0088127857 0.0176255714 0.9911872
[18,] 0.0249667377 0.0499334754 0.9750333
[19,] 0.0631751008 0.1263502016 0.9368249
[20,] 0.0460755161 0.0921510321 0.9539245
[21,] 0.0506429218 0.1012858437 0.9493571
[22,] 0.1003520524 0.2007041049 0.8996479
[23,] 0.1167604439 0.2335208877 0.8832396
[24,] 0.1313307559 0.2626615117 0.8686692
[25,] 0.1016064398 0.2032128796 0.8983936
[26,] 0.0967916324 0.1935832648 0.9032084
[27,] 0.2080658777 0.4161317555 0.7919341
[28,] 0.1559286214 0.3118572428 0.8440714
[29,] 0.1356197888 0.2712395776 0.8643802
[30,] 0.1267821487 0.2535642974 0.8732179
[31,] 0.1718924182 0.3437848363 0.8281076
[32,] 0.1769007528 0.3538015057 0.8230992
[33,] 0.1416571529 0.2833143057 0.8583428
[34,] 0.1482938964 0.2965877927 0.8517061
[35,] 0.1447843796 0.2895687591 0.8552156
[36,] 0.5181580627 0.9636838747 0.4818419
[37,] 0.5743006133 0.8513987735 0.4256994
[38,] 0.4729092178 0.9458184356 0.5270908
[39,] 0.4052691652 0.8105383304 0.5947308
[40,] 0.6274859819 0.7450280362 0.3725140
> postscript(file="/var/www/html/rcomp/tmp/12raq1258723450.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/296o21258723450.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/3ehyo1258723450.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/4l1o01258723450.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/5s8ky1258723450.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 = 73
Frequency = 1
1 2 3 4 5 6
-0.341764900 0.673947119 -0.248615378 -0.083376398 0.561341734 -0.374432771
7 8 9 10 11 12
0.661941020 -0.062855549 0.382262107 0.518029274 0.071476594 -0.959622667
13 14 15 16 17 18
0.557095500 -0.675971895 -0.636350545 -0.689398341 -0.944680209 -0.403623836
19 20 21 22 23 24
0.593735230 -0.989605441 0.395769540 0.904727843 1.449653075 1.169774400
25 26 27 28 29 30
0.569404362 -0.351457172 1.233303847 1.661969275 -0.505518454 1.518449714
31 32 33 34 35 36
1.906044092 0.161718146 0.009276974 0.758448573 -1.596626195 -0.465541610
37 38 39 40 41 42
-1.157345531 -2.374567321 -0.531262199 -1.212361460 -1.482290361 -0.561961936
43 44 45 46 47 48
-1.613382284 -0.458906802 -0.836958267 -1.103632272 -1.666030557 0.287024578
49 50 51 52 53 54
-1.121867548 -0.080545236 -1.033600370 -1.604934943 0.839783189 0.632060149
55 56 57 58 59 60
-0.475463130 2.097299129 0.509482975 -0.116933703 2.510903323 1.220061388
61 62 63 64 65 66
1.927014867 2.808594505 1.216524645 1.928101868 1.531364102 -0.810491320
67 68 69 70 71 72
-1.072874928 -0.747649482 -0.459833329 -0.960639715 -0.769376241 -1.251696088
73
-0.432536749
> postscript(file="/var/www/html/rcomp/tmp/6kakk1258723450.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.341764900 NA
1 0.673947119 -0.341764900
2 -0.248615378 0.673947119
3 -0.083376398 -0.248615378
4 0.561341734 -0.083376398
5 -0.374432771 0.561341734
6 0.661941020 -0.374432771
7 -0.062855549 0.661941020
8 0.382262107 -0.062855549
9 0.518029274 0.382262107
10 0.071476594 0.518029274
11 -0.959622667 0.071476594
12 0.557095500 -0.959622667
13 -0.675971895 0.557095500
14 -0.636350545 -0.675971895
15 -0.689398341 -0.636350545
16 -0.944680209 -0.689398341
17 -0.403623836 -0.944680209
18 0.593735230 -0.403623836
19 -0.989605441 0.593735230
20 0.395769540 -0.989605441
21 0.904727843 0.395769540
22 1.449653075 0.904727843
23 1.169774400 1.449653075
24 0.569404362 1.169774400
25 -0.351457172 0.569404362
26 1.233303847 -0.351457172
27 1.661969275 1.233303847
28 -0.505518454 1.661969275
29 1.518449714 -0.505518454
30 1.906044092 1.518449714
31 0.161718146 1.906044092
32 0.009276974 0.161718146
33 0.758448573 0.009276974
34 -1.596626195 0.758448573
35 -0.465541610 -1.596626195
36 -1.157345531 -0.465541610
37 -2.374567321 -1.157345531
38 -0.531262199 -2.374567321
39 -1.212361460 -0.531262199
40 -1.482290361 -1.212361460
41 -0.561961936 -1.482290361
42 -1.613382284 -0.561961936
43 -0.458906802 -1.613382284
44 -0.836958267 -0.458906802
45 -1.103632272 -0.836958267
46 -1.666030557 -1.103632272
47 0.287024578 -1.666030557
48 -1.121867548 0.287024578
49 -0.080545236 -1.121867548
50 -1.033600370 -0.080545236
51 -1.604934943 -1.033600370
52 0.839783189 -1.604934943
53 0.632060149 0.839783189
54 -0.475463130 0.632060149
55 2.097299129 -0.475463130
56 0.509482975 2.097299129
57 -0.116933703 0.509482975
58 2.510903323 -0.116933703
59 1.220061388 2.510903323
60 1.927014867 1.220061388
61 2.808594505 1.927014867
62 1.216524645 2.808594505
63 1.928101868 1.216524645
64 1.531364102 1.928101868
65 -0.810491320 1.531364102
66 -1.072874928 -0.810491320
67 -0.747649482 -1.072874928
68 -0.459833329 -0.747649482
69 -0.960639715 -0.459833329
70 -0.769376241 -0.960639715
71 -1.251696088 -0.769376241
72 -0.432536749 -1.251696088
73 NA -0.432536749
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.673947119 -0.341764900
[2,] -0.248615378 0.673947119
[3,] -0.083376398 -0.248615378
[4,] 0.561341734 -0.083376398
[5,] -0.374432771 0.561341734
[6,] 0.661941020 -0.374432771
[7,] -0.062855549 0.661941020
[8,] 0.382262107 -0.062855549
[9,] 0.518029274 0.382262107
[10,] 0.071476594 0.518029274
[11,] -0.959622667 0.071476594
[12,] 0.557095500 -0.959622667
[13,] -0.675971895 0.557095500
[14,] -0.636350545 -0.675971895
[15,] -0.689398341 -0.636350545
[16,] -0.944680209 -0.689398341
[17,] -0.403623836 -0.944680209
[18,] 0.593735230 -0.403623836
[19,] -0.989605441 0.593735230
[20,] 0.395769540 -0.989605441
[21,] 0.904727843 0.395769540
[22,] 1.449653075 0.904727843
[23,] 1.169774400 1.449653075
[24,] 0.569404362 1.169774400
[25,] -0.351457172 0.569404362
[26,] 1.233303847 -0.351457172
[27,] 1.661969275 1.233303847
[28,] -0.505518454 1.661969275
[29,] 1.518449714 -0.505518454
[30,] 1.906044092 1.518449714
[31,] 0.161718146 1.906044092
[32,] 0.009276974 0.161718146
[33,] 0.758448573 0.009276974
[34,] -1.596626195 0.758448573
[35,] -0.465541610 -1.596626195
[36,] -1.157345531 -0.465541610
[37,] -2.374567321 -1.157345531
[38,] -0.531262199 -2.374567321
[39,] -1.212361460 -0.531262199
[40,] -1.482290361 -1.212361460
[41,] -0.561961936 -1.482290361
[42,] -1.613382284 -0.561961936
[43,] -0.458906802 -1.613382284
[44,] -0.836958267 -0.458906802
[45,] -1.103632272 -0.836958267
[46,] -1.666030557 -1.103632272
[47,] 0.287024578 -1.666030557
[48,] -1.121867548 0.287024578
[49,] -0.080545236 -1.121867548
[50,] -1.033600370 -0.080545236
[51,] -1.604934943 -1.033600370
[52,] 0.839783189 -1.604934943
[53,] 0.632060149 0.839783189
[54,] -0.475463130 0.632060149
[55,] 2.097299129 -0.475463130
[56,] 0.509482975 2.097299129
[57,] -0.116933703 0.509482975
[58,] 2.510903323 -0.116933703
[59,] 1.220061388 2.510903323
[60,] 1.927014867 1.220061388
[61,] 2.808594505 1.927014867
[62,] 1.216524645 2.808594505
[63,] 1.928101868 1.216524645
[64,] 1.531364102 1.928101868
[65,] -0.810491320 1.531364102
[66,] -1.072874928 -0.810491320
[67,] -0.747649482 -1.072874928
[68,] -0.459833329 -0.747649482
[69,] -0.960639715 -0.459833329
[70,] -0.769376241 -0.960639715
[71,] -1.251696088 -0.769376241
[72,] -0.432536749 -1.251696088
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.673947119 -0.341764900
2 -0.248615378 0.673947119
3 -0.083376398 -0.248615378
4 0.561341734 -0.083376398
5 -0.374432771 0.561341734
6 0.661941020 -0.374432771
7 -0.062855549 0.661941020
8 0.382262107 -0.062855549
9 0.518029274 0.382262107
10 0.071476594 0.518029274
11 -0.959622667 0.071476594
12 0.557095500 -0.959622667
13 -0.675971895 0.557095500
14 -0.636350545 -0.675971895
15 -0.689398341 -0.636350545
16 -0.944680209 -0.689398341
17 -0.403623836 -0.944680209
18 0.593735230 -0.403623836
19 -0.989605441 0.593735230
20 0.395769540 -0.989605441
21 0.904727843 0.395769540
22 1.449653075 0.904727843
23 1.169774400 1.449653075
24 0.569404362 1.169774400
25 -0.351457172 0.569404362
26 1.233303847 -0.351457172
27 1.661969275 1.233303847
28 -0.505518454 1.661969275
29 1.518449714 -0.505518454
30 1.906044092 1.518449714
31 0.161718146 1.906044092
32 0.009276974 0.161718146
33 0.758448573 0.009276974
34 -1.596626195 0.758448573
35 -0.465541610 -1.596626195
36 -1.157345531 -0.465541610
37 -2.374567321 -1.157345531
38 -0.531262199 -2.374567321
39 -1.212361460 -0.531262199
40 -1.482290361 -1.212361460
41 -0.561961936 -1.482290361
42 -1.613382284 -0.561961936
43 -0.458906802 -1.613382284
44 -0.836958267 -0.458906802
45 -1.103632272 -0.836958267
46 -1.666030557 -1.103632272
47 0.287024578 -1.666030557
48 -1.121867548 0.287024578
49 -0.080545236 -1.121867548
50 -1.033600370 -0.080545236
51 -1.604934943 -1.033600370
52 0.839783189 -1.604934943
53 0.632060149 0.839783189
54 -0.475463130 0.632060149
55 2.097299129 -0.475463130
56 0.509482975 2.097299129
57 -0.116933703 0.509482975
58 2.510903323 -0.116933703
59 1.220061388 2.510903323
60 1.927014867 1.220061388
61 2.808594505 1.927014867
62 1.216524645 2.808594505
63 1.928101868 1.216524645
64 1.531364102 1.928101868
65 -0.810491320 1.531364102
66 -1.072874928 -0.810491320
67 -0.747649482 -1.072874928
68 -0.459833329 -0.747649482
69 -0.960639715 -0.459833329
70 -0.769376241 -0.960639715
71 -1.251696088 -0.769376241
72 -0.432536749 -1.251696088
> 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/70bcd1258723450.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/8grre1258723450.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/9blhv1258723450.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/10a4lz1258723450.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/11sf9e1258723450.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/12lxbh1258723450.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/13rys91258723450.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/14k6rz1258723451.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/15wxcw1258723451.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/16a8kr1258723451.tab")
+ }
>
> system("convert tmp/12raq1258723450.ps tmp/12raq1258723450.png")
> system("convert tmp/296o21258723450.ps tmp/296o21258723450.png")
> system("convert tmp/3ehyo1258723450.ps tmp/3ehyo1258723450.png")
> system("convert tmp/4l1o01258723450.ps tmp/4l1o01258723450.png")
> system("convert tmp/5s8ky1258723450.ps tmp/5s8ky1258723450.png")
> system("convert tmp/6kakk1258723450.ps tmp/6kakk1258723450.png")
> system("convert tmp/70bcd1258723450.ps tmp/70bcd1258723450.png")
> system("convert tmp/8grre1258723450.ps tmp/8grre1258723450.png")
> system("convert tmp/9blhv1258723450.ps tmp/9blhv1258723450.png")
> system("convert tmp/10a4lz1258723450.ps tmp/10a4lz1258723450.png")
>
>
> proc.time()
user system elapsed
2.563 1.578 3.794