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(107.1 + ,0 + ,96.3 + ,87.0 + ,96.8 + ,115.2 + ,0 + ,107.1 + ,96.3 + ,87.0 + ,106.1 + ,0 + ,115.2 + ,107.1 + ,96.3 + ,89.5 + ,0 + ,106.1 + ,115.2 + ,107.1 + ,91.3 + ,0 + ,89.5 + ,106.1 + ,115.2 + ,97.6 + ,0 + ,91.3 + ,89.5 + ,106.1 + ,100.7 + ,0 + ,97.6 + ,91.3 + ,89.5 + ,104.6 + ,0 + ,100.7 + ,97.6 + ,91.3 + ,94.7 + ,0 + ,104.6 + ,100.7 + ,97.6 + ,101.8 + ,0 + ,94.7 + ,104.6 + ,100.7 + ,102.5 + ,0 + ,101.8 + ,94.7 + ,104.6 + ,105.3 + ,0 + ,102.5 + ,101.8 + ,94.7 + ,110.3 + ,0 + ,105.3 + ,102.5 + ,101.8 + ,109.8 + ,0 + ,110.3 + ,105.3 + ,102.5 + ,117.3 + ,0 + ,109.8 + ,110.3 + ,105.3 + ,118.8 + ,0 + ,117.3 + ,109.8 + ,110.3 + ,131.3 + ,0 + ,118.8 + ,117.3 + ,109.8 + ,125.9 + ,0 + ,131.3 + ,118.8 + ,117.3 + ,133.1 + ,0 + ,125.9 + ,131.3 + ,118.8 + ,147.0 + ,0 + ,133.1 + ,125.9 + ,131.3 + ,145.8 + ,0 + ,147.0 + ,133.1 + ,125.9 + ,164.4 + ,0 + ,145.8 + ,147.0 + ,133.1 + ,149.8 + ,0 + ,164.4 + ,145.8 + ,147.0 + ,137.7 + ,0 + ,149.8 + ,164.4 + ,145.8 + ,151.7 + ,0 + ,137.7 + ,149.8 + ,164.4 + ,156.8 + ,0 + ,151.7 + ,137.7 + ,149.8 + ,180.0 + ,0 + ,156.8 + ,151.7 + ,137.7 + ,180.4 + ,0 + ,180.0 + ,156.8 + ,151.7 + ,170.4 + ,0 + ,180.4 + ,180.0 + ,156.8 + ,191.6 + ,0 + ,170.4 + ,180.4 + ,180.0 + ,199.5 + ,0 + ,191.6 + ,170.4 + ,180.4 + ,218.2 + ,0 + ,199.5 + ,191.6 + ,170.4 + ,217.5 + ,0 + ,218.2 + ,199.5 + ,191.6 + ,205.0 + ,0 + ,217.5 + ,218.2 + ,199.5 + ,194.0 + ,0 + ,205.0 + ,217.5 + ,218.2 + ,199.3 + ,0 + ,194.0 + ,205.0 + ,217.5 + ,219.3 + ,0 + ,199.3 + ,194.0 + ,205.0 + ,211.1 + ,0 + ,219.3 + ,199.3 + ,194.0 + ,215.2 + ,0 + ,211.1 + ,219.3 + ,199.3 + ,240.2 + ,0 + ,215.2 + ,211.1 + ,219.3 + ,242.2 + ,0 + ,240.2 + ,215.2 + ,211.1 + ,240.7 + ,0 + ,242.2 + ,240.2 + ,215.2 + ,255.4 + ,0 + ,240.7 + ,242.2 + ,240.2 + ,253.0 + ,0 + ,255.4 + ,240.7 + ,242.2 + ,218.2 + ,0 + ,253.0 + ,255.4 + ,240.7 + ,203.7 + ,0 + ,218.2 + ,253.0 + ,255.4 + ,205.6 + ,0 + ,203.7 + ,218.2 + ,253.0 + ,215.6 + ,0 + ,205.6 + ,203.7 + ,218.2 + ,188.5 + ,0 + ,215.6 + ,205.6 + ,203.7 + ,202.9 + ,0 + ,188.5 + ,215.6 + ,205.6 + ,214.0 + ,0 + ,202.9 + ,188.5 + ,215.6 + ,230.3 + ,0 + ,214.0 + ,202.9 + ,188.5 + ,230.0 + ,0 + ,230.3 + ,214.0 + ,202.9 + ,241.0 + ,0 + ,230.0 + ,230.3 + ,214.0 + ,259.6 + ,1 + ,241.0 + ,230.0 + ,230.3 + ,247.8 + ,1 + ,259.6 + ,241.0 + ,230.0 + ,270.3 + ,1 + ,247.8 + ,259.6 + ,241.0 + ,289.7 + ,1 + ,270.3 + ,247.8 + ,259.6 + ,322.7 + ,1 + ,289.7 + ,270.3 + ,247.8 + ,315.0 + ,1 + ,322.7 + ,289.7 + ,270.3 + ,320.2 + ,1 + ,315.0 + ,322.7 + ,289.7 + ,329.5 + ,1 + ,320.2 + ,315.0 + ,322.7 + ,360.6 + ,1 + ,329.5 + ,320.2 + ,315.0 + ,382.2 + ,1 + ,360.6 + ,329.5 + ,320.2 + ,435.4 + ,1 + ,382.2 + ,360.6 + ,329.5 + ,464.0 + ,1 + ,435.4 + ,382.2 + ,360.6 + ,468.8 + ,1 + ,464.0 + ,435.4 + ,382.2 + ,403.0 + ,1 + ,468.8 + ,464.0 + ,435.4 + ,351.6 + ,1 + ,403.0 + ,468.8 + ,464.0 + ,252.0 + ,1 + ,351.6 + ,403.0 + ,468.8 + ,188.0 + ,1 + ,252.0 + ,351.6 + ,403.0 + ,146.5 + ,1 + ,188.0 + ,252.0 + ,351.6 + ,152.9 + ,1 + ,146.5 + ,188.0 + ,252.0 + ,148.1 + ,1 + ,152.9 + ,146.5 + ,188.0 + ,165.1 + ,1 + ,148.1 + ,152.9 + ,146.5 + ,177.0 + ,1 + ,165.1 + ,148.1 + ,152.9 + ,206.1 + ,1 + ,177.0 + ,165.1 + ,148.1 + ,244.9 + ,1 + ,206.1 + ,177.0 + ,165.1 + ,228.6 + ,1 + ,244.9 + ,206.1 + ,177.0 + ,253.4 + ,1 + ,228.6 + ,244.9 + ,206.1 + ,241.1 + ,1 + ,253.4 + ,228.6 + ,244.9) + ,dim=c(5 + ,81) + ,dimnames=list(c('Y' + ,'D' + ,'Y1' + ,'Y2' + ,'Y3') + ,1:81)) > y <- array(NA,dim=c(5,81),dimnames=list(c('Y','D','Y1','Y2','Y3'),1:81)) > for (i in 1:dim(x)[1]) + { + for (j in 1:dim(x)[2]) + { + y[i,j] <- as.numeric(x[i,j]) + } + } > par3 = 'Linear Trend' > par2 = 'Include Monthly Dummies' > par1 = '1' > #'GNU S' R Code compiled by R2WASP v. 1.0.44 () > #Author: Prof. Dr. P. Wessa > #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/ > #Source of accompanying publication: Office for Research, Development, and Education > #Technical description: Write here your technical program description (don't use hard returns!) > library(lattice) > library(lmtest) Loading required package: zoo 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 Y D Y1 Y2 Y3 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t 1 107.1 0 96.3 87.0 96.8 1 0 0 0 0 0 0 0 0 0 0 1 2 115.2 0 107.1 96.3 87.0 0 1 0 0 0 0 0 0 0 0 0 2 3 106.1 0 115.2 107.1 96.3 0 0 1 0 0 0 0 0 0 0 0 3 4 89.5 0 106.1 115.2 107.1 0 0 0 1 0 0 0 0 0 0 0 4 5 91.3 0 89.5 106.1 115.2 0 0 0 0 1 0 0 0 0 0 0 5 6 97.6 0 91.3 89.5 106.1 0 0 0 0 0 1 0 0 0 0 0 6 7 100.7 0 97.6 91.3 89.5 0 0 0 0 0 0 1 0 0 0 0 7 8 104.6 0 100.7 97.6 91.3 0 0 0 0 0 0 0 1 0 0 0 8 9 94.7 0 104.6 100.7 97.6 0 0 0 0 0 0 0 0 1 0 0 9 10 101.8 0 94.7 104.6 100.7 0 0 0 0 0 0 0 0 0 1 0 10 11 102.5 0 101.8 94.7 104.6 0 0 0 0 0 0 0 0 0 0 1 11 12 105.3 0 102.5 101.8 94.7 0 0 0 0 0 0 0 0 0 0 0 12 13 110.3 0 105.3 102.5 101.8 1 0 0 0 0 0 0 0 0 0 0 13 14 109.8 0 110.3 105.3 102.5 0 1 0 0 0 0 0 0 0 0 0 14 15 117.3 0 109.8 110.3 105.3 0 0 1 0 0 0 0 0 0 0 0 15 16 118.8 0 117.3 109.8 110.3 0 0 0 1 0 0 0 0 0 0 0 16 17 131.3 0 118.8 117.3 109.8 0 0 0 0 1 0 0 0 0 0 0 17 18 125.9 0 131.3 118.8 117.3 0 0 0 0 0 1 0 0 0 0 0 18 19 133.1 0 125.9 131.3 118.8 0 0 0 0 0 0 1 0 0 0 0 19 20 147.0 0 133.1 125.9 131.3 0 0 0 0 0 0 0 1 0 0 0 20 21 145.8 0 147.0 133.1 125.9 0 0 0 0 0 0 0 0 1 0 0 21 22 164.4 0 145.8 147.0 133.1 0 0 0 0 0 0 0 0 0 1 0 22 23 149.8 0 164.4 145.8 147.0 0 0 0 0 0 0 0 0 0 0 1 23 24 137.7 0 149.8 164.4 145.8 0 0 0 0 0 0 0 0 0 0 0 24 25 151.7 0 137.7 149.8 164.4 1 0 0 0 0 0 0 0 0 0 0 25 26 156.8 0 151.7 137.7 149.8 0 1 0 0 0 0 0 0 0 0 0 26 27 180.0 0 156.8 151.7 137.7 0 0 1 0 0 0 0 0 0 0 0 27 28 180.4 0 180.0 156.8 151.7 0 0 0 1 0 0 0 0 0 0 0 28 29 170.4 0 180.4 180.0 156.8 0 0 0 0 1 0 0 0 0 0 0 29 30 191.6 0 170.4 180.4 180.0 0 0 0 0 0 1 0 0 0 0 0 30 31 199.5 0 191.6 170.4 180.4 0 0 0 0 0 0 1 0 0 0 0 31 32 218.2 0 199.5 191.6 170.4 0 0 0 0 0 0 0 1 0 0 0 32 33 217.5 0 218.2 199.5 191.6 0 0 0 0 0 0 0 0 1 0 0 33 34 205.0 0 217.5 218.2 199.5 0 0 0 0 0 0 0 0 0 1 0 34 35 194.0 0 205.0 217.5 218.2 0 0 0 0 0 0 0 0 0 0 1 35 36 199.3 0 194.0 205.0 217.5 0 0 0 0 0 0 0 0 0 0 0 36 37 219.3 0 199.3 194.0 205.0 1 0 0 0 0 0 0 0 0 0 0 37 38 211.1 0 219.3 199.3 194.0 0 1 0 0 0 0 0 0 0 0 0 38 39 215.2 0 211.1 219.3 199.3 0 0 1 0 0 0 0 0 0 0 0 39 40 240.2 0 215.2 211.1 219.3 0 0 0 1 0 0 0 0 0 0 0 40 41 242.2 0 240.2 215.2 211.1 0 0 0 0 1 0 0 0 0 0 0 41 42 240.7 0 242.2 240.2 215.2 0 0 0 0 0 1 0 0 0 0 0 42 43 255.4 0 240.7 242.2 240.2 0 0 0 0 0 0 1 0 0 0 0 43 44 253.0 0 255.4 240.7 242.2 0 0 0 0 0 0 0 1 0 0 0 44 45 218.2 0 253.0 255.4 240.7 0 0 0 0 0 0 0 0 1 0 0 45 46 203.7 0 218.2 253.0 255.4 0 0 0 0 0 0 0 0 0 1 0 46 47 205.6 0 203.7 218.2 253.0 0 0 0 0 0 0 0 0 0 0 1 47 48 215.6 0 205.6 203.7 218.2 0 0 0 0 0 0 0 0 0 0 0 48 49 188.5 0 215.6 205.6 203.7 1 0 0 0 0 0 0 0 0 0 0 49 50 202.9 0 188.5 215.6 205.6 0 1 0 0 0 0 0 0 0 0 0 50 51 214.0 0 202.9 188.5 215.6 0 0 1 0 0 0 0 0 0 0 0 51 52 230.3 0 214.0 202.9 188.5 0 0 0 1 0 0 0 0 0 0 0 52 53 230.0 0 230.3 214.0 202.9 0 0 0 0 1 0 0 0 0 0 0 53 54 241.0 0 230.0 230.3 214.0 0 0 0 0 0 1 0 0 0 0 0 54 55 259.6 1 241.0 230.0 230.3 0 0 0 0 0 0 1 0 0 0 0 55 56 247.8 1 259.6 241.0 230.0 0 0 0 0 0 0 0 1 0 0 0 56 57 270.3 1 247.8 259.6 241.0 0 0 0 0 0 0 0 0 1 0 0 57 58 289.7 1 270.3 247.8 259.6 0 0 0 0 0 0 0 0 0 1 0 58 59 322.7 1 289.7 270.3 247.8 0 0 0 0 0 0 0 0 0 0 1 59 60 315.0 1 322.7 289.7 270.3 0 0 0 0 0 0 0 0 0 0 0 60 61 320.2 1 315.0 322.7 289.7 1 0 0 0 0 0 0 0 0 0 0 61 62 329.5 1 320.2 315.0 322.7 0 1 0 0 0 0 0 0 0 0 0 62 63 360.6 1 329.5 320.2 315.0 0 0 1 0 0 0 0 0 0 0 0 63 64 382.2 1 360.6 329.5 320.2 0 0 0 1 0 0 0 0 0 0 0 64 65 435.4 1 382.2 360.6 329.5 0 0 0 0 1 0 0 0 0 0 0 65 66 464.0 1 435.4 382.2 360.6 0 0 0 0 0 1 0 0 0 0 0 66 67 468.8 1 464.0 435.4 382.2 0 0 0 0 0 0 1 0 0 0 0 67 68 403.0 1 468.8 464.0 435.4 0 0 0 0 0 0 0 1 0 0 0 68 69 351.6 1 403.0 468.8 464.0 0 0 0 0 0 0 0 0 1 0 0 69 70 252.0 1 351.6 403.0 468.8 0 0 0 0 0 0 0 0 0 1 0 70 71 188.0 1 252.0 351.6 403.0 0 0 0 0 0 0 0 0 0 0 1 71 72 146.5 1 188.0 252.0 351.6 0 0 0 0 0 0 0 0 0 0 0 72 73 152.9 1 146.5 188.0 252.0 1 0 0 0 0 0 0 0 0 0 0 73 74 148.1 1 152.9 146.5 188.0 0 1 0 0 0 0 0 0 0 0 0 74 75 165.1 1 148.1 152.9 146.5 0 0 1 0 0 0 0 0 0 0 0 75 76 177.0 1 165.1 148.1 152.9 0 0 0 1 0 0 0 0 0 0 0 76 77 206.1 1 177.0 165.1 148.1 0 0 0 0 1 0 0 0 0 0 0 77 78 244.9 1 206.1 177.0 165.1 0 0 0 0 0 1 0 0 0 0 0 78 79 228.6 1 244.9 206.1 177.0 0 0 0 0 0 0 1 0 0 0 0 79 80 253.4 1 228.6 244.9 206.1 0 0 0 0 0 0 0 1 0 0 0 80 81 241.1 1 253.4 228.6 244.9 0 0 0 0 0 0 0 0 1 0 0 81 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) D Y1 Y2 Y3 M1 12.06300 4.10249 1.26715 -0.08184 -0.30610 8.12828 M2 M3 M4 M5 M6 M7 1.82766 8.39072 3.13299 6.63553 8.96124 -0.41942 M8 M9 M10 M11 t -5.56120 -10.17032 -3.70759 1.00729 0.24820 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -47.271 -10.564 1.054 14.138 42.533 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 12.06300 9.99800 1.207 0.2321 D 4.10249 7.90593 0.519 0.6056 Y1 1.26715 0.12015 10.547 1.24e-15 *** Y2 -0.08184 0.20153 -0.406 0.6860 Y3 -0.30610 0.12411 -2.466 0.0163 * M1 8.12828 10.63704 0.764 0.4476 M2 1.82766 10.71927 0.171 0.8652 M3 8.39072 10.82373 0.775 0.4411 M4 3.13299 10.86559 0.288 0.7740 M5 6.63553 10.97206 0.605 0.5475 M6 8.96124 10.93353 0.820 0.4155 M7 -0.41942 11.02992 -0.038 0.9698 M8 -5.56120 11.00080 -0.506 0.6149 M9 -10.17032 10.81095 -0.941 0.3504 M10 -3.70759 11.08029 -0.335 0.7390 M11 1.00729 11.02754 0.091 0.9275 t 0.24820 0.18239 1.361 0.1784 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 19.06 on 64 degrees of freedom Multiple R-squared: 0.9621, Adjusted R-squared: 0.9526 F-statistic: 101.5 on 16 and 64 DF, p-value: < 2.2e-16 > 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,] 1.255494e-01 2.510988e-01 0.8744506 [2,] 4.962979e-02 9.925958e-02 0.9503702 [3,] 2.166906e-02 4.333811e-02 0.9783309 [4,] 2.053025e-02 4.106050e-02 0.9794697 [5,] 8.521634e-03 1.704327e-02 0.9914784 [6,] 3.234566e-03 6.469132e-03 0.9967654 [7,] 1.835594e-03 3.671189e-03 0.9981644 [8,] 3.201580e-03 6.403160e-03 0.9967984 [9,] 1.535532e-03 3.071064e-03 0.9984645 [10,] 1.544244e-03 3.088488e-03 0.9984558 [11,] 3.348590e-03 6.697179e-03 0.9966514 [12,] 1.512819e-03 3.025637e-03 0.9984872 [13,] 1.018165e-03 2.036331e-03 0.9989818 [14,] 4.442311e-04 8.884622e-04 0.9995558 [15,] 6.721527e-04 1.344305e-03 0.9993278 [16,] 3.348096e-04 6.696192e-04 0.9996652 [17,] 1.450170e-04 2.900340e-04 0.9998550 [18,] 6.738141e-05 1.347628e-04 0.9999326 [19,] 6.811581e-05 1.362316e-04 0.9999319 [20,] 3.690293e-05 7.380586e-05 0.9999631 [21,] 6.178587e-05 1.235717e-04 0.9999382 [22,] 3.790298e-05 7.580595e-05 0.9999621 [23,] 3.500385e-05 7.000770e-05 0.9999650 [24,] 2.194252e-05 4.388503e-05 0.9999781 [25,] 2.465947e-05 4.931893e-05 0.9999753 [26,] 1.268840e-04 2.537680e-04 0.9998731 [27,] 1.750714e-04 3.501428e-04 0.9998249 [28,] 1.141319e-04 2.282638e-04 0.9998859 [29,] 9.971797e-05 1.994359e-04 0.9999003 [30,] 1.603178e-03 3.206357e-03 0.9983968 [31,] 1.249825e-03 2.499651e-03 0.9987502 [32,] 8.355854e-04 1.671171e-03 0.9991644 [33,] 6.045635e-04 1.209127e-03 0.9993954 [34,] 3.221877e-04 6.443754e-04 0.9996778 [35,] 1.396249e-04 2.792498e-04 0.9998604 [36,] 5.532833e-05 1.106567e-04 0.9999447 [37,] 2.593172e-04 5.186344e-04 0.9997407 [38,] 1.943270e-03 3.886540e-03 0.9980567 [39,] 9.591264e-04 1.918253e-03 0.9990409 [40,] 9.734355e-04 1.946871e-03 0.9990266 [41,] 4.983528e-04 9.967055e-04 0.9995016 [42,] 1.086304e-02 2.172609e-02 0.9891370 > postscript(file="/var/www/html/rcomp/tmp/1ij9i1260834280.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/html/rcomp/tmp/2vm5f1260834280.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/html/rcomp/tmp/3wfjo1260834280.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/html/rcomp/tmp/4vg9b1260834280.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/html/rcomp/tmp/5nvyf1260834280.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 = 81 Frequency = 1 1 2 3 4 5 6 1.3847238 -0.3867996 -22.8313560 -18.9219554 1.8967224 -0.8021058 7 8 9 10 11 12 -1.4867110 4.4452712 -3.8535361 10.3484113 -2.5278185 -2.3051181 13 14 15 16 17 18 -6.9989942 -7.3388991 -4.7503000 -6.2547998 1.0544721 -20.3402714 19 20 21 22 23 24 4.3169278 17.3714067 1.8552150 18.6063606 -20.3690664 -12.0547295 25 26 27 28 29 30 13.4000110 1.3529799 8.7211270 -10.5644015 -21.3622291 17.0697180 31 32 33 34 35 36 6.5426602 18.7996758 5.9008189 -8.4745166 -2.9313562 15.8291386 37 38 39 40 41 42 16.0102454 -14.4137350 -3.4752536 26.7899679 -8.8140379 -12.1212722 43 44 45 46 47 48 21.4281915 5.7841187 -20.8699370 6.3192783 18.0472756 14.5597087 49 50 51 52 53 54 -37.8712914 18.3208792 5.2058995 5.3331002 -14.0558829 -0.5179315 55 56 57 58 59 60 14.1383326 -15.5286933 31.1739305 20.0799749 21.7635071 -18.5183647 61 62 63 64 65 66 -3.2987251 14.9358059 25.5086084 15.0625964 42.5333415 12.4345832 67 68 69 70 71 72 1.0921665 -47.2712744 -1.7844726 -46.8795085 -13.9825417 2.4893649 73 74 75 76 77 78 17.3740306 -12.4702314 -8.3787253 -11.4445078 -1.2523862 4.2772796 79 80 81 -46.0315676 16.3994953 -12.4220187 > postscript(file="/var/www/html/rcomp/tmp/6rpwm1260834280.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 = 81 Frequency = 1 lag(myerror, k = 1) myerror 0 1.3847238 NA 1 -0.3867996 1.3847238 2 -22.8313560 -0.3867996 3 -18.9219554 -22.8313560 4 1.8967224 -18.9219554 5 -0.8021058 1.8967224 6 -1.4867110 -0.8021058 7 4.4452712 -1.4867110 8 -3.8535361 4.4452712 9 10.3484113 -3.8535361 10 -2.5278185 10.3484113 11 -2.3051181 -2.5278185 12 -6.9989942 -2.3051181 13 -7.3388991 -6.9989942 14 -4.7503000 -7.3388991 15 -6.2547998 -4.7503000 16 1.0544721 -6.2547998 17 -20.3402714 1.0544721 18 4.3169278 -20.3402714 19 17.3714067 4.3169278 20 1.8552150 17.3714067 21 18.6063606 1.8552150 22 -20.3690664 18.6063606 23 -12.0547295 -20.3690664 24 13.4000110 -12.0547295 25 1.3529799 13.4000110 26 8.7211270 1.3529799 27 -10.5644015 8.7211270 28 -21.3622291 -10.5644015 29 17.0697180 -21.3622291 30 6.5426602 17.0697180 31 18.7996758 6.5426602 32 5.9008189 18.7996758 33 -8.4745166 5.9008189 34 -2.9313562 -8.4745166 35 15.8291386 -2.9313562 36 16.0102454 15.8291386 37 -14.4137350 16.0102454 38 -3.4752536 -14.4137350 39 26.7899679 -3.4752536 40 -8.8140379 26.7899679 41 -12.1212722 -8.8140379 42 21.4281915 -12.1212722 43 5.7841187 21.4281915 44 -20.8699370 5.7841187 45 6.3192783 -20.8699370 46 18.0472756 6.3192783 47 14.5597087 18.0472756 48 -37.8712914 14.5597087 49 18.3208792 -37.8712914 50 5.2058995 18.3208792 51 5.3331002 5.2058995 52 -14.0558829 5.3331002 53 -0.5179315 -14.0558829 54 14.1383326 -0.5179315 55 -15.5286933 14.1383326 56 31.1739305 -15.5286933 57 20.0799749 31.1739305 58 21.7635071 20.0799749 59 -18.5183647 21.7635071 60 -3.2987251 -18.5183647 61 14.9358059 -3.2987251 62 25.5086084 14.9358059 63 15.0625964 25.5086084 64 42.5333415 15.0625964 65 12.4345832 42.5333415 66 1.0921665 12.4345832 67 -47.2712744 1.0921665 68 -1.7844726 -47.2712744 69 -46.8795085 -1.7844726 70 -13.9825417 -46.8795085 71 2.4893649 -13.9825417 72 17.3740306 2.4893649 73 -12.4702314 17.3740306 74 -8.3787253 -12.4702314 75 -11.4445078 -8.3787253 76 -1.2523862 -11.4445078 77 4.2772796 -1.2523862 78 -46.0315676 4.2772796 79 16.3994953 -46.0315676 80 -12.4220187 16.3994953 81 NA -12.4220187 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -0.3867996 1.3847238 [2,] -22.8313560 -0.3867996 [3,] -18.9219554 -22.8313560 [4,] 1.8967224 -18.9219554 [5,] -0.8021058 1.8967224 [6,] -1.4867110 -0.8021058 [7,] 4.4452712 -1.4867110 [8,] -3.8535361 4.4452712 [9,] 10.3484113 -3.8535361 [10,] -2.5278185 10.3484113 [11,] -2.3051181 -2.5278185 [12,] -6.9989942 -2.3051181 [13,] -7.3388991 -6.9989942 [14,] -4.7503000 -7.3388991 [15,] -6.2547998 -4.7503000 [16,] 1.0544721 -6.2547998 [17,] -20.3402714 1.0544721 [18,] 4.3169278 -20.3402714 [19,] 17.3714067 4.3169278 [20,] 1.8552150 17.3714067 [21,] 18.6063606 1.8552150 [22,] -20.3690664 18.6063606 [23,] -12.0547295 -20.3690664 [24,] 13.4000110 -12.0547295 [25,] 1.3529799 13.4000110 [26,] 8.7211270 1.3529799 [27,] -10.5644015 8.7211270 [28,] -21.3622291 -10.5644015 [29,] 17.0697180 -21.3622291 [30,] 6.5426602 17.0697180 [31,] 18.7996758 6.5426602 [32,] 5.9008189 18.7996758 [33,] -8.4745166 5.9008189 [34,] -2.9313562 -8.4745166 [35,] 15.8291386 -2.9313562 [36,] 16.0102454 15.8291386 [37,] -14.4137350 16.0102454 [38,] -3.4752536 -14.4137350 [39,] 26.7899679 -3.4752536 [40,] -8.8140379 26.7899679 [41,] -12.1212722 -8.8140379 [42,] 21.4281915 -12.1212722 [43,] 5.7841187 21.4281915 [44,] -20.8699370 5.7841187 [45,] 6.3192783 -20.8699370 [46,] 18.0472756 6.3192783 [47,] 14.5597087 18.0472756 [48,] -37.8712914 14.5597087 [49,] 18.3208792 -37.8712914 [50,] 5.2058995 18.3208792 [51,] 5.3331002 5.2058995 [52,] -14.0558829 5.3331002 [53,] -0.5179315 -14.0558829 [54,] 14.1383326 -0.5179315 [55,] -15.5286933 14.1383326 [56,] 31.1739305 -15.5286933 [57,] 20.0799749 31.1739305 [58,] 21.7635071 20.0799749 [59,] -18.5183647 21.7635071 [60,] -3.2987251 -18.5183647 [61,] 14.9358059 -3.2987251 [62,] 25.5086084 14.9358059 [63,] 15.0625964 25.5086084 [64,] 42.5333415 15.0625964 [65,] 12.4345832 42.5333415 [66,] 1.0921665 12.4345832 [67,] -47.2712744 1.0921665 [68,] -1.7844726 -47.2712744 [69,] -46.8795085 -1.7844726 [70,] -13.9825417 -46.8795085 [71,] 2.4893649 -13.9825417 [72,] 17.3740306 2.4893649 [73,] -12.4702314 17.3740306 [74,] -8.3787253 -12.4702314 [75,] -11.4445078 -8.3787253 [76,] -1.2523862 -11.4445078 [77,] 4.2772796 -1.2523862 [78,] -46.0315676 4.2772796 [79,] 16.3994953 -46.0315676 [80,] -12.4220187 16.3994953 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -0.3867996 1.3847238 2 -22.8313560 -0.3867996 3 -18.9219554 -22.8313560 4 1.8967224 -18.9219554 5 -0.8021058 1.8967224 6 -1.4867110 -0.8021058 7 4.4452712 -1.4867110 8 -3.8535361 4.4452712 9 10.3484113 -3.8535361 10 -2.5278185 10.3484113 11 -2.3051181 -2.5278185 12 -6.9989942 -2.3051181 13 -7.3388991 -6.9989942 14 -4.7503000 -7.3388991 15 -6.2547998 -4.7503000 16 1.0544721 -6.2547998 17 -20.3402714 1.0544721 18 4.3169278 -20.3402714 19 17.3714067 4.3169278 20 1.8552150 17.3714067 21 18.6063606 1.8552150 22 -20.3690664 18.6063606 23 -12.0547295 -20.3690664 24 13.4000110 -12.0547295 25 1.3529799 13.4000110 26 8.7211270 1.3529799 27 -10.5644015 8.7211270 28 -21.3622291 -10.5644015 29 17.0697180 -21.3622291 30 6.5426602 17.0697180 31 18.7996758 6.5426602 32 5.9008189 18.7996758 33 -8.4745166 5.9008189 34 -2.9313562 -8.4745166 35 15.8291386 -2.9313562 36 16.0102454 15.8291386 37 -14.4137350 16.0102454 38 -3.4752536 -14.4137350 39 26.7899679 -3.4752536 40 -8.8140379 26.7899679 41 -12.1212722 -8.8140379 42 21.4281915 -12.1212722 43 5.7841187 21.4281915 44 -20.8699370 5.7841187 45 6.3192783 -20.8699370 46 18.0472756 6.3192783 47 14.5597087 18.0472756 48 -37.8712914 14.5597087 49 18.3208792 -37.8712914 50 5.2058995 18.3208792 51 5.3331002 5.2058995 52 -14.0558829 5.3331002 53 -0.5179315 -14.0558829 54 14.1383326 -0.5179315 55 -15.5286933 14.1383326 56 31.1739305 -15.5286933 57 20.0799749 31.1739305 58 21.7635071 20.0799749 59 -18.5183647 21.7635071 60 -3.2987251 -18.5183647 61 14.9358059 -3.2987251 62 25.5086084 14.9358059 63 15.0625964 25.5086084 64 42.5333415 15.0625964 65 12.4345832 42.5333415 66 1.0921665 12.4345832 67 -47.2712744 1.0921665 68 -1.7844726 -47.2712744 69 -46.8795085 -1.7844726 70 -13.9825417 -46.8795085 71 2.4893649 -13.9825417 72 17.3740306 2.4893649 73 -12.4702314 17.3740306 74 -8.3787253 -12.4702314 75 -11.4445078 -8.3787253 76 -1.2523862 -11.4445078 77 4.2772796 -1.2523862 78 -46.0315676 4.2772796 79 16.3994953 -46.0315676 80 -12.4220187 16.3994953 > 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/7363e1260834281.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/html/rcomp/tmp/8mnh01260834281.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/html/rcomp/tmp/9egev1260834281.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/html/rcomp/tmp/102kwe1260834281.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/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/11ih8q1260834281.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/12op3z1260834281.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/13pzku1260834281.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/14hbs91260834281.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/15skxu1260834281.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/16fi3r1260834281.tab") + } > > try(system("convert tmp/1ij9i1260834280.ps tmp/1ij9i1260834280.png",intern=TRUE)) character(0) > try(system("convert tmp/2vm5f1260834280.ps tmp/2vm5f1260834280.png",intern=TRUE)) character(0) > try(system("convert tmp/3wfjo1260834280.ps tmp/3wfjo1260834280.png",intern=TRUE)) character(0) > try(system("convert tmp/4vg9b1260834280.ps tmp/4vg9b1260834280.png",intern=TRUE)) character(0) > try(system("convert tmp/5nvyf1260834280.ps tmp/5nvyf1260834280.png",intern=TRUE)) character(0) > try(system("convert tmp/6rpwm1260834280.ps tmp/6rpwm1260834280.png",intern=TRUE)) character(0) > try(system("convert tmp/7363e1260834281.ps tmp/7363e1260834281.png",intern=TRUE)) character(0) > try(system("convert tmp/8mnh01260834281.ps tmp/8mnh01260834281.png",intern=TRUE)) character(0) > try(system("convert tmp/9egev1260834281.ps tmp/9egev1260834281.png",intern=TRUE)) character(0) > try(system("convert tmp/102kwe1260834281.ps tmp/102kwe1260834281.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.687 1.571 3.487