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(2.0 + ,0.0 + ,2.0 + ,1.7 + ,1.6 + ,1.4 + ,2.1 + ,0.0 + ,2.0 + ,2.0 + ,1.7 + ,1.6 + ,2.5 + ,0.0 + ,2.1 + ,2.0 + ,2.0 + ,1.7 + ,2.5 + ,0.0 + ,2.5 + ,2.1 + ,2.0 + ,2.0 + ,2.6 + ,0.0 + ,2.5 + ,2.5 + ,2.1 + ,2.0 + ,2.7 + ,0.0 + ,2.6 + ,2.5 + ,2.5 + ,2.1 + ,3.7 + ,0.0 + ,2.7 + ,2.6 + ,2.5 + ,2.5 + ,4.0 + ,0.0 + ,3.7 + ,2.7 + ,2.6 + ,2.5 + ,5.0 + ,0.0 + ,4.0 + ,3.7 + ,2.7 + ,2.6 + ,5.1 + ,0.0 + ,5.0 + ,4.0 + ,3.7 + ,2.7 + ,5.1 + ,0.0 + ,5.1 + ,5.0 + ,4.0 + ,3.7 + ,5.0 + ,0.0 + ,5.1 + ,5.1 + ,5.0 + ,4.0 + ,5.1 + ,0.0 + ,5.0 + ,5.1 + ,5.1 + ,5.0 + ,4.7 + ,0.0 + ,5.1 + ,5.0 + ,5.1 + ,5.1 + ,4.5 + ,0.0 + ,4.7 + ,5.1 + ,5.0 + ,5.1 + ,4.5 + ,0.0 + ,4.5 + ,4.7 + ,5.1 + ,5.0 + ,4.6 + ,0.0 + ,4.5 + ,4.5 + ,4.7 + ,5.1 + ,4.6 + ,0.0 + ,4.6 + ,4.5 + ,4.5 + ,4.7 + ,4.6 + ,0.0 + ,4.6 + ,4.6 + ,4.5 + ,4.5 + ,4.6 + ,0.0 + ,4.6 + ,4.6 + ,4.6 + ,4.5 + ,5.3 + ,0.0 + ,4.6 + ,4.6 + ,4.6 + ,4.6 + ,5.4 + ,0.0 + ,5.3 + ,4.6 + ,4.6 + ,4.6 + ,5.3 + ,0.0 + ,5.4 + ,5.3 + ,4.6 + ,4.6 + ,5.2 + ,0.0 + ,5.3 + ,5.4 + ,5.3 + ,4.6 + ,5.0 + ,0.0 + ,5.2 + ,5.3 + ,5.4 + ,5.3 + ,4.2 + ,0.0 + ,5.0 + ,5.2 + ,5.3 + ,5.4 + ,4.3 + ,0.0 + ,4.2 + ,5.0 + ,5.2 + ,5.3 + ,4.3 + ,0.0 + ,4.3 + ,4.2 + ,5.0 + ,5.2 + ,4.3 + ,0.0 + ,4.3 + ,4.3 + ,4.2 + ,5.0 + ,4.0 + ,0.0 + ,4.3 + ,4.3 + ,4.3 + ,4.2 + ,4.0 + ,0.0 + ,4.0 + ,4.3 + ,4.3 + ,4.3 + ,4.1 + ,0.0 + ,4.0 + ,4.0 + ,4.3 + ,4.3 + ,4.4 + ,0.0 + ,4.1 + ,4.0 + ,4.0 + ,4.3 + ,3.6 + ,0.0 + ,4.4 + ,4.1 + ,4.0 + ,4.0 + ,3.7 + ,0.0 + ,3.6 + ,4.4 + ,4.1 + ,4.0 + ,3.8 + ,0.0 + ,3.7 + ,3.6 + ,4.4 + ,4.1 + ,3.3 + ,0.0 + ,3.8 + ,3.7 + ,3.6 + ,4.4 + ,3.3 + ,0.0 + ,3.3 + ,3.8 + ,3.7 + ,3.6 + ,3.3 + ,0.0 + ,3.3 + ,3.3 + ,3.8 + ,3.7 + ,3.5 + ,0.0 + ,3.3 + ,3.3 + ,3.3 + ,3.8 + ,3.3 + ,0.0 + ,3.5 + ,3.3 + ,3.3 + ,3.3 + ,3.3 + ,0.0 + ,3.3 + ,3.5 + ,3.3 + ,3.3 + ,3.4 + ,0.0 + ,3.3 + ,3.3 + ,3.5 + ,3.3 + ,3.4 + ,0.0 + ,3.4 + ,3.3 + ,3.3 + ,3.5 + ,5.2 + ,0.0 + ,3.4 + ,3.4 + ,3.3 + ,3.3 + ,5.3 + ,0.0 + ,5.2 + ,3.4 + ,3.4 + ,3.3 + ,4.8 + ,1.0 + ,5.3 + ,5.2 + ,3.4 + ,3.4 + ,5.0 + ,1.0 + ,4.8 + ,5.3 + ,5.2 + ,3.4 + ,4.6 + ,1.0 + ,5.0 + ,4.8 + ,5.3 + ,5.2 + ,4.6 + ,1.0 + ,4.6 + ,5.0 + ,4.8 + ,5.3 + ,3.5 + ,1.0 + ,4.6 + ,4.6 + ,5.0 + ,4.8 + ,3.5 + ,1.0 + ,3.5 + ,4.6 + ,4.6 + ,5.0) + ,dim=c(6 + ,52) + ,dimnames=list(c('IndGez' + ,'InvlMex' + ,'Yt-1' + ,'Yt-2' + ,'Yt-3' + ,'Yt-4') + ,1:52)) > y <- array(NA,dim=c(6,52),dimnames=list(c('IndGez','InvlMex','Yt-1','Yt-2','Yt-3','Yt-4'),1:52)) > 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 IndGez InvlMex Yt-1 Yt-2 Yt-3 Yt-4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t 1 2.0 0 2.0 1.7 1.6 1.4 1 0 0 0 0 0 0 0 0 0 0 1 2 2.1 0 2.0 2.0 1.7 1.6 0 1 0 0 0 0 0 0 0 0 0 2 3 2.5 0 2.1 2.0 2.0 1.7 0 0 1 0 0 0 0 0 0 0 0 3 4 2.5 0 2.5 2.1 2.0 2.0 0 0 0 1 0 0 0 0 0 0 0 4 5 2.6 0 2.5 2.5 2.1 2.0 0 0 0 0 1 0 0 0 0 0 0 5 6 2.7 0 2.6 2.5 2.5 2.1 0 0 0 0 0 1 0 0 0 0 0 6 7 3.7 0 2.7 2.6 2.5 2.5 0 0 0 0 0 0 1 0 0 0 0 7 8 4.0 0 3.7 2.7 2.6 2.5 0 0 0 0 0 0 0 1 0 0 0 8 9 5.0 0 4.0 3.7 2.7 2.6 0 0 0 0 0 0 0 0 1 0 0 9 10 5.1 0 5.0 4.0 3.7 2.7 0 0 0 0 0 0 0 0 0 1 0 10 11 5.1 0 5.1 5.0 4.0 3.7 0 0 0 0 0 0 0 0 0 0 1 11 12 5.0 0 5.1 5.1 5.0 4.0 0 0 0 0 0 0 0 0 0 0 0 12 13 5.1 0 5.0 5.1 5.1 5.0 1 0 0 0 0 0 0 0 0 0 0 13 14 4.7 0 5.1 5.0 5.1 5.1 0 1 0 0 0 0 0 0 0 0 0 14 15 4.5 0 4.7 5.1 5.0 5.1 0 0 1 0 0 0 0 0 0 0 0 15 16 4.5 0 4.5 4.7 5.1 5.0 0 0 0 1 0 0 0 0 0 0 0 16 17 4.6 0 4.5 4.5 4.7 5.1 0 0 0 0 1 0 0 0 0 0 0 17 18 4.6 0 4.6 4.5 4.5 4.7 0 0 0 0 0 1 0 0 0 0 0 18 19 4.6 0 4.6 4.6 4.5 4.5 0 0 0 0 0 0 1 0 0 0 0 19 20 4.6 0 4.6 4.6 4.6 4.5 0 0 0 0 0 0 0 1 0 0 0 20 21 5.3 0 4.6 4.6 4.6 4.6 0 0 0 0 0 0 0 0 1 0 0 21 22 5.4 0 5.3 4.6 4.6 4.6 0 0 0 0 0 0 0 0 0 1 0 22 23 5.3 0 5.4 5.3 4.6 4.6 0 0 0 0 0 0 0 0 0 0 1 23 24 5.2 0 5.3 5.4 5.3 4.6 0 0 0 0 0 0 0 0 0 0 0 24 25 5.0 0 5.2 5.3 5.4 5.3 1 0 0 0 0 0 0 0 0 0 0 25 26 4.2 0 5.0 5.2 5.3 5.4 0 1 0 0 0 0 0 0 0 0 0 26 27 4.3 0 4.2 5.0 5.2 5.3 0 0 1 0 0 0 0 0 0 0 0 27 28 4.3 0 4.3 4.2 5.0 5.2 0 0 0 1 0 0 0 0 0 0 0 28 29 4.3 0 4.3 4.3 4.2 5.0 0 0 0 0 1 0 0 0 0 0 0 29 30 4.0 0 4.3 4.3 4.3 4.2 0 0 0 0 0 1 0 0 0 0 0 30 31 4.0 0 4.0 4.3 4.3 4.3 0 0 0 0 0 0 1 0 0 0 0 31 32 4.1 0 4.0 4.0 4.3 4.3 0 0 0 0 0 0 0 1 0 0 0 32 33 4.4 0 4.1 4.0 4.0 4.3 0 0 0 0 0 0 0 0 1 0 0 33 34 3.6 0 4.4 4.1 4.0 4.0 0 0 0 0 0 0 0 0 0 1 0 34 35 3.7 0 3.6 4.4 4.1 4.0 0 0 0 0 0 0 0 0 0 0 1 35 36 3.8 0 3.7 3.6 4.4 4.1 0 0 0 0 0 0 0 0 0 0 0 36 37 3.3 0 3.8 3.7 3.6 4.4 1 0 0 0 0 0 0 0 0 0 0 37 38 3.3 0 3.3 3.8 3.7 3.6 0 1 0 0 0 0 0 0 0 0 0 38 39 3.3 0 3.3 3.3 3.8 3.7 0 0 1 0 0 0 0 0 0 0 0 39 40 3.5 0 3.3 3.3 3.3 3.8 0 0 0 1 0 0 0 0 0 0 0 40 41 3.3 0 3.5 3.3 3.3 3.3 0 0 0 0 1 0 0 0 0 0 0 41 42 3.3 0 3.3 3.5 3.3 3.3 0 0 0 0 0 1 0 0 0 0 0 42 43 3.4 0 3.3 3.3 3.5 3.3 0 0 0 0 0 0 1 0 0 0 0 43 44 3.4 0 3.4 3.3 3.3 3.5 0 0 0 0 0 0 0 1 0 0 0 44 45 5.2 0 3.4 3.4 3.3 3.3 0 0 0 0 0 0 0 0 1 0 0 45 46 5.3 0 5.2 3.4 3.4 3.3 0 0 0 0 0 0 0 0 0 1 0 46 47 4.8 1 5.3 5.2 3.4 3.4 0 0 0 0 0 0 0 0 0 0 1 47 48 5.0 1 4.8 5.3 5.2 3.4 0 0 0 0 0 0 0 0 0 0 0 48 49 4.6 1 5.0 4.8 5.3 5.2 1 0 0 0 0 0 0 0 0 0 0 49 50 4.6 1 4.6 5.0 4.8 5.3 0 1 0 0 0 0 0 0 0 0 0 50 51 3.5 1 4.6 4.6 5.0 4.8 0 0 1 0 0 0 0 0 0 0 0 51 52 3.5 1 3.5 4.6 4.6 5.0 0 0 0 1 0 0 0 0 0 0 0 52 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) InvlMex `Yt-1` `Yt-2` `Yt-3` `Yt-4` 0.680001 -0.109917 0.829084 0.040628 0.102037 -0.109824 M1 M2 M3 M4 M5 M6 -0.199445 -0.252552 -0.236224 -0.022686 -0.074551 -0.161684 M7 M8 M9 M10 M11 t 0.160658 0.042935 0.906708 -0.040783 -0.050514 -0.002753 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -0.72902 -0.15564 0.02206 0.13708 0.80584 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.680001 0.419206 1.622 0.11402 InvlMex -0.109917 0.225445 -0.488 0.62899 `Yt-1` 0.829084 0.169489 4.892 2.37e-05 *** `Yt-2` 0.040628 0.247835 0.164 0.87076 `Yt-3` 0.102037 0.258943 0.394 0.69600 `Yt-4` -0.109824 0.189244 -0.580 0.56552 M1 -0.199445 0.307899 -0.648 0.52149 M2 -0.252552 0.310744 -0.813 0.42203 M3 -0.236224 0.295468 -0.799 0.42956 M4 -0.022686 0.330934 -0.069 0.94575 M5 -0.074551 0.344289 -0.217 0.82986 M6 -0.161684 0.302225 -0.535 0.59615 M7 0.160658 0.303948 0.529 0.60054 M8 0.042935 0.309732 0.139 0.89057 M9 0.906708 0.318013 2.851 0.00735 ** M10 -0.040783 0.322610 -0.126 0.90015 M11 -0.050514 0.349013 -0.145 0.88578 t -0.002753 0.004832 -0.570 0.57263 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3424 on 34 degrees of freedom Multiple R-squared: 0.9034, Adjusted R-squared: 0.8551 F-statistic: 18.7 on 17 and 34 DF, p-value: 1.683e-12 > 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.236818770 0.47363754 0.7631812 [2,] 0.163536222 0.32707244 0.8364638 [3,] 0.157996965 0.31599393 0.8420030 [4,] 0.121213101 0.24242620 0.8787869 [5,] 0.070886125 0.14177225 0.9291139 [6,] 0.121426373 0.24285275 0.8785736 [7,] 0.079069832 0.15813966 0.9209302 [8,] 0.044777221 0.08955444 0.9552228 [9,] 0.028842586 0.05768517 0.9711574 [10,] 0.011151380 0.02230276 0.9888486 [11,] 0.006911011 0.01382202 0.9930890 > postscript(file="/var/www/html/rcomp/tmp/10s101258733609.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/2orun1258733609.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/32inb1258733609.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/4rqik1258733609.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/5c5xs1258733609.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 = 52 Frequency = 1 1 2 3 4 5 6 -0.214545121 -0.059112476 0.224775251 -0.288758817 -0.160596187 -0.083451202 7 8 9 10 11 12 0.553918084 0.131043138 -0.018552088 0.099364859 0.067525388 -0.153388439 13 14 15 16 17 18 0.331337697 -0.080665556 0.043533749 -0.006369306 0.208171051 0.191626243 19 20 21 22 23 24 -0.153990514 -0.043718746 -0.193757105 0.276128067 0.077264051 -0.063077484 25 26 27 28 29 30 0.092764233 -0.460310139 0.296728722 0.044962804 0.155182373 -0.152994785 31 32 33 34 35 36 -0.212876352 0.019787519 -0.593530520 -0.729021798 0.024337199 0.006541298 37 38 39 40 41 42 -0.263655634 0.104620506 0.112138030 0.163354046 -0.202757237 0.044819744 43 44 45 46 47 48 -0.187051218 -0.107111911 0.805839713 0.353528872 -0.169126639 0.209924625 49 50 51 52 0.054098826 0.495467666 -0.677175751 0.086811273 > postscript(file="/var/www/html/rcomp/tmp/6blim1258733609.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 = 52 Frequency = 1 lag(myerror, k = 1) myerror 0 -0.214545121 NA 1 -0.059112476 -0.214545121 2 0.224775251 -0.059112476 3 -0.288758817 0.224775251 4 -0.160596187 -0.288758817 5 -0.083451202 -0.160596187 6 0.553918084 -0.083451202 7 0.131043138 0.553918084 8 -0.018552088 0.131043138 9 0.099364859 -0.018552088 10 0.067525388 0.099364859 11 -0.153388439 0.067525388 12 0.331337697 -0.153388439 13 -0.080665556 0.331337697 14 0.043533749 -0.080665556 15 -0.006369306 0.043533749 16 0.208171051 -0.006369306 17 0.191626243 0.208171051 18 -0.153990514 0.191626243 19 -0.043718746 -0.153990514 20 -0.193757105 -0.043718746 21 0.276128067 -0.193757105 22 0.077264051 0.276128067 23 -0.063077484 0.077264051 24 0.092764233 -0.063077484 25 -0.460310139 0.092764233 26 0.296728722 -0.460310139 27 0.044962804 0.296728722 28 0.155182373 0.044962804 29 -0.152994785 0.155182373 30 -0.212876352 -0.152994785 31 0.019787519 -0.212876352 32 -0.593530520 0.019787519 33 -0.729021798 -0.593530520 34 0.024337199 -0.729021798 35 0.006541298 0.024337199 36 -0.263655634 0.006541298 37 0.104620506 -0.263655634 38 0.112138030 0.104620506 39 0.163354046 0.112138030 40 -0.202757237 0.163354046 41 0.044819744 -0.202757237 42 -0.187051218 0.044819744 43 -0.107111911 -0.187051218 44 0.805839713 -0.107111911 45 0.353528872 0.805839713 46 -0.169126639 0.353528872 47 0.209924625 -0.169126639 48 0.054098826 0.209924625 49 0.495467666 0.054098826 50 -0.677175751 0.495467666 51 0.086811273 -0.677175751 52 NA 0.086811273 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -0.059112476 -0.214545121 [2,] 0.224775251 -0.059112476 [3,] -0.288758817 0.224775251 [4,] -0.160596187 -0.288758817 [5,] -0.083451202 -0.160596187 [6,] 0.553918084 -0.083451202 [7,] 0.131043138 0.553918084 [8,] -0.018552088 0.131043138 [9,] 0.099364859 -0.018552088 [10,] 0.067525388 0.099364859 [11,] -0.153388439 0.067525388 [12,] 0.331337697 -0.153388439 [13,] -0.080665556 0.331337697 [14,] 0.043533749 -0.080665556 [15,] -0.006369306 0.043533749 [16,] 0.208171051 -0.006369306 [17,] 0.191626243 0.208171051 [18,] -0.153990514 0.191626243 [19,] -0.043718746 -0.153990514 [20,] -0.193757105 -0.043718746 [21,] 0.276128067 -0.193757105 [22,] 0.077264051 0.276128067 [23,] -0.063077484 0.077264051 [24,] 0.092764233 -0.063077484 [25,] -0.460310139 0.092764233 [26,] 0.296728722 -0.460310139 [27,] 0.044962804 0.296728722 [28,] 0.155182373 0.044962804 [29,] -0.152994785 0.155182373 [30,] -0.212876352 -0.152994785 [31,] 0.019787519 -0.212876352 [32,] -0.593530520 0.019787519 [33,] -0.729021798 -0.593530520 [34,] 0.024337199 -0.729021798 [35,] 0.006541298 0.024337199 [36,] -0.263655634 0.006541298 [37,] 0.104620506 -0.263655634 [38,] 0.112138030 0.104620506 [39,] 0.163354046 0.112138030 [40,] -0.202757237 0.163354046 [41,] 0.044819744 -0.202757237 [42,] -0.187051218 0.044819744 [43,] -0.107111911 -0.187051218 [44,] 0.805839713 -0.107111911 [45,] 0.353528872 0.805839713 [46,] -0.169126639 0.353528872 [47,] 0.209924625 -0.169126639 [48,] 0.054098826 0.209924625 [49,] 0.495467666 0.054098826 [50,] -0.677175751 0.495467666 [51,] 0.086811273 -0.677175751 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -0.059112476 -0.214545121 2 0.224775251 -0.059112476 3 -0.288758817 0.224775251 4 -0.160596187 -0.288758817 5 -0.083451202 -0.160596187 6 0.553918084 -0.083451202 7 0.131043138 0.553918084 8 -0.018552088 0.131043138 9 0.099364859 -0.018552088 10 0.067525388 0.099364859 11 -0.153388439 0.067525388 12 0.331337697 -0.153388439 13 -0.080665556 0.331337697 14 0.043533749 -0.080665556 15 -0.006369306 0.043533749 16 0.208171051 -0.006369306 17 0.191626243 0.208171051 18 -0.153990514 0.191626243 19 -0.043718746 -0.153990514 20 -0.193757105 -0.043718746 21 0.276128067 -0.193757105 22 0.077264051 0.276128067 23 -0.063077484 0.077264051 24 0.092764233 -0.063077484 25 -0.460310139 0.092764233 26 0.296728722 -0.460310139 27 0.044962804 0.296728722 28 0.155182373 0.044962804 29 -0.152994785 0.155182373 30 -0.212876352 -0.152994785 31 0.019787519 -0.212876352 32 -0.593530520 0.019787519 33 -0.729021798 -0.593530520 34 0.024337199 -0.729021798 35 0.006541298 0.024337199 36 -0.263655634 0.006541298 37 0.104620506 -0.263655634 38 0.112138030 0.104620506 39 0.163354046 0.112138030 40 -0.202757237 0.163354046 41 0.044819744 -0.202757237 42 -0.187051218 0.044819744 43 -0.107111911 -0.187051218 44 0.805839713 -0.107111911 45 0.353528872 0.805839713 46 -0.169126639 0.353528872 47 0.209924625 -0.169126639 48 0.054098826 0.209924625 49 0.495467666 0.054098826 50 -0.677175751 0.495467666 51 0.086811273 -0.677175751 > 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/74be11258733609.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/8t5bm1258733609.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/9qo611258733609.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/10j90z1258733609.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/11btuk1258733609.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/12c27j1258733609.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/13ijp71258733609.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/14dryp1258733609.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/15s8kg1258733609.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/16y9ap1258733609.tab") + } > > system("convert tmp/10s101258733609.ps tmp/10s101258733609.png") > system("convert tmp/2orun1258733609.ps tmp/2orun1258733609.png") > system("convert tmp/32inb1258733609.ps tmp/32inb1258733609.png") > system("convert tmp/4rqik1258733609.ps tmp/4rqik1258733609.png") > system("convert tmp/5c5xs1258733609.ps tmp/5c5xs1258733609.png") > system("convert tmp/6blim1258733609.ps tmp/6blim1258733609.png") > system("convert tmp/74be11258733609.ps tmp/74be11258733609.png") > system("convert tmp/8t5bm1258733609.ps tmp/8t5bm1258733609.png") > system("convert tmp/9qo611258733609.ps tmp/9qo611258733609.png") > system("convert tmp/10j90z1258733609.ps tmp/10j90z1258733609.png") > > > proc.time() user system elapsed 2.337 1.569 5.318