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 + ,1 + ,4 + ,0 + ,2 + ,1 + ,1 + ,0 + ,0 + ,2 + ,0 + ,1 + ,4 + ,1 + ,1.5 + ,0 + ,0 + ,0 + ,0 + ,0 + ,1 + ,1 + ,0 + ,1 + ,1 + ,1 + ,1 + ,0 + ,1 + ,2 + ,1 + ,1 + ,0 + ,1 + ,2 + ,0 + ,1 + ,0 + ,1 + ,1 + ,0 + ,1 + ,4 + ,1 + ,2 + ,1 + ,1 + ,1 + ,0 + ,2 + ,0 + ,0 + ,4 + ,0 + ,2 + ,0 + ,1 + ,0 + ,1 + ,0 + ,0 + ,1 + ,2 + ,1 + ,0 + ,0 + ,1 + ,0 + ,0 + ,2 + ,0 + ,0 + ,0 + ,NA + ,NA + ,1 + ,1 + ,0 + ,1 + ,2 + ,1 + ,1 + ,1 + ,0 + ,2 + ,1 + ,1 + ,0 + ,1 + ,0.5 + ,0 + ,1 + ,0 + ,1 + ,2 + ,0 + ,0 + ,2 + ,1 + ,0 + ,1 + ,1 + ,2 + ,1 + ,2 + ,1 + ,1 + ,1 + ,0 + ,0 + ,0 + ,0 + ,2 + ,NA + ,NA + ,1 + ,0 + ,0 + ,NA + ,NA + ,1 + ,1 + ,3 + ,1 + ,2 + ,1 + ,0 + ,0 + ,1 + ,0 + ,1 + ,1 + ,0 + ,NA + ,NA + ,0 + ,0 + ,0 + ,NA + ,NA + ,0 + ,0 + ,1 + ,0 + ,2 + ,1 + ,1 + ,0 + ,1 + ,1 + ,1 + ,0 + ,0 + ,0 + ,0.5 + ,1 + ,1 + ,4 + ,0 + ,2 + ,0 + ,0 + ,0 + ,1 + ,0.5 + ,0 + ,0 + ,1 + ,NA + ,NA + ,0 + ,0 + ,0 + ,1 + ,0.5 + ,1 + ,1 + ,0 + ,NA + ,NA + ,1 + ,1 + ,4 + ,0 + ,2 + ,0 + ,1 + ,1 + ,1 + ,0 + ,0 + ,1 + ,0 + ,1 + ,1 + ,1 + ,1 + ,4 + ,1 + ,2 + ,1 + ,1 + ,0 + ,1 + ,1 + ,1 + ,1 + ,4 + ,1 + ,2 + ,1 + ,1 + ,0 + ,0 + ,0 + ,1 + ,1 + ,0 + ,1 + ,0.5 + ,0 + ,0 + ,0 + ,1 + ,0 + ,0 + ,1 + ,4 + ,1 + ,2 + ,0 + ,1 + ,0 + ,0 + ,0 + ,1 + ,1 + ,0 + ,0 + ,1 + ,1 + ,1 + ,4 + ,1 + ,2 + ,0 + ,0 + ,4 + ,0 + ,0.5 + ,0 + ,1 + ,0 + ,1 + ,2 + ,1 + ,1 + ,1 + ,1 + ,2 + ,0 + ,1 + ,0 + ,1 + ,2 + ,0 + ,0 + ,4 + ,NA + ,NA + ,0 + ,1 + ,0 + ,0 + ,0 + ,0 + ,1 + ,2 + ,1 + ,0 + ,0 + ,1 + ,0 + ,1 + ,0.5 + ,0 + ,1 + ,4 + ,NA + ,NA + ,0 + ,0 + ,4 + ,0 + ,2 + ,0 + ,0 + ,0 + ,NA + ,NA + ,0 + ,1 + ,0 + ,1 + ,0 + ,1 + ,1 + ,4 + ,1 + ,2 + ,1 + ,1 + ,0 + ,1 + ,1 + ,1 + ,0 + ,0 + ,1 + ,0 + ,0 + ,0 + ,2 + ,1 + ,2 + ,0 + ,1 + ,0 + ,0 + ,1 + ,0 + ,1 + ,0 + ,1 + ,2 + ,0 + ,0 + ,0 + ,0 + ,0 + ,1 + ,1 + ,4 + ,1 + ,1 + ,1 + ,1 + ,4 + ,1 + ,2 + ,0 + ,1 + ,2 + ,0 + ,0 + ,0 + ,1 + ,0 + ,0 + ,0 + ,0 + ,1 + ,0 + ,0 + ,0 + ,0 + ,1 + ,4 + ,0 + ,0 + ,1 + ,1 + ,0 + ,1 + ,2 + ,1 + ,0 + ,0 + ,1 + ,2 + ,0 + ,0 + ,1 + ,1 + ,2 + ,1 + ,1 + ,2 + ,1 + ,2 + ,1 + ,0 + ,0 + ,1 + ,2 + ,1 + ,1 + ,2 + ,1 + ,2 + ,0 + ,0 + ,0 + ,1 + ,2 + ,0 + ,0 + ,4 + ,1 + ,2 + ,0 + ,0 + ,4 + ,1 + ,2 + ,1 + ,0 + ,0 + ,1 + ,2 + ,0 + ,0 + ,0 + ,NA + ,NA + ,0 + ,0 + ,4 + ,1 + ,2 + ,1 + ,0 + ,0 + ,NA + ,NA + ,1 + ,1 + ,4 + ,1 + ,2 + ,0 + ,0 + ,2 + ,1 + ,2 + ,0 + ,0 + ,2 + ,NA + ,NA + ,1 + ,1 + ,0 + ,0 + ,0 + ,1 + ,1 + ,0 + ,1 + ,2 + ,1 + ,1 + ,4 + ,NA + ,NA + ,0 + ,1 + ,0 + ,1 + ,2 + ,1 + ,1 + ,0 + ,1 + ,2 + ,1 + ,1 + ,0 + ,1 + ,2 + ,1 + ,1 + ,4 + ,1 + ,2 + ,1 + ,1 + ,4 + ,1 + ,2 + ,0 + ,0 + ,0 + ,NA + ,NA + ,0 + ,0 + ,0 + ,0 + ,0 + ,1 + ,1 + ,2 + ,0 + ,0 + ,0 + ,0 + ,1 + ,1 + ,2 + ,0 + ,0 + ,0 + ,0 + ,0 + ,0 + ,0 + ,2 + ,1 + ,2 + ,0 + ,1 + ,1 + ,0 + ,0) + ,dim=c(5 + ,105) + ,dimnames=list(c('pre' + ,'post1' + ,'post2' + ,'post3' + ,'post4') + ,1:105)) > y <- array(NA,dim=c(5,105),dimnames=list(c('pre','post1','post2','post3','post4'),1:105)) > 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 = '5' > par3 <- 'No Linear Trend' > par2 <- 'Do not include Seasonal Dummies' > par1 <- '5' > #'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, 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 post4 pre post1 post2 post3 1 2.0 1 1 4 0 2 2.0 1 1 0 0 3 1.5 0 1 4 1 4 0.0 0 0 0 0 5 1.0 1 1 0 1 6 2.0 1 1 0 1 7 2.0 1 1 0 1 8 1.0 0 1 0 1 9 2.0 0 1 4 1 10 2.0 1 1 1 0 11 2.0 0 0 4 0 12 0.0 0 1 0 1 13 0.0 0 1 2 1 14 2.0 0 1 0 0 15 NA 0 0 0 NA 16 2.0 1 1 0 1 17 2.0 1 1 1 0 18 0.5 1 1 0 1 19 2.0 0 1 0 1 20 0.0 0 0 2 1 21 2.0 1 1 2 1 22 0.0 1 1 1 0 23 NA 0 0 2 NA 24 NA 1 0 0 NA 25 2.0 1 1 3 1 26 0.0 1 0 0 1 27 NA 1 1 0 NA 28 NA 0 0 0 NA 29 2.0 0 0 1 0 30 1.0 1 1 0 1 31 0.5 1 0 0 0 32 2.0 1 1 4 0 33 0.5 0 0 0 1 34 NA 0 0 1 NA 35 0.5 0 0 0 1 36 NA 1 1 0 NA 37 2.0 1 1 4 0 38 0.0 0 1 1 1 39 1.0 0 1 0 1 40 2.0 1 1 4 1 41 1.0 1 1 0 1 42 2.0 1 1 4 1 43 0.0 1 1 0 0 44 0.5 1 1 0 1 45 0.0 0 0 0 1 46 2.0 0 1 4 1 47 0.0 0 1 0 0 48 1.0 1 1 0 0 49 2.0 1 1 4 1 50 0.5 0 0 4 0 51 2.0 0 1 0 1 52 2.0 1 1 1 1 53 2.0 0 1 0 1 54 NA 0 0 4 NA 55 0.0 0 1 0 0 56 0.0 0 1 2 1 57 0.5 0 1 0 1 58 NA 0 1 4 NA 59 2.0 0 0 4 0 60 NA 0 0 0 NA 61 0.0 0 1 0 1 62 2.0 1 1 4 1 63 1.0 1 1 0 1 64 0.0 1 0 0 1 65 2.0 0 0 2 1 66 1.0 0 1 0 0 67 2.0 0 1 0 1 68 0.0 0 0 0 0 69 1.0 1 1 4 1 70 2.0 1 1 4 1 71 0.0 0 1 2 0 72 0.0 0 1 0 0 73 0.0 0 1 0 0 74 0.0 0 1 4 0 75 2.0 1 1 0 1 76 2.0 1 0 0 1 77 2.0 0 0 1 1 78 2.0 1 1 2 1 79 2.0 1 0 0 1 80 2.0 1 1 2 1 81 2.0 0 0 0 1 82 2.0 0 0 4 1 83 2.0 0 0 4 1 84 2.0 1 0 0 1 85 NA 0 0 0 NA 86 2.0 0 0 4 1 87 NA 1 0 0 NA 88 2.0 1 1 4 1 89 2.0 0 0 2 1 90 NA 0 0 2 NA 91 0.0 1 1 0 0 92 2.0 1 1 0 1 93 NA 1 1 4 NA 94 2.0 0 1 0 1 95 2.0 1 1 0 1 96 2.0 1 1 0 1 97 2.0 1 1 4 1 98 2.0 1 1 4 1 99 NA 0 0 0 NA 100 0.0 0 0 0 0 101 0.0 1 1 2 0 102 2.0 0 0 1 1 103 0.0 0 0 0 0 104 2.0 0 0 2 1 105 0.0 0 1 1 0 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) pre post1 post2 post3 0.4638 0.4263 -0.1021 0.1751 0.5997 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -1.4897 -0.4637 -0.0246 0.6059 1.6383 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.46375 0.19655 2.360 0.020591 * pre 0.42625 0.17240 2.472 0.015410 * post1 -0.10206 0.18695 -0.546 0.586537 post2 0.17506 0.04836 3.620 0.000499 *** post3 0.59966 0.17277 3.471 0.000818 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.7614 on 85 degrees of freedom (15 observations deleted due to missingness) Multiple R-squared: 0.2923, Adjusted R-squared: 0.259 F-statistic: 8.776 on 4 and 85 DF, p-value: 5.582e-06 > 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 + } Error in if (gqarr[mypoint - kp3 + 1, 2] < 0.01) numsignificant1 <- numsignificant1 + : missing value where TRUE/FALSE needed Execution halted