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(6000 + ,10800 + ,10100 + ,16100 + ,17700 + ,13900 + ,17700 + ,6000 + ,10900 + ,10000 + ,15800 + ,17700 + ,13500 + ,19800 + ,6000 + ,11000 + ,10000 + ,16900 + ,17700 + ,13900 + ,19400 + ,6000 + ,11000 + ,10000 + ,17800 + ,17700 + ,13700 + ,18500 + ,6000 + ,11100 + ,10600 + ,17600 + ,17400 + ,13800 + ,18400 + ,6000 + ,11000 + ,12200 + ,18300 + ,17800 + ,15100 + ,18200 + ,6000 + ,11000 + ,12400 + ,18000 + ,17800 + ,15100 + ,18300 + ,6000 + ,11100 + ,13400 + ,15700 + ,17800 + ,14500 + ,19100 + ,6100 + ,11100 + ,13000 + ,14500 + ,17800 + ,13000 + ,18500 + ,6100 + ,11100 + ,10500 + ,14000 + ,18100 + ,12900 + ,18100 + ,6100 + ,11100 + ,10000 + ,15500 + ,18400 + ,14400 + ,18300 + ,6100 + ,11100 + ,10000 + ,15800 + ,18000 + ,14600 + ,17900 + ,6100 + ,11200 + ,10100 + ,15800 + ,17800 + ,15000 + ,18000 + ,6100 + ,11100 + ,10200 + ,15900 + ,17600 + ,13900 + ,18200 + ,6200 + ,11100 + ,10600 + ,18000 + ,17400 + ,14800 + ,18800 + ,6200 + ,11200 + ,10900 + ,19900 + ,17200 + ,15200 + ,20100 + ,6200 + ,11200 + ,10900 + ,20600 + ,17300 + ,16800 + ,19700 + ,6300 + ,11100 + ,11500 + ,20600 + ,17700 + ,17400 + ,19200 + ,6300 + ,11200 + ,12500 + ,20800 + ,18100 + ,17200 + ,19800 + ,6300 + ,11100 + ,13700 + ,20000 + ,18300 + ,17400 + ,20200 + ,6300 + ,11100 + ,15100 + ,18500 + ,18700 + ,18300 + ,19000 + ,6300 + ,11000 + ,13500 + ,17700 + ,18900 + ,19900 + ,19400 + ,6300 + ,11000 + ,13200 + ,17000 + ,18200 + ,18500 + ,19600 + ,6400 + ,11000 + ,13000 + ,16600 + ,17900 + ,16800 + ,18400 + ,6300 + ,11100 + ,13900 + ,16700 + ,17900 + ,16200 + ,18700 + ,6300 + ,11000 + ,14000 + ,17300 + ,18200 + ,16200 + ,18400 + ,6300 + ,11000 + ,13900 + ,19100 + ,18200 + ,16400 + ,20700 + ,6300 + ,10900 + ,14200 + ,20200 + ,18100 + ,15900 + ,20800 + ,6300 + ,11000 + ,14400 + ,20700 + ,18100 + ,16300 + ,21400 + ,6300 + ,11000 + ,14400 + ,21500 + ,17800 + ,16800 + ,21500 + ,6400 + ,11100 + ,14500 + ,21000 + ,18000 + ,15900 + ,20500 + ,6400 + ,11300 + ,13900 + ,16800 + ,17900 + ,15400 + ,20500 + ,6400 + ,11300 + ,14800 + ,16800 + ,18300 + ,15100 + ,19500 + ,6500 + ,11300 + ,13200 + ,16500 + ,18200 + ,15000 + ,20200 + ,6500 + ,11300 + ,12900 + ,17200 + ,18000 + ,17100 + ,20200 + ,6500 + ,11400 + ,13100 + ,17300 + ,18200 + ,16000 + ,18800 + ,6500 + ,11400 + ,12700 + ,17600 + ,18400 + ,15500 + ,19600 + ,6500 + ,11400 + ,13800 + ,18400 + ,18200 + ,16300 + ,19300 + ,6500 + ,11500 + ,13800 + ,19900 + ,18100 + ,16400 + ,20300 + ,6500 + ,11500 + ,14500 + ,20500 + ,17900 + ,16800 + ,21000 + ,6500 + ,11500 + ,15000 + ,21200 + ,18700 + ,17200 + ,19500 + ,6500 + ,11500 + ,16300 + ,21300 + ,18900 + ,17600 + ,20700 + ,6600 + ,11500 + ,17300 + ,20800 + ,19200 + ,18400 + ,20900 + ,6600 + ,11500 + ,18400 + ,18800 + ,19000 + ,18900 + ,20100 + ,6600 + ,11400 + ,17500 + ,18100 + ,19100 + ,18600 + ,19200 + ,6500 + ,11400 + ,13400 + ,18100 + ,19500 + ,18100 + ,19900 + ,6500 + ,11400 + ,13600 + ,18800 + ,20400 + ,18300 + ,21100 + ,6500 + ,11300 + ,13300 + ,18700 + ,19900 + ,17200 + ,20000 + ,6500 + ,11200 + ,13700 + ,18700 + ,19400 + ,15900 + ,20900 + ,6500 + ,11300 + ,13900 + ,19000 + ,19300 + ,16600 + ,20400 + ,6500 + ,11300 + ,14000 + ,20100 + ,18900 + ,15900 + ,20900 + ,6500 + ,11300 + ,14000 + ,20500 + ,18700 + ,16000 + ,20900 + ,6600 + ,11200 + ,14300 + ,21600 + ,18900 + ,15600 + ,21300 + ,6700 + ,11300 + ,15200 + ,21800 + ,19000 + ,16000 + ,21300 + ,6600 + ,11200 + ,15400 + ,21500 + ,19300 + ,16200 + ,21700 + ,6700 + ,11200 + ,18500 + ,21200 + ,19400 + ,16000 + ,21300 + ,6600 + ,11100 + ,18300 + ,20400 + ,18800 + ,16000 + ,20000 + ,6600 + ,11100 + ,12900 + ,20400 + ,18900 + ,16800 + ,20500 + ,6600 + ,11100 + ,12000 + ,20600 + ,19200 + ,17700 + ,20800 + ,6600 + ,11100 + ,12000 + ,19300 + ,19100 + ,17500 + ,20700 + ,7100 + ,11400 + ,12100 + ,18600 + ,18900 + ,17600 + ,21200 + ,7400 + ,11500 + ,12100 + ,19400 + ,18900 + ,18900 + ,21300 + ,7500 + ,11500 + ,11900 + ,23500 + ,19800 + ,18800 + ,21600 + ,7500 + ,11600 + ,11800 + ,24600 + ,20200 + ,19000 + ,22500 + ,7500 + ,11500 + ,11700 + ,25900 + ,20200 + ,19100 + ,22600 + ,7500 + ,11600 + ,12200 + ,26600 + ,19900 + ,19100 + ,23900 + ,7000 + ,11300 + ,12500 + ,24100 + ,19700 + ,18400 + ,23600 + ,6900 + ,11300 + ,13000 + ,21800 + ,19600 + ,16900 + ,22600 + ,6900 + ,11200 + ,13300 + ,21300 + ,19500 + ,16100 + ,22600 + ,6800 + ,11200 + ,11800 + ,21100 + ,19800 + ,16700 + ,22700 + ,6800 + ,11100 + ,11800 + ,21200 + ,20000 + ,18400 + ,22900 + ,6800 + ,11100 + ,11900 + ,21600 + ,20000 + ,18400 + ,22100) + ,dim=c(7 + ,72) + ,dimnames=list(c('Mineraalwater' + ,'Vruchtesappen' + ,'Appelen' + ,'Sinaasappelen' + ,'Citroenen' + ,'Pompelmoezen' + ,'Bananen') + ,1:72)) > y <- array(NA,dim=c(7,72),dimnames=list(c('Mineraalwater','Vruchtesappen','Appelen','Sinaasappelen','Citroenen','Pompelmoezen','Bananen'),1:72)) > 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 = '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 Vruchtesappen Mineraalwater Appelen Sinaasappelen Citroenen Pompelmoezen 1 10800 6000 10100 16100 17700 13900 2 10900 6000 10000 15800 17700 13500 3 11000 6000 10000 16900 17700 13900 4 11000 6000 10000 17800 17700 13700 5 11100 6000 10600 17600 17400 13800 6 11000 6000 12200 18300 17800 15100 7 11000 6000 12400 18000 17800 15100 8 11100 6000 13400 15700 17800 14500 9 11100 6100 13000 14500 17800 13000 10 11100 6100 10500 14000 18100 12900 11 11100 6100 10000 15500 18400 14400 12 11100 6100 10000 15800 18000 14600 13 11200 6100 10100 15800 17800 15000 14 11100 6100 10200 15900 17600 13900 15 11100 6200 10600 18000 17400 14800 16 11200 6200 10900 19900 17200 15200 17 11200 6200 10900 20600 17300 16800 18 11100 6300 11500 20600 17700 17400 19 11200 6300 12500 20800 18100 17200 20 11100 6300 13700 20000 18300 17400 21 11100 6300 15100 18500 18700 18300 22 11000 6300 13500 17700 18900 19900 23 11000 6300 13200 17000 18200 18500 24 11000 6400 13000 16600 17900 16800 25 11100 6300 13900 16700 17900 16200 26 11000 6300 14000 17300 18200 16200 27 11000 6300 13900 19100 18200 16400 28 10900 6300 14200 20200 18100 15900 29 11000 6300 14400 20700 18100 16300 30 11000 6300 14400 21500 17800 16800 31 11100 6400 14500 21000 18000 15900 32 11300 6400 13900 16800 17900 15400 33 11300 6400 14800 16800 18300 15100 34 11300 6500 13200 16500 18200 15000 35 11300 6500 12900 17200 18000 17100 36 11400 6500 13100 17300 18200 16000 37 11400 6500 12700 17600 18400 15500 38 11400 6500 13800 18400 18200 16300 39 11500 6500 13800 19900 18100 16400 40 11500 6500 14500 20500 17900 16800 41 11500 6500 15000 21200 18700 17200 42 11500 6500 16300 21300 18900 17600 43 11500 6600 17300 20800 19200 18400 44 11500 6600 18400 18800 19000 18900 45 11400 6600 17500 18100 19100 18600 46 11400 6500 13400 18100 19500 18100 47 11400 6500 13600 18800 20400 18300 48 11300 6500 13300 18700 19900 17200 49 11200 6500 13700 18700 19400 15900 50 11300 6500 13900 19000 19300 16600 51 11300 6500 14000 20100 18900 15900 52 11300 6500 14000 20500 18700 16000 53 11200 6600 14300 21600 18900 15600 54 11300 6700 15200 21800 19000 16000 55 11200 6600 15400 21500 19300 16200 56 11200 6700 18500 21200 19400 16000 57 11100 6600 18300 20400 18800 16000 58 11100 6600 12900 20400 18900 16800 59 11100 6600 12000 20600 19200 17700 60 11100 6600 12000 19300 19100 17500 61 11400 7100 12100 18600 18900 17600 62 11500 7400 12100 19400 18900 18900 63 11500 7500 11900 23500 19800 18800 64 11600 7500 11800 24600 20200 19000 65 11500 7500 11700 25900 20200 19100 66 11600 7500 12200 26600 19900 19100 67 11300 7000 12500 24100 19700 18400 68 11300 6900 13000 21800 19600 16900 69 11200 6900 13300 21300 19500 16100 70 11200 6800 11800 21100 19800 16700 71 11100 6800 11800 21200 20000 18400 72 11100 6800 11900 21600 20000 18400 Bananen M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t 1 17700 1 0 0 0 0 0 0 0 0 0 0 1 2 19800 0 1 0 0 0 0 0 0 0 0 0 2 3 19400 0 0 1 0 0 0 0 0 0 0 0 3 4 18500 0 0 0 1 0 0 0 0 0 0 0 4 5 18400 0 0 0 0 1 0 0 0 0 0 0 5 6 18200 0 0 0 0 0 1 0 0 0 0 0 6 7 18300 0 0 0 0 0 0 1 0 0 0 0 7 8 19100 0 0 0 0 0 0 0 1 0 0 0 8 9 18500 0 0 0 0 0 0 0 0 1 0 0 9 10 18100 0 0 0 0 0 0 0 0 0 1 0 10 11 18300 0 0 0 0 0 0 0 0 0 0 1 11 12 17900 0 0 0 0 0 0 0 0 0 0 0 12 13 18000 1 0 0 0 0 0 0 0 0 0 0 13 14 18200 0 1 0 0 0 0 0 0 0 0 0 14 15 18800 0 0 1 0 0 0 0 0 0 0 0 15 16 20100 0 0 0 1 0 0 0 0 0 0 0 16 17 19700 0 0 0 0 1 0 0 0 0 0 0 17 18 19200 0 0 0 0 0 1 0 0 0 0 0 18 19 19800 0 0 0 0 0 0 1 0 0 0 0 19 20 20200 0 0 0 0 0 0 0 1 0 0 0 20 21 19000 0 0 0 0 0 0 0 0 1 0 0 21 22 19400 0 0 0 0 0 0 0 0 0 1 0 22 23 19600 0 0 0 0 0 0 0 0 0 0 1 23 24 18400 0 0 0 0 0 0 0 0 0 0 0 24 25 18700 1 0 0 0 0 0 0 0 0 0 0 25 26 18400 0 1 0 0 0 0 0 0 0 0 0 26 27 20700 0 0 1 0 0 0 0 0 0 0 0 27 28 20800 0 0 0 1 0 0 0 0 0 0 0 28 29 21400 0 0 0 0 1 0 0 0 0 0 0 29 30 21500 0 0 0 0 0 1 0 0 0 0 0 30 31 20500 0 0 0 0 0 0 1 0 0 0 0 31 32 20500 0 0 0 0 0 0 0 1 0 0 0 32 33 19500 0 0 0 0 0 0 0 0 1 0 0 33 34 20200 0 0 0 0 0 0 0 0 0 1 0 34 35 20200 0 0 0 0 0 0 0 0 0 0 1 35 36 18800 0 0 0 0 0 0 0 0 0 0 0 36 37 19600 1 0 0 0 0 0 0 0 0 0 0 37 38 19300 0 1 0 0 0 0 0 0 0 0 0 38 39 20300 0 0 1 0 0 0 0 0 0 0 0 39 40 21000 0 0 0 1 0 0 0 0 0 0 0 40 41 19500 0 0 0 0 1 0 0 0 0 0 0 41 42 20700 0 0 0 0 0 1 0 0 0 0 0 42 43 20900 0 0 0 0 0 0 1 0 0 0 0 43 44 20100 0 0 0 0 0 0 0 1 0 0 0 44 45 19200 0 0 0 0 0 0 0 0 1 0 0 45 46 19900 0 0 0 0 0 0 0 0 0 1 0 46 47 21100 0 0 0 0 0 0 0 0 0 0 1 47 48 20000 0 0 0 0 0 0 0 0 0 0 0 48 49 20900 1 0 0 0 0 0 0 0 0 0 0 49 50 20400 0 1 0 0 0 0 0 0 0 0 0 50 51 20900 0 0 1 0 0 0 0 0 0 0 0 51 52 20900 0 0 0 1 0 0 0 0 0 0 0 52 53 21300 0 0 0 0 1 0 0 0 0 0 0 53 54 21300 0 0 0 0 0 1 0 0 0 0 0 54 55 21700 0 0 0 0 0 0 1 0 0 0 0 55 56 21300 0 0 0 0 0 0 0 1 0 0 0 56 57 20000 0 0 0 0 0 0 0 0 1 0 0 57 58 20500 0 0 0 0 0 0 0 0 0 1 0 58 59 20800 0 0 0 0 0 0 0 0 0 0 1 59 60 20700 0 0 0 0 0 0 0 0 0 0 0 60 61 21200 1 0 0 0 0 0 0 0 0 0 0 61 62 21300 0 1 0 0 0 0 0 0 0 0 0 62 63 21600 0 0 1 0 0 0 0 0 0 0 0 63 64 22500 0 0 0 1 0 0 0 0 0 0 0 64 65 22600 0 0 0 0 1 0 0 0 0 0 0 65 66 23900 0 0 0 0 0 1 0 0 0 0 0 66 67 23600 0 0 0 0 0 0 1 0 0 0 0 67 68 22600 0 0 0 0 0 0 0 1 0 0 0 68 69 22600 0 0 0 0 0 0 0 0 1 0 0 69 70 22700 0 0 0 0 0 0 0 0 0 1 0 70 71 22900 0 0 0 0 0 0 0 0 0 0 1 71 72 22100 0 0 0 0 0 0 0 0 0 0 0 72 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) Mineraalwater Appelen Sinaasappelen Citroenen 9.132e+03 4.541e-01 7.711e-03 -6.558e-02 3.405e-02 Pompelmoezen Bananen M1 M2 M3 1.926e-02 -4.908e-02 5.741e+01 7.775e+01 2.505e+02 M4 M5 M6 M7 M8 3.442e+02 3.560e+02 3.598e+02 3.029e+02 1.926e+02 M9 M10 M11 t 4.963e+01 3.993e+01 4.494e+01 4.556e+00 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -250.738 -95.685 -4.094 87.860 251.314 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 9.132e+03 9.919e+02 9.206 1.40e-12 *** Mineraalwater 4.541e-01 1.124e-01 4.041 0.000173 *** Appelen 7.711e-03 1.322e-02 0.583 0.562241 Sinaasappelen -6.558e-02 2.230e-02 -2.941 0.004845 ** Citroenen 3.405e-02 4.111e-02 0.828 0.411251 Pompelmoezen 1.926e-02 1.795e-02 1.073 0.288291 Bananen -4.908e-02 2.823e-02 -1.738 0.087974 . M1 5.741e+01 7.916e+01 0.725 0.471510 M2 7.775e+01 8.150e+01 0.954 0.344426 M3 2.505e+02 9.727e+01 2.575 0.012846 * M4 3.442e+02 1.105e+02 3.115 0.002964 ** M5 3.560e+02 1.157e+02 3.076 0.003313 ** M6 3.598e+02 1.230e+02 2.924 0.005071 ** M7 3.029e+02 1.142e+02 2.652 0.010524 * M8 1.926e+02 1.016e+02 1.896 0.063364 . M9 4.963e+01 9.273e+01 0.535 0.594735 M10 3.993e+01 7.775e+01 0.514 0.609705 M11 4.494e+01 7.995e+01 0.562 0.576481 t 4.556e+00 2.876e+00 1.584 0.119109 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 130 on 53 degrees of freedom Multiple R-squared: 0.6387, Adjusted R-squared: 0.516 F-statistic: 5.206 on 18 and 53 DF, p-value: 1.309e-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 + } [,1] [,2] [,3] [1,] 0.008290463 0.016580926 0.991709537 [2,] 0.007356980 0.014713960 0.992643020 [3,] 0.001887408 0.003774816 0.998112592 [4,] 0.002749855 0.005499711 0.997250145 [5,] 0.001260149 0.002520298 0.998739851 [6,] 0.023191675 0.046383350 0.976808325 [7,] 0.151391135 0.302782269 0.848608865 [8,] 0.150820988 0.301641976 0.849179012 [9,] 0.429895436 0.859790872 0.570104564 [10,] 0.652309434 0.695381131 0.347690566 [11,] 0.791665405 0.416669191 0.208334595 [12,] 0.795101446 0.409797108 0.204898554 [13,] 0.893025053 0.213949895 0.106974947 [14,] 0.960245263 0.079509474 0.039754737 [15,] 0.985135162 0.029729676 0.014864838 [16,] 0.979917958 0.040164083 0.020082042 [17,] 0.982945935 0.034108130 0.017054065 [18,] 0.988506008 0.022987984 0.011493992 [19,] 0.990855213 0.018289574 0.009144787 [20,] 0.998558070 0.002883859 0.001441930 [21,] 0.998355351 0.003289298 0.001644649 [22,] 0.998594832 0.002810336 0.001405168 [23,] 0.996576475 0.006847050 0.003423525 [24,] 0.993300133 0.013399734 0.006699867 [25,] 0.988161277 0.023677446 0.011838723 [26,] 0.971068839 0.057862322 0.028931161 [27,] 0.968801759 0.062396483 0.031198241 [28,] 0.974047417 0.051905167 0.025952583 [29,] 0.944797668 0.110404665 0.055202332 > postscript(file="/var/www/html/rcomp/tmp/1ir3r1292353954.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/rcomp/tmp/2tiku1292353954.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/rcomp/tmp/3tiku1292353954.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/rcomp/tmp/4tiku1292353954.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/rcomp/tmp/54r1f1292353954.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 = 72 Frequency = 1 1 2 3 4 5 6 -141.9391715 25.0220063 -7.4666255 -87.0630234 -17.7802358 -141.0062182 7 8 9 10 11 12 -104.9191983 -6.9471517 9.8693690 -26.4208772 36.9612464 87.1539572 13 14 15 16 17 18 128.4356664 47.1322567 -22.0070289 164.8920416 140.5962313 -67.4783734 19 20 21 22 23 24 110.0001717 62.9576162 2.3431388 -150.6323044 -143.1765556 -188.8364334 25 26 27 28 29 30 -79.4948537 -190.7545404 -140.2045893 -250.7382355 -114.1019852 -64.4661515 31 32 33 34 35 36 -29.5968180 18.3228547 92.8626249 84.9465740 89.9767269 181.0362202 37 38 39 40 41 42 183.9125653 179.6792062 251.3140504 220.4244362 137.5600267 170.1472780 43 44 45 46 47 48 120.8265671 44.8069712 2.4419301 114.9715012 174.1689680 94.5274716 49 50 51 52 53 54 15.7105505 74.3305586 20.0441767 -47.1410881 -118.5546458 -77.2201233 55 56 57 58 59 60 -95.0689557 -97.5326409 -108.0162446 -55.5045876 -57.8341437 -100.3643778 61 62 63 64 65 66 -106.6247571 -135.4094875 -101.6799834 -0.3741308 -27.7193912 180.0235884 67 68 69 70 71 72 -1.2417668 -21.6076496 0.4991818 32.6396940 -100.0962421 -73.5168377 > postscript(file="/var/www/html/rcomp/tmp/64r1f1292353954.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 = 72 Frequency = 1 lag(myerror, k = 1) myerror 0 -141.9391715 NA 1 25.0220063 -141.9391715 2 -7.4666255 25.0220063 3 -87.0630234 -7.4666255 4 -17.7802358 -87.0630234 5 -141.0062182 -17.7802358 6 -104.9191983 -141.0062182 7 -6.9471517 -104.9191983 8 9.8693690 -6.9471517 9 -26.4208772 9.8693690 10 36.9612464 -26.4208772 11 87.1539572 36.9612464 12 128.4356664 87.1539572 13 47.1322567 128.4356664 14 -22.0070289 47.1322567 15 164.8920416 -22.0070289 16 140.5962313 164.8920416 17 -67.4783734 140.5962313 18 110.0001717 -67.4783734 19 62.9576162 110.0001717 20 2.3431388 62.9576162 21 -150.6323044 2.3431388 22 -143.1765556 -150.6323044 23 -188.8364334 -143.1765556 24 -79.4948537 -188.8364334 25 -190.7545404 -79.4948537 26 -140.2045893 -190.7545404 27 -250.7382355 -140.2045893 28 -114.1019852 -250.7382355 29 -64.4661515 -114.1019852 30 -29.5968180 -64.4661515 31 18.3228547 -29.5968180 32 92.8626249 18.3228547 33 84.9465740 92.8626249 34 89.9767269 84.9465740 35 181.0362202 89.9767269 36 183.9125653 181.0362202 37 179.6792062 183.9125653 38 251.3140504 179.6792062 39 220.4244362 251.3140504 40 137.5600267 220.4244362 41 170.1472780 137.5600267 42 120.8265671 170.1472780 43 44.8069712 120.8265671 44 2.4419301 44.8069712 45 114.9715012 2.4419301 46 174.1689680 114.9715012 47 94.5274716 174.1689680 48 15.7105505 94.5274716 49 74.3305586 15.7105505 50 20.0441767 74.3305586 51 -47.1410881 20.0441767 52 -118.5546458 -47.1410881 53 -77.2201233 -118.5546458 54 -95.0689557 -77.2201233 55 -97.5326409 -95.0689557 56 -108.0162446 -97.5326409 57 -55.5045876 -108.0162446 58 -57.8341437 -55.5045876 59 -100.3643778 -57.8341437 60 -106.6247571 -100.3643778 61 -135.4094875 -106.6247571 62 -101.6799834 -135.4094875 63 -0.3741308 -101.6799834 64 -27.7193912 -0.3741308 65 180.0235884 -27.7193912 66 -1.2417668 180.0235884 67 -21.6076496 -1.2417668 68 0.4991818 -21.6076496 69 32.6396940 0.4991818 70 -100.0962421 32.6396940 71 -73.5168377 -100.0962421 72 NA -73.5168377 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] 25.0220063 -141.9391715 [2,] -7.4666255 25.0220063 [3,] -87.0630234 -7.4666255 [4,] -17.7802358 -87.0630234 [5,] -141.0062182 -17.7802358 [6,] -104.9191983 -141.0062182 [7,] -6.9471517 -104.9191983 [8,] 9.8693690 -6.9471517 [9,] -26.4208772 9.8693690 [10,] 36.9612464 -26.4208772 [11,] 87.1539572 36.9612464 [12,] 128.4356664 87.1539572 [13,] 47.1322567 128.4356664 [14,] -22.0070289 47.1322567 [15,] 164.8920416 -22.0070289 [16,] 140.5962313 164.8920416 [17,] -67.4783734 140.5962313 [18,] 110.0001717 -67.4783734 [19,] 62.9576162 110.0001717 [20,] 2.3431388 62.9576162 [21,] -150.6323044 2.3431388 [22,] -143.1765556 -150.6323044 [23,] -188.8364334 -143.1765556 [24,] -79.4948537 -188.8364334 [25,] -190.7545404 -79.4948537 [26,] -140.2045893 -190.7545404 [27,] -250.7382355 -140.2045893 [28,] -114.1019852 -250.7382355 [29,] -64.4661515 -114.1019852 [30,] -29.5968180 -64.4661515 [31,] 18.3228547 -29.5968180 [32,] 92.8626249 18.3228547 [33,] 84.9465740 92.8626249 [34,] 89.9767269 84.9465740 [35,] 181.0362202 89.9767269 [36,] 183.9125653 181.0362202 [37,] 179.6792062 183.9125653 [38,] 251.3140504 179.6792062 [39,] 220.4244362 251.3140504 [40,] 137.5600267 220.4244362 [41,] 170.1472780 137.5600267 [42,] 120.8265671 170.1472780 [43,] 44.8069712 120.8265671 [44,] 2.4419301 44.8069712 [45,] 114.9715012 2.4419301 [46,] 174.1689680 114.9715012 [47,] 94.5274716 174.1689680 [48,] 15.7105505 94.5274716 [49,] 74.3305586 15.7105505 [50,] 20.0441767 74.3305586 [51,] -47.1410881 20.0441767 [52,] -118.5546458 -47.1410881 [53,] -77.2201233 -118.5546458 [54,] -95.0689557 -77.2201233 [55,] -97.5326409 -95.0689557 [56,] -108.0162446 -97.5326409 [57,] -55.5045876 -108.0162446 [58,] -57.8341437 -55.5045876 [59,] -100.3643778 -57.8341437 [60,] -106.6247571 -100.3643778 [61,] -135.4094875 -106.6247571 [62,] -101.6799834 -135.4094875 [63,] -0.3741308 -101.6799834 [64,] -27.7193912 -0.3741308 [65,] 180.0235884 -27.7193912 [66,] -1.2417668 180.0235884 [67,] -21.6076496 -1.2417668 [68,] 0.4991818 -21.6076496 [69,] 32.6396940 0.4991818 [70,] -100.0962421 32.6396940 [71,] -73.5168377 -100.0962421 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 25.0220063 -141.9391715 2 -7.4666255 25.0220063 3 -87.0630234 -7.4666255 4 -17.7802358 -87.0630234 5 -141.0062182 -17.7802358 6 -104.9191983 -141.0062182 7 -6.9471517 -104.9191983 8 9.8693690 -6.9471517 9 -26.4208772 9.8693690 10 36.9612464 -26.4208772 11 87.1539572 36.9612464 12 128.4356664 87.1539572 13 47.1322567 128.4356664 14 -22.0070289 47.1322567 15 164.8920416 -22.0070289 16 140.5962313 164.8920416 17 -67.4783734 140.5962313 18 110.0001717 -67.4783734 19 62.9576162 110.0001717 20 2.3431388 62.9576162 21 -150.6323044 2.3431388 22 -143.1765556 -150.6323044 23 -188.8364334 -143.1765556 24 -79.4948537 -188.8364334 25 -190.7545404 -79.4948537 26 -140.2045893 -190.7545404 27 -250.7382355 -140.2045893 28 -114.1019852 -250.7382355 29 -64.4661515 -114.1019852 30 -29.5968180 -64.4661515 31 18.3228547 -29.5968180 32 92.8626249 18.3228547 33 84.9465740 92.8626249 34 89.9767269 84.9465740 35 181.0362202 89.9767269 36 183.9125653 181.0362202 37 179.6792062 183.9125653 38 251.3140504 179.6792062 39 220.4244362 251.3140504 40 137.5600267 220.4244362 41 170.1472780 137.5600267 42 120.8265671 170.1472780 43 44.8069712 120.8265671 44 2.4419301 44.8069712 45 114.9715012 2.4419301 46 174.1689680 114.9715012 47 94.5274716 174.1689680 48 15.7105505 94.5274716 49 74.3305586 15.7105505 50 20.0441767 74.3305586 51 -47.1410881 20.0441767 52 -118.5546458 -47.1410881 53 -77.2201233 -118.5546458 54 -95.0689557 -77.2201233 55 -97.5326409 -95.0689557 56 -108.0162446 -97.5326409 57 -55.5045876 -108.0162446 58 -57.8341437 -55.5045876 59 -100.3643778 -57.8341437 60 -106.6247571 -100.3643778 61 -135.4094875 -106.6247571 62 -101.6799834 -135.4094875 63 -0.3741308 -101.6799834 64 -27.7193912 -0.3741308 65 180.0235884 -27.7193912 66 -1.2417668 180.0235884 67 -21.6076496 -1.2417668 68 0.4991818 -21.6076496 69 32.6396940 0.4991818 70 -100.0962421 32.6396940 71 -73.5168377 -100.0962421 > 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/7wjii1292353954.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/rcomp/tmp/8wjii1292353954.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/rcomp/tmp/9psi31292353954.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/html/rcomp/tmp/10psi31292353954.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/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/11sayq1292353954.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/12ebxe1292353954.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/13s3vn1292353954.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/14dlbb1292353954.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/15z4rh1292353954.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/16k4qn1292353954.tab") + } > > try(system("convert tmp/1ir3r1292353954.ps tmp/1ir3r1292353954.png",intern=TRUE)) character(0) > try(system("convert tmp/2tiku1292353954.ps tmp/2tiku1292353954.png",intern=TRUE)) character(0) > try(system("convert tmp/3tiku1292353954.ps tmp/3tiku1292353954.png",intern=TRUE)) character(0) > try(system("convert tmp/4tiku1292353954.ps tmp/4tiku1292353954.png",intern=TRUE)) character(0) > try(system("convert tmp/54r1f1292353954.ps tmp/54r1f1292353954.png",intern=TRUE)) character(0) > try(system("convert tmp/64r1f1292353954.ps tmp/64r1f1292353954.png",intern=TRUE)) character(0) > try(system("convert tmp/7wjii1292353954.ps tmp/7wjii1292353954.png",intern=TRUE)) character(0) > try(system("convert tmp/8wjii1292353954.ps tmp/8wjii1292353954.png",intern=TRUE)) character(0) > try(system("convert tmp/9psi31292353954.ps tmp/9psi31292353954.png",intern=TRUE)) character(0) > try(system("convert tmp/10psi31292353954.ps tmp/10psi31292353954.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.675 1.703 16.526