R version 2.12.0 (2010-10-15) Copyright (C) 2010 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(-999.0 + ,38.6 + ,6654.000 + ,5712.000 + ,645.0 + ,3 + ,5 + ,3 + ,6.3 + ,4.5 + ,1.000 + ,6.600 + ,42.0 + ,3 + ,1 + ,3 + ,-999.0 + ,14.0 + ,3.385 + ,44.500 + ,60.0 + ,1 + ,1 + ,1 + ,-999.0 + ,-999.0 + ,0.920 + ,5.700 + ,25.0 + ,5 + ,2 + ,3 + ,2.1 + ,69.0 + ,2547.000 + ,4603.000 + ,624.0 + ,3 + ,5 + ,4 + ,9.1 + ,27.0 + ,10.550 + ,179.500 + ,180.0 + ,4 + ,4 + ,4 + ,15.8 + ,19.0 + ,0.023 + ,0.300 + ,35.0 + ,1 + ,1 + ,1 + ,5.2 + ,30.4 + ,160.000 + ,169.000 + ,392.0 + ,4 + ,5 + ,4 + ,10.9 + ,28.0 + ,3.300 + ,25.600 + ,63.0 + ,1 + ,2 + ,1 + ,8.3 + ,50.0 + ,52.160 + ,440.000 + ,230.0 + ,1 + ,1 + ,1 + ,11.0 + ,7.0 + ,0.425 + ,6.400 + ,112.0 + ,5 + ,4 + ,4 + ,3.2 + ,30.0 + ,465.000 + ,423.000 + ,281.0 + ,5 + ,5 + ,5 + ,7.6 + ,-999.0 + ,0.550 + ,2.400 + ,-999.0 + ,2 + ,1 + ,2 + ,-999.0 + ,40.0 + ,187.100 + ,419.000 + ,365.0 + ,5 + ,5 + ,5 + ,6.3 + ,0.075 + ,1.200 + ,42.0 + ,1 + ,1 + ,1 + ,8.6 + ,50.0 + ,3.000 + ,25.000 + ,28.0 + ,2 + ,2 + ,2 + ,6.6 + ,6.0 + ,0.785 + ,3.500 + ,42.0 + ,2 + ,2 + ,2 + ,9.5 + ,10.4 + ,0.200 + ,5.000 + ,120.0 + ,2 + ,2 + ,2 + ,4.8 + ,34.0 + ,1.410 + ,17.500 + ,-999.0 + ,1 + ,2 + ,1 + ,12.0 + ,7.0 + ,60.000 + ,81.000 + ,-999.0 + ,1 + ,1 + ,1 + ,-999.0 + ,28.0 + ,529.000 + ,680.000 + ,400.0 + ,5 + ,5 + ,5 + ,3.3 + ,20.0 + ,27.660 + ,115.000 + ,148.0 + ,5 + ,5 + ,5 + ,11.0 + ,3.9 + ,0.120 + ,1.000 + ,16.0 + ,3 + ,1 + ,2 + ,-999.0 + ,39.3 + ,207.000 + ,406.000 + ,252.0 + ,1 + ,4 + ,1 + ,4.7 + ,41.0 + ,85.000 + ,325.000 + ,310.0 + ,1 + ,3 + ,1 + ,-999.0 + ,16.2 + ,36.330 + ,119.500 + ,63.0 + ,1 + ,1 + ,1 + ,10.4 + ,9.0 + ,0.101 + ,4.000 + ,28.0 + ,5 + ,1 + ,3 + ,7.4 + ,7.6 + ,1.040 + ,5.500 + ,68.0 + ,5 + ,3 + ,4 + ,2.1 + ,46.0 + ,521.000 + ,655.000 + ,336.0 + ,5 + ,5 + ,5 + ,-999.0 + ,22.4 + ,100.000 + ,157.000 + ,100.0 + ,1 + ,1 + ,1 + ,-999.0 + ,16.3 + ,35.000 + ,56.000 + ,33.0 + ,3 + ,5 + ,4 + ,7.7 + ,2.6 + ,0.005 + ,0.140 + ,21.5 + ,5 + ,2 + ,4 + ,17.9 + ,24.0 + ,0.010 + ,0.250 + ,50.0 + ,1 + ,1 + ,1 + ,6.1 + ,100.0 + ,62.000 + ,1320.000 + ,267.0 + ,1 + ,1 + ,1 + ,8.2 + ,-999.0 + ,0.122 + ,3.000 + ,30.0 + ,2 + ,1 + ,1 + ,8.4 + ,-999.0 + ,1.350 + ,8.100 + ,45.0 + ,3 + ,1 + ,3 + ,11.9 + ,3.2 + ,0.023 + ,0.400 + ,19.0 + ,4 + ,1 + ,3 + ,10.8 + ,2.0 + ,0.048 + ,0.330 + ,30.0 + ,4 + ,1 + ,3 + ,13.8 + ,5.0 + ,1.700 + ,6.300 + ,12.0 + ,2 + ,1 + ,1 + ,14.3 + ,6.5 + ,3.500 + ,10.800 + ,120.0 + ,2 + ,1 + ,1 + ,-999.0 + ,23.6 + ,250.000 + ,490.000 + ,440.0 + ,5 + ,5 + ,5 + ,15.2 + ,12.0 + ,0.480 + ,15.500 + ,140.0 + ,2 + ,2 + ,2 + ,10.0 + ,20.2 + ,10.000 + ,115.000 + ,170.0 + ,4 + ,4 + ,4 + ,11.9 + ,13.0 + ,1.620 + ,11.400 + ,17.0 + ,2 + ,1 + ,2 + ,6.5 + ,27.0 + ,192.000 + ,180.000 + ,115.0 + ,4 + ,4 + ,4 + ,7.5 + ,18.0 + ,2.500 + ,12.100 + ,31.0 + ,5 + ,5 + ,5 + ,-999.0 + ,13.7 + ,4.288 + ,39.200 + ,63.0 + ,2 + ,2 + ,2 + ,10.6 + ,4.7 + ,0.280 + ,1.900 + ,21.0 + ,3 + ,1 + ,3 + ,7.4 + ,9.8 + ,4.235 + ,50.400 + ,52.0 + ,1 + ,1 + ,1 + ,8.4 + ,29.0 + ,6.800 + ,179.000 + ,164.0 + ,2 + ,3 + ,2 + ,5.7 + ,7.0 + ,0.750 + ,12.300 + ,225.0 + ,2 + ,2 + ,2 + ,4.9 + ,6.0 + ,3.600 + ,21.000 + ,225.0 + ,3 + ,2 + ,3 + ,-999.0 + ,17.0 + ,14.830 + ,98.200 + ,150.0 + ,5 + ,5 + ,5 + ,3.2 + ,20.0 + ,55.500 + ,175.000 + ,151.0 + ,5 + ,5 + ,5 + ,-999.0 + ,12.7 + ,1.400 + ,12.500 + ,90.0 + ,2 + ,2 + ,2 + ,8.1 + ,3.5 + ,0.060 + ,1.000 + ,-999.0 + ,3 + ,1 + ,1 + ,11.0 + ,4.5 + ,0.900 + ,2.600 + ,38.0 + ,2 + ,1 + ,2 + ,-999.0 + ,7.5 + ,2.000 + ,17.000 + ,200.0 + ,3 + ,1 + ,3 + ,13.2 + ,2.3 + ,0.104 + ,2.500 + ,46.0 + ,3 + ,2 + ,2 + ,9.7 + ,13.0 + ,4.050 + ,58.000 + ,210.0 + ,4 + ,3 + ,4 + ,12.8 + ,3.0 + ,3.500 + ,3.900 + ,14.0 + ,2 + ,1 + ,1 + ,4.9 + ,24.0 + ,4.190 + ,12.300 + ,60.0 + ,3 + ,1 + ,2) + ,dim=c(8 + ,62) + ,dimnames=list(c('SWS' + ,'L' + ,'wb' + ,'wbr' + ,'tg' + ,'P' + ,'S' + ,'D ') + ,1:62)) > y <- array(NA,dim=c(8,62),dimnames=list(c('SWS','L','wb','wbr','tg','P','S','D '),1:62)) > 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 = '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 SWS L wb wbr tg P S D\r 1 -999.0 38.600 6654.000 5712.0 645 3 5 3.0 2 6.3 4.500 1.000 6.6 42 3 1 3.0 3 -999.0 14.000 3.385 44.5 60 1 1 1.0 4 -999.0 -999.000 0.920 5.7 25 5 2 3.0 5 2.1 69.000 2547.000 4603.0 624 3 5 4.0 6 9.1 27.000 10.550 179.5 180 4 4 4.0 7 15.8 19.000 0.023 0.3 35 1 1 1.0 8 5.2 30.400 160.000 169.0 392 4 5 4.0 9 10.9 28.000 3.300 25.6 63 1 2 1.0 10 8.3 50.000 52.160 440.0 230 1 1 1.0 11 11.0 7.000 0.425 6.4 112 5 4 4.0 12 3.2 30.000 465.000 423.0 281 5 5 5.0 13 7.6 -999.000 0.550 2.4 -999 2 1 2.0 14 -999.0 40.000 187.100 419.0 365 5 5 5.0 15 6.3 0.075 1.200 42.0 1 1 1 8.6 16 50.0 3.000 25.000 28.0 2 2 2 6.6 17 6.0 0.785 3.500 42.0 2 2 2 9.5 18 10.4 0.200 5.000 120.0 2 2 2 4.8 19 34.0 1.410 17.500 -999.0 1 2 1 12.0 20 7.0 60.000 81.000 -999.0 1 1 1 -999.0 21 28.0 529.000 680.000 400.0 5 5 5 3.3 22 20.0 27.660 115.000 148.0 5 5 5 11.0 23 3.9 0.120 1.000 16.0 3 1 2 -999.0 24 39.3 207.000 406.000 252.0 1 4 1 4.7 25 41.0 85.000 325.000 310.0 1 3 1 -999.0 26 16.2 36.330 119.500 63.0 1 1 1 10.4 27 9.0 0.101 4.000 28.0 5 1 3 7.4 28 7.6 1.040 5.500 68.0 5 3 4 2.1 29 46.0 521.000 655.000 336.0 5 5 5 -999.0 30 22.4 100.000 157.000 100.0 1 1 1 -999.0 31 16.3 35.000 56.000 33.0 3 5 4 7.7 32 2.6 0.005 0.140 21.5 5 2 4 17.9 33 24.0 0.010 0.250 50.0 1 1 1 6.1 34 100.0 62.000 1320.000 267.0 1 1 1 8.2 35 -999.0 0.122 3.000 30.0 2 1 1 8.4 36 -999.0 1.350 8.100 45.0 3 1 3 11.9 37 3.2 0.023 0.400 19.0 4 1 3 10.8 38 2.0 0.048 0.330 30.0 4 1 3 13.8 39 5.0 1.700 6.300 12.0 2 1 1 14.3 40 6.5 3.500 10.800 120.0 2 1 1 -999.0 41 23.6 250.000 490.000 440.0 5 5 5 15.2 42 12.0 0.480 15.500 140.0 2 2 2 10.0 43 20.2 10.000 115.000 170.0 4 4 4 11.9 44 13.0 1.620 11.400 17.0 2 1 2 6.5 45 27.0 192.000 180.000 115.0 4 4 4 7.5 46 18.0 2.500 12.100 31.0 5 5 5 -999.0 47 13.7 4.288 39.200 63.0 2 2 2 10.6 48 4.7 0.280 1.900 21.0 3 1 3 7.4 49 9.8 4.235 50.400 52.0 1 1 1 8.4 50 29.0 6.800 179.000 164.0 2 3 2 5.7 51 7.0 0.750 12.300 225.0 2 2 2 4.9 52 6.0 3.600 21.000 225.0 3 2 3 -999.0 53 17.0 14.830 98.200 150.0 5 5 5 3.2 54 20.0 55.500 175.000 151.0 5 5 5 -999.0 55 12.7 1.400 12.500 90.0 2 2 2 8.1 56 3.5 0.060 1.000 -999.0 3 1 1 11.0 57 4.5 0.900 2.600 38.0 2 1 2 -999.0 58 7.5 2.000 17.000 200.0 3 1 3 13.2 59 2.3 0.104 2.500 46.0 3 2 2 9.7 60 13.0 4.050 58.000 210.0 4 3 4 12.8 61 3.0 3.500 3.900 14.0 2 1 1 4.9 62 24.0 4.190 12.300 60.0 3 1 2 -999.0 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) L wb wbr tg P -59.38472 0.66987 -0.15338 0.11272 -0.67714 1.93900 S `D\r` -2.07827 -0.04529 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -940.74 12.58 63.24 77.95 317.97 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -59.38472 76.26155 -0.779 0.43956 L 0.66987 0.20058 3.340 0.00153 ** wb -0.15338 0.09240 -1.660 0.10274 wbr 0.11272 0.09734 1.158 0.25193 tg -0.67714 0.26466 -2.559 0.01335 * P 1.93900 34.87417 0.056 0.95587 S -2.07827 39.64169 -0.052 0.95838 `D\r` -0.04529 0.09202 -0.492 0.62461 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 271.8 on 54 degrees of freedom Multiple R-squared: 0.2871, Adjusted R-squared: 0.1947 F-statistic: 3.106 on 7 and 54 DF, p-value: 0.00793 > 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.9885825 2.283503e-02 1.141752e-02 [2,] 0.9772781 4.544384e-02 2.272192e-02 [3,] 0.9903321 1.933570e-02 9.667851e-03 [4,] 0.9997814 4.372155e-04 2.186077e-04 [5,] 0.9994820 1.035986e-03 5.179930e-04 [6,] 0.9988317 2.336555e-03 1.168277e-03 [7,] 0.9974947 5.010615e-03 2.505308e-03 [8,] 0.9949739 1.005225e-02 5.026127e-03 [9,] 0.9915409 1.691821e-02 8.459107e-03 [10,] 0.9896695 2.066099e-02 1.033050e-02 [11,] 0.9864446 2.711088e-02 1.355544e-02 [12,] 0.9766600 4.668000e-02 2.334000e-02 [13,] 0.9634815 7.303691e-02 3.651846e-02 [14,] 0.9436012 1.127976e-01 5.639881e-02 [15,] 0.9150052 1.699896e-01 8.499481e-02 [16,] 0.8796686 2.406629e-01 1.203314e-01 [17,] 0.8356917 3.286166e-01 1.643083e-01 [18,] 0.7791208 4.417584e-01 2.208792e-01 [19,] 0.7327079 5.345842e-01 2.672921e-01 [20,] 0.6593896 6.812208e-01 3.406104e-01 [21,] 0.5834770 8.330461e-01 4.165230e-01 [22,] 0.5095824 9.808352e-01 4.904176e-01 [23,] 0.4320472 8.640944e-01 5.679528e-01 [24,] 0.3868188 7.736375e-01 6.131812e-01 [25,] 0.9821061 3.578786e-02 1.789393e-02 [26,] 1.0000000 1.324723e-28 6.623616e-29 [27,] 1.0000000 7.849893e-27 3.924947e-27 [28,] 1.0000000 4.546449e-25 2.273224e-25 [29,] 1.0000000 2.436430e-23 1.218215e-23 [30,] 1.0000000 1.173188e-21 5.865939e-22 [31,] 1.0000000 8.727041e-21 4.363520e-21 [32,] 1.0000000 4.682625e-19 2.341312e-19 [33,] 1.0000000 2.489244e-17 1.244622e-17 [34,] 1.0000000 8.734445e-16 4.367222e-16 [35,] 1.0000000 1.219173e-14 6.095865e-15 [36,] 1.0000000 4.997260e-13 2.498630e-13 [37,] 1.0000000 1.903925e-11 9.519625e-12 [38,] 1.0000000 1.004985e-09 5.024923e-10 [39,] 1.0000000 4.885821e-08 2.442910e-08 [40,] 0.9999992 1.623292e-06 8.116458e-07 [41,] 0.9999653 6.945045e-05 3.472523e-05 > postscript(file="/var/www/rcomp/tmp/13pt61292869839.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/www/rcomp/tmp/23pt61292869839.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/www/rcomp/tmp/3wgsr1292869839.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/www/rcomp/tmp/4wgsr1292869839.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/www/rcomp/tmp/5wgsr1292869839.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 = 62 Frequency = 1 1 2 3 4 5 -147.30571345 86.91671712 -912.67742455 -259.39333369 314.34712725 6 7 8 9 10 154.40625354 86.31138852 317.96621337 94.07152003 148.52066164 11 12 13 14 15 139.67826344 257.32650138 57.82451592 -736.86485531 62.29009405 16 17 18 19 20 109.98502930 62.72442143 58.74124435 207.15385159 106.80121766 21 22 23 24 25 -203.53627196 66.39188596 20.56304831 -10.90057520 10.04799065 26 27 28 29 30 63.76283831 73.79098793 65.44344886 -222.18709032 -16.81774264 31 32 33 34 35 58.10683074 68.21071954 78.87294293 291.40083837 -940.74457255 36 37 38 39 40 -937.48359276 67.98240653 65.65084369 65.00068840 7.92354455 41 42 43 44 45 -54.15475272 59.74518805 75.16621732 74.99793454 -23.97975128 46 47 48 49 50 32.91351340 71.23607507 68.48376253 69.41329039 92.74961729 51 52 53 54 55 44.26124133 -0.02007406 68.83089258 10.86897717 64.91881118 56 57 58 59 60 178.27545564 17.72894266 52.53308979 59.56253551 60.68037471 61 62 60.77569866 34.71009731 > postscript(file="/var/www/rcomp/tmp/67q9c1292869839.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 = 62 Frequency = 1 lag(myerror, k = 1) myerror 0 -147.30571345 NA 1 86.91671712 -147.30571345 2 -912.67742455 86.91671712 3 -259.39333369 -912.67742455 4 314.34712725 -259.39333369 5 154.40625354 314.34712725 6 86.31138852 154.40625354 7 317.96621337 86.31138852 8 94.07152003 317.96621337 9 148.52066164 94.07152003 10 139.67826344 148.52066164 11 257.32650138 139.67826344 12 57.82451592 257.32650138 13 -736.86485531 57.82451592 14 62.29009405 -736.86485531 15 109.98502930 62.29009405 16 62.72442143 109.98502930 17 58.74124435 62.72442143 18 207.15385159 58.74124435 19 106.80121766 207.15385159 20 -203.53627196 106.80121766 21 66.39188596 -203.53627196 22 20.56304831 66.39188596 23 -10.90057520 20.56304831 24 10.04799065 -10.90057520 25 63.76283831 10.04799065 26 73.79098793 63.76283831 27 65.44344886 73.79098793 28 -222.18709032 65.44344886 29 -16.81774264 -222.18709032 30 58.10683074 -16.81774264 31 68.21071954 58.10683074 32 78.87294293 68.21071954 33 291.40083837 78.87294293 34 -940.74457255 291.40083837 35 -937.48359276 -940.74457255 36 67.98240653 -937.48359276 37 65.65084369 67.98240653 38 65.00068840 65.65084369 39 7.92354455 65.00068840 40 -54.15475272 7.92354455 41 59.74518805 -54.15475272 42 75.16621732 59.74518805 43 74.99793454 75.16621732 44 -23.97975128 74.99793454 45 32.91351340 -23.97975128 46 71.23607507 32.91351340 47 68.48376253 71.23607507 48 69.41329039 68.48376253 49 92.74961729 69.41329039 50 44.26124133 92.74961729 51 -0.02007406 44.26124133 52 68.83089258 -0.02007406 53 10.86897717 68.83089258 54 64.91881118 10.86897717 55 178.27545564 64.91881118 56 17.72894266 178.27545564 57 52.53308979 17.72894266 58 59.56253551 52.53308979 59 60.68037471 59.56253551 60 60.77569866 60.68037471 61 34.71009731 60.77569866 62 NA 34.71009731 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] 86.91671712 -147.30571345 [2,] -912.67742455 86.91671712 [3,] -259.39333369 -912.67742455 [4,] 314.34712725 -259.39333369 [5,] 154.40625354 314.34712725 [6,] 86.31138852 154.40625354 [7,] 317.96621337 86.31138852 [8,] 94.07152003 317.96621337 [9,] 148.52066164 94.07152003 [10,] 139.67826344 148.52066164 [11,] 257.32650138 139.67826344 [12,] 57.82451592 257.32650138 [13,] -736.86485531 57.82451592 [14,] 62.29009405 -736.86485531 [15,] 109.98502930 62.29009405 [16,] 62.72442143 109.98502930 [17,] 58.74124435 62.72442143 [18,] 207.15385159 58.74124435 [19,] 106.80121766 207.15385159 [20,] -203.53627196 106.80121766 [21,] 66.39188596 -203.53627196 [22,] 20.56304831 66.39188596 [23,] -10.90057520 20.56304831 [24,] 10.04799065 -10.90057520 [25,] 63.76283831 10.04799065 [26,] 73.79098793 63.76283831 [27,] 65.44344886 73.79098793 [28,] -222.18709032 65.44344886 [29,] -16.81774264 -222.18709032 [30,] 58.10683074 -16.81774264 [31,] 68.21071954 58.10683074 [32,] 78.87294293 68.21071954 [33,] 291.40083837 78.87294293 [34,] -940.74457255 291.40083837 [35,] -937.48359276 -940.74457255 [36,] 67.98240653 -937.48359276 [37,] 65.65084369 67.98240653 [38,] 65.00068840 65.65084369 [39,] 7.92354455 65.00068840 [40,] -54.15475272 7.92354455 [41,] 59.74518805 -54.15475272 [42,] 75.16621732 59.74518805 [43,] 74.99793454 75.16621732 [44,] -23.97975128 74.99793454 [45,] 32.91351340 -23.97975128 [46,] 71.23607507 32.91351340 [47,] 68.48376253 71.23607507 [48,] 69.41329039 68.48376253 [49,] 92.74961729 69.41329039 [50,] 44.26124133 92.74961729 [51,] -0.02007406 44.26124133 [52,] 68.83089258 -0.02007406 [53,] 10.86897717 68.83089258 [54,] 64.91881118 10.86897717 [55,] 178.27545564 64.91881118 [56,] 17.72894266 178.27545564 [57,] 52.53308979 17.72894266 [58,] 59.56253551 52.53308979 [59,] 60.68037471 59.56253551 [60,] 60.77569866 60.68037471 [61,] 34.71009731 60.77569866 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 86.91671712 -147.30571345 2 -912.67742455 86.91671712 3 -259.39333369 -912.67742455 4 314.34712725 -259.39333369 5 154.40625354 314.34712725 6 86.31138852 154.40625354 7 317.96621337 86.31138852 8 94.07152003 317.96621337 9 148.52066164 94.07152003 10 139.67826344 148.52066164 11 257.32650138 139.67826344 12 57.82451592 257.32650138 13 -736.86485531 57.82451592 14 62.29009405 -736.86485531 15 109.98502930 62.29009405 16 62.72442143 109.98502930 17 58.74124435 62.72442143 18 207.15385159 58.74124435 19 106.80121766 207.15385159 20 -203.53627196 106.80121766 21 66.39188596 -203.53627196 22 20.56304831 66.39188596 23 -10.90057520 20.56304831 24 10.04799065 -10.90057520 25 63.76283831 10.04799065 26 73.79098793 63.76283831 27 65.44344886 73.79098793 28 -222.18709032 65.44344886 29 -16.81774264 -222.18709032 30 58.10683074 -16.81774264 31 68.21071954 58.10683074 32 78.87294293 68.21071954 33 291.40083837 78.87294293 34 -940.74457255 291.40083837 35 -937.48359276 -940.74457255 36 67.98240653 -937.48359276 37 65.65084369 67.98240653 38 65.00068840 65.65084369 39 7.92354455 65.00068840 40 -54.15475272 7.92354455 41 59.74518805 -54.15475272 42 75.16621732 59.74518805 43 74.99793454 75.16621732 44 -23.97975128 74.99793454 45 32.91351340 -23.97975128 46 71.23607507 32.91351340 47 68.48376253 71.23607507 48 69.41329039 68.48376253 49 92.74961729 69.41329039 50 44.26124133 92.74961729 51 -0.02007406 44.26124133 52 68.83089258 -0.02007406 53 10.86897717 68.83089258 54 64.91881118 10.86897717 55 178.27545564 64.91881118 56 17.72894266 178.27545564 57 52.53308979 17.72894266 58 59.56253551 52.53308979 59 60.68037471 59.56253551 60 60.77569866 60.68037471 61 34.71009731 60.77569866 > 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/rcomp/tmp/77q9c1292869839.ps",horizontal=F,onefile=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/rcomp/tmp/8hz9x1292869839.ps",horizontal=F,onefile=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/rcomp/tmp/9hz9x1292869839.ps",horizontal=F,onefile=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/rcomp/tmp/10s88z1292869839.ps",horizontal=F,onefile=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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/11vqon1292869839.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/rcomp/tmp/12oi5q1292869839.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/rcomp/tmp/13d1321292869839.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/rcomp/tmp/14ns251292869839.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/rcomp/tmp/15rtib1292869839.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/rcomp/tmp/16nkyj1292869839.tab") + } > > try(system("convert tmp/13pt61292869839.ps tmp/13pt61292869839.png",intern=TRUE)) character(0) > try(system("convert tmp/23pt61292869839.ps tmp/23pt61292869839.png",intern=TRUE)) character(0) > try(system("convert tmp/3wgsr1292869839.ps tmp/3wgsr1292869839.png",intern=TRUE)) character(0) > try(system("convert tmp/4wgsr1292869839.ps tmp/4wgsr1292869839.png",intern=TRUE)) character(0) > try(system("convert tmp/5wgsr1292869839.ps tmp/5wgsr1292869839.png",intern=TRUE)) character(0) > try(system("convert tmp/67q9c1292869839.ps tmp/67q9c1292869839.png",intern=TRUE)) character(0) > try(system("convert tmp/77q9c1292869839.ps tmp/77q9c1292869839.png",intern=TRUE)) character(0) > try(system("convert tmp/8hz9x1292869839.ps tmp/8hz9x1292869839.png",intern=TRUE)) character(0) > try(system("convert tmp/9hz9x1292869839.ps tmp/9hz9x1292869839.png",intern=TRUE)) character(0) > try(system("convert tmp/10s88z1292869839.ps tmp/10s88z1292869839.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.170 0.860 4.018