R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-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(3.3,76.1,3.7,63.1,6.8,70.0,9.8,51.3,13.6,66.5,16.2,71.8,18.4,73.5,18.0,79.3,14.9,68.9,11.1,74.5,6.8,76.4,3.9,81.0),dim=c(2,12),dimnames=list(c('temperatuur','neerslag'),1:12)) > y <- array(NA,dim=c(2,12),dimnames=list(c('temperatuur','neerslag'),1:12)) > 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 > 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 temperatuur neerslag M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t 1 3.3 76.1 1 0 0 0 0 0 0 0 0 0 0 1 2 3.7 63.1 0 1 0 0 0 0 0 0 0 0 0 2 3 6.8 70.0 0 0 1 0 0 0 0 0 0 0 0 3 4 9.8 51.3 0 0 0 1 0 0 0 0 0 0 0 4 5 13.6 66.5 0 0 0 0 1 0 0 0 0 0 0 5 6 16.2 71.8 0 0 0 0 0 1 0 0 0 0 0 6 7 18.4 73.5 0 0 0 0 0 0 1 0 0 0 0 7 8 18.0 79.3 0 0 0 0 0 0 0 1 0 0 0 8 9 14.9 68.9 0 0 0 0 0 0 0 0 1 0 0 9 10 11.1 74.5 0 0 0 0 0 0 0 0 0 1 0 10 11 6.8 76.4 0 0 0 0 0 0 0 0 0 0 1 11 12 3.9 81.0 0 0 0 0 0 0 0 0 0 0 0 12 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) neerslag M1 M2 M3 M4 54.9652 -0.6304 -3.6891 -11.4848 -4.0348 -12.8239 M5 M6 M7 M8 M9 M10 0.5587 6.5000 9.7717 13.0283 3.3717 3.1022 M11 t NA NA > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: ALL 12 residuals are 0: no residual degrees of freedom! Coefficients: (2 not defined because of singularities) Estimate Std. Error t value Pr(>|t|) (Intercept) 54.9652 NA NA NA neerslag -0.6304 NA NA NA M1 -3.6891 NA NA NA M2 -11.4848 NA NA NA M3 -4.0348 NA NA NA M4 -12.8239 NA NA NA M5 0.5587 NA NA NA M6 6.5000 NA NA NA M7 9.7717 NA NA NA M8 13.0283 NA NA NA M9 3.3717 NA NA NA M10 3.1022 NA NA NA M11 NA NA NA NA t NA NA NA NA Residual standard error: NaN on 0 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: NaN F-statistic: NaN on 11 and 0 DF, p-value: NA > 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 + } > postscript(file="/var/wessaorg/rcomp/tmp/1dsdr1321956535.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/28gld1321956535.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/3jdwm1321956535.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/4gtxi1321956535.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/5efb81321956535.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 = 12 Frequency = 1 1 2 3 4 5 6 7 8 9 10 11 12 0 0 0 0 0 0 0 0 0 0 0 0 > postscript(file="/var/wessaorg/rcomp/tmp/63n6v1321956535.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 = 12 Frequency = 1 lag(myerror, k = 1) myerror 0 0 NA 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0 10 0 0 11 0 0 12 NA 0 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] 0 0 [2,] 0 0 [3,] 0 0 [4,] 0 0 [5,] 0 0 [6,] 0 0 [7,] 0 0 [8,] 0 0 [9,] 0 0 [10,] 0 0 [11,] 0 0 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0 10 0 0 11 0 0 > 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)) Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : 'a' and 'b' must be finite Calls: abline -> int_abline Execution halted