R version 2.11.1 (2010-05-31) Copyright (C) 2010 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(1 + ,28 + ,6 + ,6 + ,6.06 + ,6.06 + ,3.53 + ,3.53 + ,48 + ,48 + ,5 + ,5 + ,1 + ,40 + ,5 + ,5 + ,8.1 + ,8.1 + ,4.52 + ,4.52 + ,63 + ,63 + ,11 + ,11 + ,1 + ,79 + ,3 + ,3 + ,79.38 + ,79.38 + ,3.72 + ,3.72 + ,113 + ,113 + ,13 + ,13 + ,1 + ,16 + ,2 + ,2 + ,26.26 + ,26.26 + ,3.17 + ,3.17 + ,104 + ,104 + ,1 + ,1 + ,1 + ,90 + ,2 + ,2 + ,39.56 + ,39.56 + ,3.39 + ,3.39 + ,89 + ,89 + ,11 + ,11 + ,1 + ,87 + ,5 + ,5 + ,65.61 + ,65.61 + ,4.15 + ,4.15 + ,97 + ,97 + ,3 + ,3 + ,1 + ,53 + ,5 + ,5 + ,80.3 + ,80.3 + ,3.09 + ,3.09 + ,114 + ,114 + ,11 + ,11 + ,1 + ,23 + ,5 + ,5 + ,34.68 + ,34.68 + ,2.76 + ,2.76 + ,57 + ,57 + ,9 + ,9 + ,1 + ,42 + ,6 + ,6 + ,7.17 + ,7.17 + ,5.14 + ,5.14 + ,127 + ,127 + ,10 + ,10 + ,1 + ,64 + ,4 + ,4 + ,65.88 + ,65.88 + ,4.78 + ,4.78 + ,64 + ,64 + ,4 + ,4 + ,1 + ,87 + ,6 + ,6 + ,42.69 + ,42.69 + ,4.22 + ,4.22 + ,91 + ,91 + ,2 + ,2 + ,1 + ,77 + ,2 + ,2 + ,54.94 + ,54.94 + ,3.93 + ,3.93 + ,127 + ,127 + ,2 + ,2 + ,1 + ,70 + ,4 + ,4 + ,89.99 + ,89.99 + ,3.01 + ,3.01 + ,45 + ,45 + ,10 + ,10 + ,1 + ,82 + ,4 + ,4 + ,72.64 + ,72.64 + ,5.12 + ,5.12 + ,40 + ,40 + ,9 + ,9 + ,1 + ,44 + ,3 + ,3 + ,24.96 + ,24.96 + ,5.82 + ,5.82 + ,33 + ,33 + ,1 + ,1 + ,1 + ,36 + ,2 + ,2 + ,57.52 + ,57.52 + ,2.83 + ,2.83 + ,60 + ,60 + ,7 + ,7 + ,0 + ,73 + ,2 + ,0 + ,71.91 + ,0 + ,5.11 + ,0 + ,50 + ,0 + ,3 + ,0 + ,0 + ,75 + ,3 + ,0 + ,65.34 + ,0 + ,5.99 + ,0 + ,128 + ,0 + ,11 + ,0 + ,0 + ,21 + ,3 + ,0 + ,34.62 + ,0 + ,3.15 + ,0 + ,52 + ,0 + ,7 + ,0 + ,0 + ,81 + ,2 + ,0 + ,60.92 + ,0 + ,3.5 + ,0 + ,40 + ,0 + ,1 + ,0 + ,0 + ,99 + ,3 + ,0 + ,56.49 + ,0 + ,4.5 + ,0 + ,29 + ,0 + ,9 + ,0 + ,0 + ,54 + ,3 + ,0 + ,56.19 + ,0 + ,3.31 + ,0 + ,36 + ,0 + ,5 + ,0 + ,0 + ,6 + ,4 + ,0 + ,61.2 + ,0 + ,5.31 + ,0 + ,49 + ,0 + ,9 + ,0 + ,0 + ,71 + ,5 + ,0 + ,58.2 + ,0 + ,4.24 + ,0 + ,57 + ,0 + ,7 + ,0 + ,0 + ,93 + ,6 + ,0 + ,75.91 + ,0 + ,5.06 + ,0 + ,82 + ,0 + ,4 + ,0 + ,0 + ,82 + ,3 + ,0 + ,73.66 + ,0 + ,4.72 + ,0 + ,34 + ,0 + ,10 + ,0 + ,0 + ,32 + ,4 + ,0 + ,73.87 + ,0 + ,4.58 + ,0 + ,36 + ,0 + ,13 + ,0 + ,0 + ,93 + ,4 + ,0 + ,87.21 + ,0 + ,5.3 + ,0 + ,89 + ,0 + ,9 + ,0 + ,0 + ,24 + ,4 + ,0 + ,64.29 + ,0 + ,5.11 + ,0 + ,69 + ,0 + ,5 + ,0 + ,0 + ,96 + ,5 + ,0 + ,71.82 + ,0 + ,4.05 + ,0 + ,35 + ,0 + ,8 + ,0 + ,0 + ,88 + ,4 + ,0 + ,89.31 + ,0 + ,4.62 + ,0 + ,65 + ,0 + ,12 + ,0 + ,0 + ,83 + ,2 + ,0 + ,1.41 + ,0 + ,4.66 + ,0 + ,70 + ,0 + ,8 + ,0 + ,0 + ,23 + ,6 + ,0 + ,35.17 + ,0 + ,4.66 + ,0 + ,60 + ,0 + ,5 + ,0 + ,0 + ,20 + ,5 + ,0 + ,41.08 + ,0 + ,5.1 + ,0 + ,127 + ,0 + ,11 + ,0 + ,0 + ,33 + ,3 + ,0 + ,30.57 + ,0 + ,4.97 + ,0 + ,96 + ,0 + ,8 + ,0 + ,0 + ,88 + ,2 + ,0 + ,68.84 + ,0 + ,2.87 + ,0 + ,61 + ,0 + ,9 + ,0 + ,0 + ,98 + ,2 + ,0 + ,71.05 + ,0 + ,4.98 + ,0 + ,36 + ,0 + ,1 + ,0 + ,0 + ,34 + ,4 + ,0 + ,23.32 + ,0 + ,4.55 + ,0 + ,55 + ,0 + ,9 + ,0 + ,0 + ,59 + ,3 + ,0 + ,61.39 + ,0 + ,5.45 + ,0 + ,75 + ,0 + ,2 + ,0 + ,0 + ,26 + ,6 + ,0 + ,8.41 + ,0 + ,4.36 + ,0 + ,42 + ,0 + ,3 + ,0 + ,0 + ,13 + ,1 + ,0 + ,64.06 + ,0 + ,4.74 + ,0 + ,83 + ,0 + ,3 + ,0 + ,0 + ,6 + ,2 + ,0 + ,26.8 + ,0 + ,5.44 + ,0 + ,56 + ,0 + ,1 + ,0 + ,0 + ,49 + ,4 + ,0 + ,12.78 + ,0 + ,5.78 + ,0 + ,114 + ,0 + ,5 + ,0 + ,0 + ,3 + ,5 + ,0 + ,23.84 + ,0 + ,2.92 + ,0 + ,33 + ,0 + ,4 + ,0 + ,0 + ,76 + ,4 + ,0 + ,5.68 + ,0 + ,3.22 + ,0 + ,80 + ,0 + ,6 + ,0 + ,0 + ,12 + ,2 + ,0 + ,45.92 + ,0 + ,3.04 + ,0 + ,115 + ,0 + ,7 + ,0 + ,0 + ,63 + ,5 + ,0 + ,18.17 + ,0 + ,3.11 + ,0 + ,127 + ,0 + ,13 + ,0 + ,0 + ,35 + ,1 + ,0 + ,29.12 + ,0 + ,3.87 + ,0 + ,45 + ,0 + ,9 + ,0 + ,0 + ,69 + ,1 + ,0 + ,40.08 + ,0 + ,3.75 + ,0 + ,74 + ,0 + ,11 + ,0 + ,0 + ,10 + ,5 + ,0 + ,1.08 + ,0 + ,4.82 + ,0 + ,105 + ,0 + ,10 + ,0) + ,dim=c(12 + ,50) + ,dimnames=list(c('klant' + ,'slaagkans' + ,'verzekeraar' + ,'verzekeraar_klant' + ,'kost' + ,'kost_klant' + ,'grootte' + ,'grootte_klant' + ,'snelheid' + ,'snelheid_klant' + ,'maand' + ,'maand_klant') + ,1:50)) > y <- array(NA,dim=c(12,50),dimnames=list(c('klant','slaagkans','verzekeraar','verzekeraar_klant','kost','kost_klant','grootte','grootte_klant','snelheid','snelheid_klant','maand','maand_klant'),1:50)) > 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 > 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 slaagkans klant verzekeraar verzekeraar_klant kost kost_klant grootte 1 28 1 6 6 6.06 6.06 3.53 2 40 1 5 5 8.10 8.10 4.52 3 79 1 3 3 79.38 79.38 3.72 4 16 1 2 2 26.26 26.26 3.17 5 90 1 2 2 39.56 39.56 3.39 6 87 1 5 5 65.61 65.61 4.15 7 53 1 5 5 80.30 80.30 3.09 8 23 1 5 5 34.68 34.68 2.76 9 42 1 6 6 7.17 7.17 5.14 10 64 1 4 4 65.88 65.88 4.78 11 87 1 6 6 42.69 42.69 4.22 12 77 1 2 2 54.94 54.94 3.93 13 70 1 4 4 89.99 89.99 3.01 14 82 1 4 4 72.64 72.64 5.12 15 44 1 3 3 24.96 24.96 5.82 16 36 1 2 2 57.52 57.52 2.83 17 73 0 2 0 71.91 0.00 5.11 18 75 0 3 0 65.34 0.00 5.99 19 21 0 3 0 34.62 0.00 3.15 20 81 0 2 0 60.92 0.00 3.50 21 99 0 3 0 56.49 0.00 4.50 22 54 0 3 0 56.19 0.00 3.31 23 6 0 4 0 61.20 0.00 5.31 24 71 0 5 0 58.20 0.00 4.24 25 93 0 6 0 75.91 0.00 5.06 26 82 0 3 0 73.66 0.00 4.72 27 32 0 4 0 73.87 0.00 4.58 28 93 0 4 0 87.21 0.00 5.30 29 24 0 4 0 64.29 0.00 5.11 30 96 0 5 0 71.82 0.00 4.05 31 88 0 4 0 89.31 0.00 4.62 32 83 0 2 0 1.41 0.00 4.66 33 23 0 6 0 35.17 0.00 4.66 34 20 0 5 0 41.08 0.00 5.10 35 33 0 3 0 30.57 0.00 4.97 36 88 0 2 0 68.84 0.00 2.87 37 98 0 2 0 71.05 0.00 4.98 38 34 0 4 0 23.32 0.00 4.55 39 59 0 3 0 61.39 0.00 5.45 40 26 0 6 0 8.41 0.00 4.36 41 13 0 1 0 64.06 0.00 4.74 42 6 0 2 0 26.80 0.00 5.44 43 49 0 4 0 12.78 0.00 5.78 44 3 0 5 0 23.84 0.00 2.92 45 76 0 4 0 5.68 0.00 3.22 46 12 0 2 0 45.92 0.00 3.04 47 63 0 5 0 18.17 0.00 3.11 48 35 0 1 0 29.12 0.00 3.87 49 69 0 1 0 40.08 0.00 3.75 50 10 0 5 0 1.08 0.00 4.82 grootte_klant snelheid snelheid_klant maand maand_klant 1 3.53 48 48 5 5 2 4.52 63 63 11 11 3 3.72 113 113 13 13 4 3.17 104 104 1 1 5 3.39 89 89 11 11 6 4.15 97 97 3 3 7 3.09 114 114 11 11 8 2.76 57 57 9 9 9 5.14 127 127 10 10 10 4.78 64 64 4 4 11 4.22 91 91 2 2 12 3.93 127 127 2 2 13 3.01 45 45 10 10 14 5.12 40 40 9 9 15 5.82 33 33 1 1 16 2.83 60 60 7 7 17 0.00 50 0 3 0 18 0.00 128 0 11 0 19 0.00 52 0 7 0 20 0.00 40 0 1 0 21 0.00 29 0 9 0 22 0.00 36 0 5 0 23 0.00 49 0 9 0 24 0.00 57 0 7 0 25 0.00 82 0 4 0 26 0.00 34 0 10 0 27 0.00 36 0 13 0 28 0.00 89 0 9 0 29 0.00 69 0 5 0 30 0.00 35 0 8 0 31 0.00 65 0 12 0 32 0.00 70 0 8 0 33 0.00 60 0 5 0 34 0.00 127 0 11 0 35 0.00 96 0 8 0 36 0.00 61 0 9 0 37 0.00 36 0 1 0 38 0.00 55 0 9 0 39 0.00 75 0 2 0 40 0.00 42 0 3 0 41 0.00 83 0 3 0 42 0.00 56 0 1 0 43 0.00 114 0 5 0 44 0.00 33 0 4 0 45 0.00 80 0 6 0 46 0.00 115 0 7 0 47 0.00 127 0 13 0 48 0.00 45 0 9 0 49 0.00 74 0 11 0 50 0.00 105 0 10 0 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) klant verzekeraar verzekeraar_klant 36.78638 -70.40610 -1.04099 1.39545 kost kost_klant grootte grootte_klant 0.55479 0.05675 -2.24036 13.59903 snelheid snelheid_klant maand maand_klant -0.04266 0.22118 0.85427 -0.61258 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -54.277 -15.175 1.977 15.978 54.106 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 36.78638 30.30401 1.214 0.2323 klant -70.40610 59.29137 -1.187 0.2424 verzekeraar -1.04099 3.52878 -0.295 0.7696 verzekeraar_klant 1.39545 6.37253 0.219 0.8278 kost 0.55479 0.21581 2.571 0.0142 * kost_klant 0.05675 0.36290 0.156 0.8766 grootte -2.24036 6.31464 -0.355 0.7247 grootte_klant 13.59903 10.75192 1.265 0.2136 snelheid -0.04266 0.19360 -0.220 0.8268 snelheid_klant 0.22118 0.30496 0.725 0.4727 maand 0.85427 1.55729 0.549 0.5865 maand_klant -0.61258 2.47010 -0.248 0.8055 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 28.8 on 38 degrees of freedom Multiple R-squared: 0.2858, Adjusted R-squared: 0.07908 F-statistic: 1.383 on 11 and 38 DF, p-value: 0.2208 > 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.67388239 0.65223522 0.3261176 [2,] 0.52253176 0.95493649 0.4774682 [3,] 0.36672276 0.73344552 0.6332772 [4,] 0.24138290 0.48276579 0.7586171 [5,] 0.15692286 0.31384573 0.8430771 [6,] 0.09239675 0.18479350 0.9076033 [7,] 0.05889686 0.11779372 0.9411031 [8,] 0.02992166 0.05984332 0.9700783 [9,] 0.02847258 0.05694515 0.9715274 [10,] 0.09671493 0.19342987 0.9032851 [11,] 0.08325173 0.16650346 0.9167483 [12,] 0.05862842 0.11725685 0.9413716 [13,] 0.09677436 0.19354872 0.9032256 [14,] 0.06433422 0.12866844 0.9356658 [15,] 0.06728468 0.13456935 0.9327153 [16,] 0.05195056 0.10390113 0.9480494 [17,] 0.03023170 0.06046340 0.9697683 [18,] 0.06943095 0.13886189 0.9305691 [19,] 0.04222474 0.08444948 0.9577753 [20,] 0.03432928 0.06865855 0.9656707 [21,] 0.01537506 0.03075012 0.9846249 > postscript(file="/var/www/rcomp/tmp/117yk1290537364.ps",horizontal=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/2ugx51290537364.ps",horizontal=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/3ugx51290537364.ps",horizontal=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/4ugx51290537364.ps",horizontal=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/5ugx51290537364.ps",horizontal=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 = 50 Frequency = 1 1 2 3 4 5 6 5.9135175 1.6474533 -2.5569342 -21.9630338 41.6654450 13.5440374 7 8 9 10 11 12 -22.3676766 -10.0616728 -14.3641874 -10.7731294 27.7238761 8.5176422 13 14 15 16 17 18 2.5291144 2.3069113 -10.9481076 -10.8132555 9.4189773 14.5699550 19 20 21 22 23 24 -28.5746612 21.1910856 37.6267543 -6.1571387 -54.2774071 13.0805990 25 26 27 28 29 30 31.7626428 10.9528846 -40.9137611 19.9765000 -34.1694889 28.3058289 31 32 33 34 35 36 8.7013103 54.1056078 -18.3240694 -26.9253825 -9.2274455 16.4474678 37 38 39 40 41 42 35.7161202 -6.7085608 4.9789152 -0.2093026 -46.6879812 -29.8505007 43 44 45 46 47 48 22.8287048 -37.2750804 45.7276473 -42.4435547 22.6180714 -13.9993792 49 50 13.1799076 -15.4452658 > postscript(file="/var/www/rcomp/tmp/6npw81290537364.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > dum <- cbind(lag(myerror,k=1),myerror) > dum Time Series: Start = 0 End = 50 Frequency = 1 lag(myerror, k = 1) myerror 0 5.9135175 NA 1 1.6474533 5.9135175 2 -2.5569342 1.6474533 3 -21.9630338 -2.5569342 4 41.6654450 -21.9630338 5 13.5440374 41.6654450 6 -22.3676766 13.5440374 7 -10.0616728 -22.3676766 8 -14.3641874 -10.0616728 9 -10.7731294 -14.3641874 10 27.7238761 -10.7731294 11 8.5176422 27.7238761 12 2.5291144 8.5176422 13 2.3069113 2.5291144 14 -10.9481076 2.3069113 15 -10.8132555 -10.9481076 16 9.4189773 -10.8132555 17 14.5699550 9.4189773 18 -28.5746612 14.5699550 19 21.1910856 -28.5746612 20 37.6267543 21.1910856 21 -6.1571387 37.6267543 22 -54.2774071 -6.1571387 23 13.0805990 -54.2774071 24 31.7626428 13.0805990 25 10.9528846 31.7626428 26 -40.9137611 10.9528846 27 19.9765000 -40.9137611 28 -34.1694889 19.9765000 29 28.3058289 -34.1694889 30 8.7013103 28.3058289 31 54.1056078 8.7013103 32 -18.3240694 54.1056078 33 -26.9253825 -18.3240694 34 -9.2274455 -26.9253825 35 16.4474678 -9.2274455 36 35.7161202 16.4474678 37 -6.7085608 35.7161202 38 4.9789152 -6.7085608 39 -0.2093026 4.9789152 40 -46.6879812 -0.2093026 41 -29.8505007 -46.6879812 42 22.8287048 -29.8505007 43 -37.2750804 22.8287048 44 45.7276473 -37.2750804 45 -42.4435547 45.7276473 46 22.6180714 -42.4435547 47 -13.9993792 22.6180714 48 13.1799076 -13.9993792 49 -15.4452658 13.1799076 50 NA -15.4452658 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] 1.6474533 5.9135175 [2,] -2.5569342 1.6474533 [3,] -21.9630338 -2.5569342 [4,] 41.6654450 -21.9630338 [5,] 13.5440374 41.6654450 [6,] -22.3676766 13.5440374 [7,] -10.0616728 -22.3676766 [8,] -14.3641874 -10.0616728 [9,] -10.7731294 -14.3641874 [10,] 27.7238761 -10.7731294 [11,] 8.5176422 27.7238761 [12,] 2.5291144 8.5176422 [13,] 2.3069113 2.5291144 [14,] -10.9481076 2.3069113 [15,] -10.8132555 -10.9481076 [16,] 9.4189773 -10.8132555 [17,] 14.5699550 9.4189773 [18,] -28.5746612 14.5699550 [19,] 21.1910856 -28.5746612 [20,] 37.6267543 21.1910856 [21,] -6.1571387 37.6267543 [22,] -54.2774071 -6.1571387 [23,] 13.0805990 -54.2774071 [24,] 31.7626428 13.0805990 [25,] 10.9528846 31.7626428 [26,] -40.9137611 10.9528846 [27,] 19.9765000 -40.9137611 [28,] -34.1694889 19.9765000 [29,] 28.3058289 -34.1694889 [30,] 8.7013103 28.3058289 [31,] 54.1056078 8.7013103 [32,] -18.3240694 54.1056078 [33,] -26.9253825 -18.3240694 [34,] -9.2274455 -26.9253825 [35,] 16.4474678 -9.2274455 [36,] 35.7161202 16.4474678 [37,] -6.7085608 35.7161202 [38,] 4.9789152 -6.7085608 [39,] -0.2093026 4.9789152 [40,] -46.6879812 -0.2093026 [41,] -29.8505007 -46.6879812 [42,] 22.8287048 -29.8505007 [43,] -37.2750804 22.8287048 [44,] 45.7276473 -37.2750804 [45,] -42.4435547 45.7276473 [46,] 22.6180714 -42.4435547 [47,] -13.9993792 22.6180714 [48,] 13.1799076 -13.9993792 [49,] -15.4452658 13.1799076 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 1.6474533 5.9135175 2 -2.5569342 1.6474533 3 -21.9630338 -2.5569342 4 41.6654450 -21.9630338 5 13.5440374 41.6654450 6 -22.3676766 13.5440374 7 -10.0616728 -22.3676766 8 -14.3641874 -10.0616728 9 -10.7731294 -14.3641874 10 27.7238761 -10.7731294 11 8.5176422 27.7238761 12 2.5291144 8.5176422 13 2.3069113 2.5291144 14 -10.9481076 2.3069113 15 -10.8132555 -10.9481076 16 9.4189773 -10.8132555 17 14.5699550 9.4189773 18 -28.5746612 14.5699550 19 21.1910856 -28.5746612 20 37.6267543 21.1910856 21 -6.1571387 37.6267543 22 -54.2774071 -6.1571387 23 13.0805990 -54.2774071 24 31.7626428 13.0805990 25 10.9528846 31.7626428 26 -40.9137611 10.9528846 27 19.9765000 -40.9137611 28 -34.1694889 19.9765000 29 28.3058289 -34.1694889 30 8.7013103 28.3058289 31 54.1056078 8.7013103 32 -18.3240694 54.1056078 33 -26.9253825 -18.3240694 34 -9.2274455 -26.9253825 35 16.4474678 -9.2274455 36 35.7161202 16.4474678 37 -6.7085608 35.7161202 38 4.9789152 -6.7085608 39 -0.2093026 4.9789152 40 -46.6879812 -0.2093026 41 -29.8505007 -46.6879812 42 22.8287048 -29.8505007 43 -37.2750804 22.8287048 44 45.7276473 -37.2750804 45 -42.4435547 45.7276473 46 22.6180714 -42.4435547 47 -13.9993792 22.6180714 48 13.1799076 -13.9993792 49 -15.4452658 13.1799076 > 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/7gzdb1290537364.ps",horizontal=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/8gzdb1290537364.ps",horizontal=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/9gzdb1290537364.ps",horizontal=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/1088ve1290537364.ps",horizontal=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/11uqt21290537364.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/12f9sq1290537364.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/13mspj1290537364.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/14fj651290537364.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/150kms1290537364.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/16wbkj1290537364.tab") + } > > try(system("convert tmp/117yk1290537364.ps tmp/117yk1290537364.png",intern=TRUE)) character(0) > try(system("convert tmp/2ugx51290537364.ps tmp/2ugx51290537364.png",intern=TRUE)) character(0) > try(system("convert tmp/3ugx51290537364.ps tmp/3ugx51290537364.png",intern=TRUE)) character(0) > try(system("convert tmp/4ugx51290537364.ps tmp/4ugx51290537364.png",intern=TRUE)) character(0) > try(system("convert tmp/5ugx51290537364.ps tmp/5ugx51290537364.png",intern=TRUE)) character(0) > try(system("convert tmp/6npw81290537364.ps tmp/6npw81290537364.png",intern=TRUE)) character(0) > try(system("convert tmp/7gzdb1290537364.ps tmp/7gzdb1290537364.png",intern=TRUE)) character(0) > try(system("convert tmp/8gzdb1290537364.ps tmp/8gzdb1290537364.png",intern=TRUE)) character(0) > try(system("convert tmp/9gzdb1290537364.ps tmp/9gzdb1290537364.png",intern=TRUE)) character(0) > try(system("convert tmp/1088ve1290537364.ps tmp/1088ve1290537364.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.580 2.090 5.634