R version 2.8.0 (2008-10-20) Copyright (C) 2008 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. Natural language support but running in an English locale 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 + ,-999.0 + ,38.6 + ,6654.000 + ,5712.000 + ,645.0 + ,3 + ,5 + ,3 + ,6.3 + ,2.0 + ,4.5 + ,1.000 + ,6.600 + ,42.0 + ,3 + ,1 + ,3 + ,-999.0 + ,-999.0 + ,14.0 + ,3.385 + ,44.500 + ,60.0 + ,1 + ,1 + ,1 + ,-999.0 + ,-999.0 + ,-999.0 + ,0.920 + ,5.700 + ,25.0 + ,5 + ,2 + ,3 + ,2.1 + ,1.8 + ,69.0 + ,2547.000 + ,4603.000 + ,624.0 + ,3 + ,5 + ,4 + ,9.1 + ,0.7 + ,27.0 + ,10.550 + ,179.500 + ,180.0 + ,4 + ,4 + ,4 + ,15.8 + ,3.9 + ,19.0 + ,0.023 + ,0.300 + ,35.0 + ,1 + ,1 + ,1 + ,5.2 + ,1.0 + ,30.4 + ,160.000 + ,169.000 + ,392.0 + ,4 + ,5 + ,4 + ,10.9 + ,36.0 + ,28.0 + ,3.300 + ,25.600 + ,63.0 + ,1 + ,2 + ,1 + ,8.3 + ,1.4 + ,50.0 + ,52.160 + ,440.000 + ,230.0 + ,1 + ,1 + ,1 + ,11.0 + ,1.5 + ,7.0 + ,0.425 + ,6.400 + ,112.0 + ,5 + ,4 + ,4 + ,3.2 + ,0.7 + ,30.0 + ,465.000 + ,423.000 + ,281.0 + ,5 + ,5 + ,5 + ,7.6 + ,2.7 + ,-999.0 + ,0.550 + ,2.400 + ,-999.0 + ,2 + ,1 + ,2 + ,-999.0 + ,-999.0 + ,40.0 + ,187.100 + ,419.000 + ,365.0 + ,5 + ,5 + ,5 + ,6.3 + ,2.1 + ,3.5 + ,0.075 + ,1.200 + ,42.0 + ,1 + ,1 + ,1 + ,8.6 + ,0.0 + ,50.0 + ,3.000 + ,25.000 + ,28.0 + ,2 + ,2 + ,2 + ,6.6 + ,4.1 + ,6.0 + ,0.785 + ,3.500 + ,42.0 + ,2 + ,2 + ,2 + ,9.5 + ,1.2 + ,10.4 + ,0.200 + ,5.000 + ,120.0 + ,2 + ,2 + ,2 + ,4.8 + ,1.3 + ,34.0 + ,1.410 + ,17.500 + ,-999.0 + ,1 + ,2 + ,1 + ,12.0 + ,6.1 + ,7.0 + ,60.000 + ,81.000 + ,-999.0 + ,1 + ,1 + ,1 + ,-999.0 + ,0.3 + ,28.0 + ,529.000 + ,680.000 + ,400.0 + ,5 + ,5 + ,5 + ,3.3 + ,0.5 + ,20.0 + ,27.660 + ,115.000 + ,148.0 + ,5 + ,5 + ,5 + ,11.0 + ,3.4 + ,3.9 + ,0.120 + ,1.000 + ,16.0 + ,3 + ,1 + ,2 + ,-999.0 + ,-999.0 + ,39.3 + ,207.000 + ,406.000 + ,252.0 + ,1 + ,4 + ,1 + ,4.7 + ,1.5 + ,41.0 + ,85.000 + ,325.000 + ,310.0 + ,1 + ,3 + ,1 + ,-999.0 + ,-999.0 + ,16.2 + ,36.330 + ,119.500 + ,63.0 + ,1 + ,1 + ,1 + ,10.4 + ,3.4 + ,9.0 + ,0.101 + ,4.000 + ,28.0 + ,5 + ,1 + ,3 + ,7.4 + ,0.8 + ,7.6 + ,1.040 + ,5.500 + ,68.0 + ,5 + ,3 + ,4 + ,2.1 + ,0.8 + ,46.0 + ,521.000 + ,655.000 + ,336.0 + ,5 + ,5 + ,5 + ,-999.0 + ,-999.0 + ,22.4 + ,100.000 + ,157.000 + ,100.0 + ,1 + ,1 + ,1 + ,-999.0 + ,-999.0 + ,16.3 + ,35.000 + ,56.000 + ,33.0 + ,3 + ,5 + ,4 + ,7.7 + ,1.4 + ,2.6 + ,0.005 + ,0.140 + ,21.5 + ,5 + ,2 + ,4 + ,17.9 + ,2.0 + ,24.0 + ,0.010 + ,0.250 + ,50.0 + ,1 + ,1 + ,1 + ,6.1 + ,1.9 + ,100.0 + ,62.000 + ,1320.000 + ,267.0 + ,1 + ,1 + ,1 + ,8.2 + ,2.4 + ,-999.0 + ,0.122 + ,3.000 + ,30.0 + ,2 + ,1 + ,1 + ,8.4 + ,2.8 + ,-999.0 + ,1.350 + ,8.100 + ,45.0 + ,3 + ,1 + ,3 + ,11.9 + ,1.3 + ,3.2 + ,0.023 + ,0.400 + ,19.0 + ,4 + ,1 + ,3 + ,10.8 + ,2.0 + ,2.0 + ,0.048 + ,0.330 + ,30.0 + ,4 + ,1 + ,3 + ,13.8 + ,5.6 + ,5.0 + ,1.700 + ,6.300 + ,12.0 + ,2 + ,1 + ,1 + ,14.3 + ,3.1 + ,6.5 + ,3.500 + ,10.800 + ,120.0 + ,2 + ,1 + ,1 + ,-999.0 + ,1.0 + ,23.6 + ,250.000 + ,490.000 + ,440.0 + ,5 + ,5 + ,5 + ,15.2 + ,1.8 + ,12.0 + ,0.480 + ,15.500 + ,140.0 + ,2 + ,2 + ,2 + ,10.0 + ,0.9 + ,20.2 + ,10.000 + ,115.000 + ,170.0 + ,4 + ,4 + ,4 + ,11.9 + ,1.8 + ,13.0 + ,1.620 + ,11.400 + ,17.0 + ,2 + ,1 + ,2 + ,6.5 + ,1.9 + ,27.0 + ,192.000 + ,180.000 + ,115.0 + ,4 + ,4 + ,4 + ,7.5 + ,0.9 + ,18.0 + ,2.500 + ,12.100 + ,31.0 + ,5 + ,5 + ,5 + ,-999.0 + ,-999.0 + ,13.7 + ,4.288 + ,39.200 + ,63.0 + ,2 + ,2 + ,2 + ,10.6 + ,2.6 + ,4.7 + ,0.280 + ,1.900 + ,21.0 + ,3 + ,1 + ,3 + ,7.4 + ,2.4 + ,9.8 + ,4.235 + ,50.400 + ,52.0 + ,1 + ,1 + ,1 + ,8.4 + ,1.2 + ,29.0 + ,6.800 + ,179.000 + ,164.0 + ,2 + ,3 + ,2 + ,5.7 + ,0.9 + ,7.0 + ,0.750 + ,12.300 + ,225.0 + ,2 + ,2 + ,2 + ,4.9 + ,0.5 + ,6.0 + ,3.600 + ,21.000 + ,225.0 + ,3 + ,2 + ,3 + ,-999.0 + ,-999.0 + ,17.0 + ,14.830 + ,98.200 + ,150.0 + ,5 + ,5 + ,5 + ,3.2 + ,0.6 + ,20.0 + ,55.500 + ,175.000 + ,151.0 + ,5 + ,5 + ,5 + ,-999.0 + ,-999.0 + ,12.7 + ,1.400 + ,12.500 + ,90.0 + ,2 + ,2 + ,2 + ,8.1 + ,2.2 + ,3.5 + ,0.060 + ,1.000 + ,-999.0 + ,3 + ,1 + ,1 + ,11.0 + ,2.3 + ,4.5 + ,0.900 + ,2.600 + ,38.0 + ,2 + ,1 + ,2 + ,-999.0 + ,0.5 + ,7.5 + ,2.000 + ,17.000 + ,200.0 + ,3 + ,1 + ,3 + ,13.2 + ,-999.0 + ,2.3 + ,0.104 + ,2.500 + ,46.0 + ,3 + ,2 + ,2 + ,9.7 + ,0.6 + ,13.0 + ,4.050 + ,58.000 + ,210.0 + ,4 + ,3 + ,4 + ,12.8 + ,6.6 + ,3.0 + ,3.500 + ,3.900 + ,14.0 + ,2 + ,1 + ,1 + ,4.9 + ,2.6 + ,24.0 + ,4.190 + ,12.300 + ,60.0 + ,3 + ,1 + ,2) + ,dim=c(9 + ,62) + ,dimnames=list(c('SWS' + ,'PS' + ,'L' + ,'wb' + ,'wbr' + ,'tg' + ,'P' + ,'S' + ,'D ') + ,1:62)) > y <- array(NA,dim=c(9,62),dimnames=list(c('SWS','PS','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 = '2' > #'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 PS SWS L wb wbr tg P S D\r 1 -999.0 -999.0 38.6 6654.000 5712.00 645.0 3 5 3 2 2.0 6.3 4.5 1.000 6.60 42.0 3 1 3 3 -999.0 -999.0 14.0 3.385 44.50 60.0 1 1 1 4 -999.0 -999.0 -999.0 0.920 5.70 25.0 5 2 3 5 1.8 2.1 69.0 2547.000 4603.00 624.0 3 5 4 6 0.7 9.1 27.0 10.550 179.50 180.0 4 4 4 7 3.9 15.8 19.0 0.023 0.30 35.0 1 1 1 8 1.0 5.2 30.4 160.000 169.00 392.0 4 5 4 9 36.0 10.9 28.0 3.300 25.60 63.0 1 2 1 10 1.4 8.3 50.0 52.160 440.00 230.0 1 1 1 11 1.5 11.0 7.0 0.425 6.40 112.0 5 4 4 12 0.7 3.2 30.0 465.000 423.00 281.0 5 5 5 13 2.7 7.6 -999.0 0.550 2.40 -999.0 2 1 2 14 -999.0 -999.0 40.0 187.100 419.00 365.0 5 5 5 15 2.1 6.3 3.5 0.075 1.20 42.0 1 1 1 16 0.0 8.6 50.0 3.000 25.00 28.0 2 2 2 17 4.1 6.6 6.0 0.785 3.50 42.0 2 2 2 18 1.2 9.5 10.4 0.200 5.00 120.0 2 2 2 19 1.3 4.8 34.0 1.410 17.50 -999.0 1 2 1 20 6.1 12.0 7.0 60.000 81.00 -999.0 1 1 1 21 0.3 -999.0 28.0 529.000 680.00 400.0 5 5 5 22 0.5 3.3 20.0 27.660 115.00 148.0 5 5 5 23 3.4 11.0 3.9 0.120 1.00 16.0 3 1 2 24 -999.0 -999.0 39.3 207.000 406.00 252.0 1 4 1 25 1.5 4.7 41.0 85.000 325.00 310.0 1 3 1 26 -999.0 -999.0 16.2 36.330 119.50 63.0 1 1 1 27 3.4 10.4 9.0 0.101 4.00 28.0 5 1 3 28 0.8 7.4 7.6 1.040 5.50 68.0 5 3 4 29 0.8 2.1 46.0 521.000 655.00 336.0 5 5 5 30 -999.0 -999.0 22.4 100.000 157.00 100.0 1 1 1 31 -999.0 -999.0 16.3 35.000 56.00 33.0 3 5 4 32 1.4 7.7 2.6 0.005 0.14 21.5 5 2 4 33 2.0 17.9 24.0 0.010 0.25 50.0 1 1 1 34 1.9 6.1 100.0 62.000 1320.00 267.0 1 1 1 35 2.4 8.2 -999.0 0.122 3.00 30.0 2 1 1 36 2.8 8.4 -999.0 1.350 8.10 45.0 3 1 3 37 1.3 11.9 3.2 0.023 0.40 19.0 4 1 3 38 2.0 10.8 2.0 0.048 0.33 30.0 4 1 3 39 5.6 13.8 5.0 1.700 6.30 12.0 2 1 1 40 3.1 14.3 6.5 3.500 10.80 120.0 2 1 1 41 1.0 -999.0 23.6 250.000 490.00 440.0 5 5 5 42 1.8 15.2 12.0 0.480 15.50 140.0 2 2 2 43 0.9 10.0 20.2 10.000 115.00 170.0 4 4 4 44 1.8 11.9 13.0 1.620 11.40 17.0 2 1 2 45 1.9 6.5 27.0 192.000 180.00 115.0 4 4 4 46 0.9 7.5 18.0 2.500 12.10 31.0 5 5 5 47 -999.0 -999.0 13.7 4.288 39.20 63.0 2 2 2 48 2.6 10.6 4.7 0.280 1.90 21.0 3 1 3 49 2.4 7.4 9.8 4.235 50.40 52.0 1 1 1 50 1.2 8.4 29.0 6.800 179.00 164.0 2 3 2 51 0.9 5.7 7.0 0.750 12.30 225.0 2 2 2 52 0.5 4.9 6.0 3.600 21.00 225.0 3 2 3 53 -999.0 -999.0 17.0 14.830 98.20 150.0 5 5 5 54 0.6 3.2 20.0 55.500 175.00 151.0 5 5 5 55 -999.0 -999.0 12.7 1.400 12.50 90.0 2 2 2 56 2.2 8.1 3.5 0.060 1.00 -999.0 3 1 1 57 2.3 11.0 4.5 0.900 2.60 38.0 2 1 2 58 0.5 -999.0 7.5 2.000 17.00 200.0 3 1 3 59 -999.0 13.2 2.3 0.104 2.50 46.0 3 2 2 60 0.6 9.7 13.0 4.050 58.00 210.0 4 3 4 61 6.6 12.8 3.0 3.500 3.90 14.0 2 1 1 62 2.6 4.9 24.0 4.190 12.30 60.0 3 1 2 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) SWS L wb wbr tg -90.747120 0.747139 0.048076 -0.058508 0.050150 0.007202 P S `D\r` -56.503995 -57.168937 138.812334 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -912.45 -53.93 -14.48 55.57 704.36 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -90.747120 72.068777 -1.259 0.2135 SWS 0.747139 0.077462 9.645 2.93e-13 *** L 0.048076 0.129588 0.371 0.7121 wb -0.058508 0.099518 -0.588 0.5591 wbr 0.050150 0.099368 0.505 0.6159 tg 0.007202 0.118331 0.061 0.9517 P -56.503995 59.041082 -0.957 0.3429 S -57.168937 37.237168 -1.535 0.1307 `D\r` 138.812334 75.848946 1.830 0.0729 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 236.7 on 53 degrees of freedom Multiple R-squared: 0.6943, Adjusted R-squared: 0.6481 F-statistic: 15.05 on 8 and 53 DF, p-value: 3.107e-11 > 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,] 2.469785e-04 4.939569e-04 0.9997530 [2,] 1.993374e-05 3.986749e-05 0.9999801 [3,] 1.088382e-06 2.176765e-06 0.9999989 [4,] 5.024930e-08 1.004986e-07 0.9999999 [5,] 2.355591e-09 4.711183e-09 1.0000000 [6,] 9.158024e-11 1.831605e-10 1.0000000 [7,] 3.950450e-12 7.900900e-12 1.0000000 [8,] 1.476013e-13 2.952025e-13 1.0000000 [9,] 5.061939e-15 1.012388e-14 1.0000000 [10,] 3.630902e-01 7.261803e-01 0.6369098 [11,] 2.900966e-01 5.801932e-01 0.7099034 [12,] 2.215534e-01 4.431067e-01 0.7784466 [13,] 1.676603e-01 3.353205e-01 0.8323397 [14,] 1.340730e-01 2.681459e-01 0.8659270 [15,] 1.126999e-01 2.253997e-01 0.8873001 [16,] 7.632347e-02 1.526469e-01 0.9236765 [17,] 4.945799e-02 9.891599e-02 0.9505420 [18,] 5.150218e-02 1.030044e-01 0.9484978 [19,] 5.909410e-02 1.181882e-01 0.9409059 [20,] 5.741794e-02 1.148359e-01 0.9425821 [21,] 3.706056e-02 7.412112e-02 0.9629394 [22,] 2.301069e-02 4.602138e-02 0.9769893 [23,] 3.432208e-02 6.864417e-02 0.9656779 [24,] 3.102232e-02 6.204464e-02 0.9689777 [25,] 1.984750e-02 3.969501e-02 0.9801525 [26,] 1.207691e-02 2.415382e-02 0.9879231 [27,] 7.397858e-03 1.479572e-02 0.9926021 [28,] 5.218207e-03 1.043641e-02 0.9947818 [29,] 4.675678e-03 9.351355e-03 0.9953243 [30,] 5.710477e-02 1.142095e-01 0.9428952 [31,] 3.920766e-02 7.841531e-02 0.9607923 [32,] 2.357352e-02 4.714704e-02 0.9764265 [33,] 1.590174e-02 3.180348e-02 0.9840983 [34,] 8.767056e-03 1.753411e-02 0.9912329 [35,] 1.831791e-02 3.663583e-02 0.9816821 [36,] 1.799898e-02 3.599796e-02 0.9820010 [37,] 2.623105e-02 5.246210e-02 0.9737689 [38,] 1.553613e-02 3.107225e-02 0.9844639 [39,] 5.841097e-02 1.168219e-01 0.9415890 > postscript(file="/var/www/html/freestat/rcomp/tmp/1beia1292871063.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/html/freestat/rcomp/tmp/2beia1292871063.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/html/freestat/rcomp/tmp/3beia1292871063.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/html/freestat/rcomp/tmp/445zd1292871063.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/html/freestat/rcomp/tmp/545zd1292871063.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 6 -26.582934 -102.507232 -190.138963 -133.824473 -98.544272 -26.888435 7 8 9 10 11 12 56.523720 41.074754 147.742252 37.732174 38.535755 -34.113675 13 14 15 16 17 18 35.455051 -302.175305 62.474210 30.359270 38.916674 32.967247 19 20 21 22 23 24 125.255270 69.048466 704.364453 -43.091706 34.638991 -27.447170 25 26 27 28 29 30 162.404923 -192.100010 8.799763 -16.274330 -42.715407 -190.819934 31 32 33 34 35 36 -263.535334 -72.283902 52.708070 -6.350431 166.053339 -55.108704 37 38 39 40 41 42 -50.405311 -48.900014 116.857913 113.014079 698.192543 28.577406 43 44 45 46 47 48 -23.759449 -25.015915 -12.686495 -41.202591 -214.966921 -104.784731 49 50 51 52 53 54 59.353489 81.407049 34.579732 -47.752374 -293.512372 -44.318711 55 56 57 58 59 60 -213.943267 181.743555 -23.186892 645.346686 -912.451072 -78.435759 61 62 118.912475 36.784779 > postscript(file="/var/www/html/freestat/rcomp/tmp/645zd1292871063.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 -26.582934 NA 1 -102.507232 -26.582934 2 -190.138963 -102.507232 3 -133.824473 -190.138963 4 -98.544272 -133.824473 5 -26.888435 -98.544272 6 56.523720 -26.888435 7 41.074754 56.523720 8 147.742252 41.074754 9 37.732174 147.742252 10 38.535755 37.732174 11 -34.113675 38.535755 12 35.455051 -34.113675 13 -302.175305 35.455051 14 62.474210 -302.175305 15 30.359270 62.474210 16 38.916674 30.359270 17 32.967247 38.916674 18 125.255270 32.967247 19 69.048466 125.255270 20 704.364453 69.048466 21 -43.091706 704.364453 22 34.638991 -43.091706 23 -27.447170 34.638991 24 162.404923 -27.447170 25 -192.100010 162.404923 26 8.799763 -192.100010 27 -16.274330 8.799763 28 -42.715407 -16.274330 29 -190.819934 -42.715407 30 -263.535334 -190.819934 31 -72.283902 -263.535334 32 52.708070 -72.283902 33 -6.350431 52.708070 34 166.053339 -6.350431 35 -55.108704 166.053339 36 -50.405311 -55.108704 37 -48.900014 -50.405311 38 116.857913 -48.900014 39 113.014079 116.857913 40 698.192543 113.014079 41 28.577406 698.192543 42 -23.759449 28.577406 43 -25.015915 -23.759449 44 -12.686495 -25.015915 45 -41.202591 -12.686495 46 -214.966921 -41.202591 47 -104.784731 -214.966921 48 59.353489 -104.784731 49 81.407049 59.353489 50 34.579732 81.407049 51 -47.752374 34.579732 52 -293.512372 -47.752374 53 -44.318711 -293.512372 54 -213.943267 -44.318711 55 181.743555 -213.943267 56 -23.186892 181.743555 57 645.346686 -23.186892 58 -912.451072 645.346686 59 -78.435759 -912.451072 60 118.912475 -78.435759 61 36.784779 118.912475 62 NA 36.784779 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -102.507232 -26.582934 [2,] -190.138963 -102.507232 [3,] -133.824473 -190.138963 [4,] -98.544272 -133.824473 [5,] -26.888435 -98.544272 [6,] 56.523720 -26.888435 [7,] 41.074754 56.523720 [8,] 147.742252 41.074754 [9,] 37.732174 147.742252 [10,] 38.535755 37.732174 [11,] -34.113675 38.535755 [12,] 35.455051 -34.113675 [13,] -302.175305 35.455051 [14,] 62.474210 -302.175305 [15,] 30.359270 62.474210 [16,] 38.916674 30.359270 [17,] 32.967247 38.916674 [18,] 125.255270 32.967247 [19,] 69.048466 125.255270 [20,] 704.364453 69.048466 [21,] -43.091706 704.364453 [22,] 34.638991 -43.091706 [23,] -27.447170 34.638991 [24,] 162.404923 -27.447170 [25,] -192.100010 162.404923 [26,] 8.799763 -192.100010 [27,] -16.274330 8.799763 [28,] -42.715407 -16.274330 [29,] -190.819934 -42.715407 [30,] -263.535334 -190.819934 [31,] -72.283902 -263.535334 [32,] 52.708070 -72.283902 [33,] -6.350431 52.708070 [34,] 166.053339 -6.350431 [35,] -55.108704 166.053339 [36,] -50.405311 -55.108704 [37,] -48.900014 -50.405311 [38,] 116.857913 -48.900014 [39,] 113.014079 116.857913 [40,] 698.192543 113.014079 [41,] 28.577406 698.192543 [42,] -23.759449 28.577406 [43,] -25.015915 -23.759449 [44,] -12.686495 -25.015915 [45,] -41.202591 -12.686495 [46,] -214.966921 -41.202591 [47,] -104.784731 -214.966921 [48,] 59.353489 -104.784731 [49,] 81.407049 59.353489 [50,] 34.579732 81.407049 [51,] -47.752374 34.579732 [52,] -293.512372 -47.752374 [53,] -44.318711 -293.512372 [54,] -213.943267 -44.318711 [55,] 181.743555 -213.943267 [56,] -23.186892 181.743555 [57,] 645.346686 -23.186892 [58,] -912.451072 645.346686 [59,] -78.435759 -912.451072 [60,] 118.912475 -78.435759 [61,] 36.784779 118.912475 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -102.507232 -26.582934 2 -190.138963 -102.507232 3 -133.824473 -190.138963 4 -98.544272 -133.824473 5 -26.888435 -98.544272 6 56.523720 -26.888435 7 41.074754 56.523720 8 147.742252 41.074754 9 37.732174 147.742252 10 38.535755 37.732174 11 -34.113675 38.535755 12 35.455051 -34.113675 13 -302.175305 35.455051 14 62.474210 -302.175305 15 30.359270 62.474210 16 38.916674 30.359270 17 32.967247 38.916674 18 125.255270 32.967247 19 69.048466 125.255270 20 704.364453 69.048466 21 -43.091706 704.364453 22 34.638991 -43.091706 23 -27.447170 34.638991 24 162.404923 -27.447170 25 -192.100010 162.404923 26 8.799763 -192.100010 27 -16.274330 8.799763 28 -42.715407 -16.274330 29 -190.819934 -42.715407 30 -263.535334 -190.819934 31 -72.283902 -263.535334 32 52.708070 -72.283902 33 -6.350431 52.708070 34 166.053339 -6.350431 35 -55.108704 166.053339 36 -50.405311 -55.108704 37 -48.900014 -50.405311 38 116.857913 -48.900014 39 113.014079 116.857913 40 698.192543 113.014079 41 28.577406 698.192543 42 -23.759449 28.577406 43 -25.015915 -23.759449 44 -12.686495 -25.015915 45 -41.202591 -12.686495 46 -214.966921 -41.202591 47 -104.784731 -214.966921 48 59.353489 -104.784731 49 81.407049 59.353489 50 34.579732 81.407049 51 -47.752374 34.579732 52 -293.512372 -47.752374 53 -44.318711 -293.512372 54 -213.943267 -44.318711 55 181.743555 -213.943267 56 -23.186892 181.743555 57 645.346686 -23.186892 58 -912.451072 645.346686 59 -78.435759 -912.451072 60 118.912475 -78.435759 61 36.784779 118.912475 > 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/freestat/rcomp/tmp/7wfyy1292871063.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/html/freestat/rcomp/tmp/8wfyy1292871063.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/html/freestat/rcomp/tmp/976x11292871063.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') Warning messages: 1: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced 2: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced > par(opar) > dev.off() null device 1 > if (n > n25) { + postscript(file="/var/www/html/freestat/rcomp/tmp/1076x11292871063.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/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11s6e71292871063.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/freestat/rcomp/tmp/12epcd1292871063.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/freestat/rcomp/tmp/13lq961292871063.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/freestat/rcomp/tmp/14vh9r1292871063.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/freestat/rcomp/tmp/15zipx1292871063.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/freestat/rcomp/tmp/16vs561292871063.tab") + } > > try(system("convert tmp/1beia1292871063.ps tmp/1beia1292871063.png",intern=TRUE)) character(0) > try(system("convert tmp/2beia1292871063.ps tmp/2beia1292871063.png",intern=TRUE)) character(0) > try(system("convert tmp/3beia1292871063.ps tmp/3beia1292871063.png",intern=TRUE)) character(0) > try(system("convert tmp/445zd1292871063.ps tmp/445zd1292871063.png",intern=TRUE)) character(0) > try(system("convert tmp/545zd1292871063.ps tmp/545zd1292871063.png",intern=TRUE)) character(0) > try(system("convert tmp/645zd1292871063.ps tmp/645zd1292871063.png",intern=TRUE)) character(0) > try(system("convert tmp/7wfyy1292871063.ps tmp/7wfyy1292871063.png",intern=TRUE)) character(0) > try(system("convert tmp/8wfyy1292871063.ps tmp/8wfyy1292871063.png",intern=TRUE)) character(0) > try(system("convert tmp/976x11292871063.ps tmp/976x11292871063.png",intern=TRUE)) character(0) > try(system("convert tmp/1076x11292871063.ps tmp/1076x11292871063.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.895 2.420 4.490