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(4716.99 + ,0 + ,4926.65 + ,0 + ,4920.10 + ,0 + ,5170.09 + ,0 + ,5246.24 + ,0 + ,5283.61 + ,0 + ,4979.05 + ,0 + ,4825.20 + ,0 + ,4695.12 + ,0 + ,4711.54 + ,0 + ,4727.22 + ,0 + ,4384.96 + ,0 + ,4378.75 + ,0 + ,4472.93 + ,0 + ,4564.07 + ,0 + ,4310.54 + ,0 + ,4171.38 + ,0 + ,4049.38 + ,0 + ,3591.37 + ,0 + ,3720.46 + ,0 + ,4107.23 + ,0 + ,4101.71 + ,0 + ,4162.34 + ,0 + ,4136.22 + ,0 + ,4125.88 + ,0 + ,4031.48 + ,0 + ,3761.36 + ,0 + ,3408.56 + ,0 + ,3228.47 + ,0 + ,3090.45 + ,0 + ,2741.14 + ,0 + ,2980.44 + ,0 + ,3104.33 + ,0 + ,3181.57 + ,0 + ,2863.86 + ,0 + ,2898.01 + ,0 + ,3112.33 + ,0 + ,3254.33 + ,0 + ,3513.47 + ,0 + ,3587.61 + ,0 + ,3727.45 + ,0 + ,3793.34 + ,0 + ,3817.58 + ,0 + ,3845.13 + ,0 + ,3931.86 + ,0 + ,4197.52 + ,0 + ,4307.13 + ,0 + ,4229.43 + ,0 + ,4362.28 + ,0 + ,4217.34 + ,0 + ,4361.28 + ,0 + ,4327.74 + ,0 + ,4417.65 + ,0 + ,4557.68 + ,0 + ,4650.35 + ,0 + ,4967.18 + ,0 + ,5123.42 + ,0 + ,5290.85 + ,0 + ,5535.66 + ,0 + ,5514.06 + ,0 + ,5493.88 + ,0 + ,5694.83 + ,0 + ,5850.41 + ,0 + ,6116.64 + ,0 + ,6175.00 + ,0 + ,6513.58 + ,0 + ,6383.78 + ,0 + ,6673.66 + ,0 + ,6936.61 + ,0 + ,7300.68 + ,0 + ,7392.93 + ,0 + ,7497.31 + ,0 + ,7584.71 + ,0 + ,7160.79 + ,0 + ,7196.19 + ,0 + ,7245.63 + ,0 + ,7347.51 + ,0 + ,7425.75 + ,0 + ,7778.51 + ,0 + ,7822.33 + ,0 + ,8181.22 + ,0 + ,8371.47 + ,0 + ,8347.71 + ,0 + ,8672.11 + ,0 + ,8802.79 + ,0 + ,9138.46 + ,0 + ,9123.29 + ,0 + ,9023.21 + ,1 + ,8850.41 + ,1 + ,8864.58 + ,1 + ,9163.74 + ,1 + ,8516.66 + ,1 + ,8553.44 + ,1 + ,7555.20 + ,1 + ,7851.22 + ,1 + ,7442.00 + ,1 + ,7992.53 + ,1 + ,8264.04 + ,1 + ,7517.39 + ,1 + ,7200.40 + ,1 + ,7193.69 + ,1 + ,6193.58 + ,1 + ,5104.21 + ,1 + ,4800.46 + ,1 + ,4461.61 + ,1 + ,4398.59 + ,1 + ,4243.63 + ,1 + ,4293.82 + ,1) + ,dim=c(2 + ,108) + ,dimnames=list(c('Y' + ,'X') + ,1:108)) > y <- array(NA,dim=c(2,108),dimnames=list(c('Y','X'),1:108)) > 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 = '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 X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 1 4716.99 0 1 0 0 0 0 0 0 0 0 0 0 2 4926.65 0 0 1 0 0 0 0 0 0 0 0 0 3 4920.10 0 0 0 1 0 0 0 0 0 0 0 0 4 5170.09 0 0 0 0 1 0 0 0 0 0 0 0 5 5246.24 0 0 0 0 0 1 0 0 0 0 0 0 6 5283.61 0 0 0 0 0 0 1 0 0 0 0 0 7 4979.05 0 0 0 0 0 0 0 1 0 0 0 0 8 4825.20 0 0 0 0 0 0 0 0 1 0 0 0 9 4695.12 0 0 0 0 0 0 0 0 0 1 0 0 10 4711.54 0 0 0 0 0 0 0 0 0 0 1 0 11 4727.22 0 0 0 0 0 0 0 0 0 0 0 1 12 4384.96 0 0 0 0 0 0 0 0 0 0 0 0 13 4378.75 0 1 0 0 0 0 0 0 0 0 0 0 14 4472.93 0 0 1 0 0 0 0 0 0 0 0 0 15 4564.07 0 0 0 1 0 0 0 0 0 0 0 0 16 4310.54 0 0 0 0 1 0 0 0 0 0 0 0 17 4171.38 0 0 0 0 0 1 0 0 0 0 0 0 18 4049.38 0 0 0 0 0 0 1 0 0 0 0 0 19 3591.37 0 0 0 0 0 0 0 1 0 0 0 0 20 3720.46 0 0 0 0 0 0 0 0 1 0 0 0 21 4107.23 0 0 0 0 0 0 0 0 0 1 0 0 22 4101.71 0 0 0 0 0 0 0 0 0 0 1 0 23 4162.34 0 0 0 0 0 0 0 0 0 0 0 1 24 4136.22 0 0 0 0 0 0 0 0 0 0 0 0 25 4125.88 0 1 0 0 0 0 0 0 0 0 0 0 26 4031.48 0 0 1 0 0 0 0 0 0 0 0 0 27 3761.36 0 0 0 1 0 0 0 0 0 0 0 0 28 3408.56 0 0 0 0 1 0 0 0 0 0 0 0 29 3228.47 0 0 0 0 0 1 0 0 0 0 0 0 30 3090.45 0 0 0 0 0 0 1 0 0 0 0 0 31 2741.14 0 0 0 0 0 0 0 1 0 0 0 0 32 2980.44 0 0 0 0 0 0 0 0 1 0 0 0 33 3104.33 0 0 0 0 0 0 0 0 0 1 0 0 34 3181.57 0 0 0 0 0 0 0 0 0 0 1 0 35 2863.86 0 0 0 0 0 0 0 0 0 0 0 1 36 2898.01 0 0 0 0 0 0 0 0 0 0 0 0 37 3112.33 0 1 0 0 0 0 0 0 0 0 0 0 38 3254.33 0 0 1 0 0 0 0 0 0 0 0 0 39 3513.47 0 0 0 1 0 0 0 0 0 0 0 0 40 3587.61 0 0 0 0 1 0 0 0 0 0 0 0 41 3727.45 0 0 0 0 0 1 0 0 0 0 0 0 42 3793.34 0 0 0 0 0 0 1 0 0 0 0 0 43 3817.58 0 0 0 0 0 0 0 1 0 0 0 0 44 3845.13 0 0 0 0 0 0 0 0 1 0 0 0 45 3931.86 0 0 0 0 0 0 0 0 0 1 0 0 46 4197.52 0 0 0 0 0 0 0 0 0 0 1 0 47 4307.13 0 0 0 0 0 0 0 0 0 0 0 1 48 4229.43 0 0 0 0 0 0 0 0 0 0 0 0 49 4362.28 0 1 0 0 0 0 0 0 0 0 0 0 50 4217.34 0 0 1 0 0 0 0 0 0 0 0 0 51 4361.28 0 0 0 1 0 0 0 0 0 0 0 0 52 4327.74 0 0 0 0 1 0 0 0 0 0 0 0 53 4417.65 0 0 0 0 0 1 0 0 0 0 0 0 54 4557.68 0 0 0 0 0 0 1 0 0 0 0 0 55 4650.35 0 0 0 0 0 0 0 1 0 0 0 0 56 4967.18 0 0 0 0 0 0 0 0 1 0 0 0 57 5123.42 0 0 0 0 0 0 0 0 0 1 0 0 58 5290.85 0 0 0 0 0 0 0 0 0 0 1 0 59 5535.66 0 0 0 0 0 0 0 0 0 0 0 1 60 5514.06 0 0 0 0 0 0 0 0 0 0 0 0 61 5493.88 0 1 0 0 0 0 0 0 0 0 0 0 62 5694.83 0 0 1 0 0 0 0 0 0 0 0 0 63 5850.41 0 0 0 1 0 0 0 0 0 0 0 0 64 6116.64 0 0 0 0 1 0 0 0 0 0 0 0 65 6175.00 0 0 0 0 0 1 0 0 0 0 0 0 66 6513.58 0 0 0 0 0 0 1 0 0 0 0 0 67 6383.78 0 0 0 0 0 0 0 1 0 0 0 0 68 6673.66 0 0 0 0 0 0 0 0 1 0 0 0 69 6936.61 0 0 0 0 0 0 0 0 0 1 0 0 70 7300.68 0 0 0 0 0 0 0 0 0 0 1 0 71 7392.93 0 0 0 0 0 0 0 0 0 0 0 1 72 7497.31 0 0 0 0 0 0 0 0 0 0 0 0 73 7584.71 0 1 0 0 0 0 0 0 0 0 0 0 74 7160.79 0 0 1 0 0 0 0 0 0 0 0 0 75 7196.19 0 0 0 1 0 0 0 0 0 0 0 0 76 7245.63 0 0 0 0 1 0 0 0 0 0 0 0 77 7347.51 0 0 0 0 0 1 0 0 0 0 0 0 78 7425.75 0 0 0 0 0 0 1 0 0 0 0 0 79 7778.51 0 0 0 0 0 0 0 1 0 0 0 0 80 7822.33 0 0 0 0 0 0 0 0 1 0 0 0 81 8181.22 0 0 0 0 0 0 0 0 0 1 0 0 82 8371.47 0 0 0 0 0 0 0 0 0 0 1 0 83 8347.71 0 0 0 0 0 0 0 0 0 0 0 1 84 8672.11 0 0 0 0 0 0 0 0 0 0 0 0 85 8802.79 0 1 0 0 0 0 0 0 0 0 0 0 86 9138.46 0 0 1 0 0 0 0 0 0 0 0 0 87 9123.29 0 0 0 1 0 0 0 0 0 0 0 0 88 9023.21 1 0 0 0 1 0 0 0 0 0 0 0 89 8850.41 1 0 0 0 0 1 0 0 0 0 0 0 90 8864.58 1 0 0 0 0 0 1 0 0 0 0 0 91 9163.74 1 0 0 0 0 0 0 1 0 0 0 0 92 8516.66 1 0 0 0 0 0 0 0 1 0 0 0 93 8553.44 1 0 0 0 0 0 0 0 0 1 0 0 94 7555.20 1 0 0 0 0 0 0 0 0 0 1 0 95 7851.22 1 0 0 0 0 0 0 0 0 0 0 1 96 7442.00 1 0 0 0 0 0 0 0 0 0 0 0 97 7992.53 1 1 0 0 0 0 0 0 0 0 0 0 98 8264.04 1 0 1 0 0 0 0 0 0 0 0 0 99 7517.39 1 0 0 1 0 0 0 0 0 0 0 0 100 7200.40 1 0 0 0 1 0 0 0 0 0 0 0 101 7193.69 1 0 0 0 0 1 0 0 0 0 0 0 102 6193.58 1 0 0 0 0 0 1 0 0 0 0 0 103 5104.21 1 0 0 0 0 0 0 1 0 0 0 0 104 4800.46 1 0 0 0 0 0 0 0 1 0 0 0 105 4461.61 1 0 0 0 0 0 0 0 0 1 0 0 106 4398.59 1 0 0 0 0 0 0 0 0 0 1 0 107 4243.63 1 0 0 0 0 0 0 0 0 0 0 1 108 4293.82 1 0 0 0 0 0 0 0 0 0 0 0 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) X M1 M2 M3 M4 5025.841 1917.674 379.988 445.623 406.368 146.944 M5 M6 M7 M8 M9 M10 143.320 78.226 -95.354 -101.822 2.991 4.579 M11 40.420 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -2740.3 -1135.6 -336.4 1420.5 3691.1 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5025.841 602.185 8.346 5.59e-13 *** X 1917.674 436.729 4.391 2.93e-05 *** M1 379.988 841.885 0.451 0.653 M2 445.623 841.885 0.529 0.598 M3 406.368 841.885 0.483 0.630 M4 146.944 840.486 0.175 0.862 M5 143.320 840.486 0.171 0.865 M6 78.226 840.486 0.093 0.926 M7 -95.354 840.486 -0.113 0.910 M8 -101.822 840.486 -0.121 0.904 M9 2.991 840.486 0.004 0.997 M10 4.579 840.486 0.005 0.996 M11 40.420 840.486 0.048 0.962 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1783 on 95 degrees of freedom Multiple R-squared: 0.1715, Adjusted R-squared: 0.0668 F-statistic: 1.638 on 12 and 95 DF, p-value: 0.0939 > 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.614937e-02 3.229874e-02 0.9838506 [2,] 1.109204e-02 2.218408e-02 0.9889080 [3,] 8.481777e-03 1.696355e-02 0.9915182 [4,] 7.215770e-03 1.443154e-02 0.9927842 [5,] 4.180421e-03 8.360841e-03 0.9958196 [6,] 1.613716e-03 3.227431e-03 0.9983863 [7,] 6.198449e-04 1.239690e-03 0.9993802 [8,] 2.276936e-04 4.553873e-04 0.9997723 [9,] 7.047087e-05 1.409417e-04 0.9999295 [10,] 2.440841e-05 4.881682e-05 0.9999756 [11,] 1.078205e-05 2.156411e-05 0.9999892 [12,] 7.622735e-06 1.524547e-05 0.9999924 [13,] 9.815857e-06 1.963171e-05 0.9999902 [14,] 1.506068e-05 3.012136e-05 0.9999849 [15,] 2.410588e-05 4.821175e-05 0.9999759 [16,] 3.444737e-05 6.889474e-05 0.9999656 [17,] 3.312993e-05 6.625986e-05 0.9999669 [18,] 3.191483e-05 6.382966e-05 0.9999681 [19,] 2.836349e-05 5.672698e-05 0.9999716 [20,] 3.972086e-05 7.944172e-05 0.9999603 [21,] 4.369560e-05 8.739120e-05 0.9999563 [22,] 5.504887e-05 1.100977e-04 0.9999450 [23,] 6.546111e-05 1.309222e-04 0.9999345 [24,] 5.821375e-05 1.164275e-04 0.9999418 [25,] 4.451170e-05 8.902339e-05 0.9999555 [26,] 3.091048e-05 6.182097e-05 0.9999691 [27,] 2.027045e-05 4.054090e-05 0.9999797 [28,] 1.256400e-05 2.512800e-05 0.9999874 [29,] 7.728095e-06 1.545619e-05 0.9999923 [30,] 4.853522e-06 9.707044e-06 0.9999951 [31,] 2.864311e-06 5.728622e-06 0.9999971 [32,] 1.786961e-06 3.573921e-06 0.9999982 [33,] 1.179374e-06 2.358747e-06 0.9999988 [34,] 9.941410e-07 1.988282e-06 0.9999990 [35,] 9.638247e-07 1.927649e-06 0.9999990 [36,] 8.800055e-07 1.760011e-06 0.9999991 [37,] 8.328729e-07 1.665746e-06 0.9999992 [38,] 8.417348e-07 1.683470e-06 0.9999992 [39,] 8.525450e-07 1.705090e-06 0.9999991 [40,] 1.052518e-06 2.105037e-06 0.9999989 [41,] 1.382809e-06 2.765618e-06 0.9999986 [42,] 1.862547e-06 3.725093e-06 0.9999981 [43,] 2.424763e-06 4.849527e-06 0.9999976 [44,] 4.007388e-06 8.014775e-06 0.9999960 [45,] 7.108015e-06 1.421603e-05 0.9999929 [46,] 1.697566e-05 3.395133e-05 0.9999830 [47,] 4.400232e-05 8.800464e-05 0.9999560 [48,] 9.823973e-05 1.964795e-04 0.9999018 [49,] 2.556371e-04 5.112742e-04 0.9997444 [50,] 6.252391e-04 1.250478e-03 0.9993748 [51,] 1.439048e-03 2.878097e-03 0.9985610 [52,] 3.036970e-03 6.073939e-03 0.9969630 [53,] 5.490076e-03 1.098015e-02 0.9945099 [54,] 9.174658e-03 1.834932e-02 0.9908253 [55,] 1.500954e-02 3.001908e-02 0.9849905 [56,] 2.236628e-02 4.473256e-02 0.9776337 [57,] 3.256660e-02 6.513320e-02 0.9674334 [58,] 4.755542e-02 9.511085e-02 0.9524446 [59,] 6.487974e-02 1.297595e-01 0.9351203 [60,] 7.702408e-02 1.540482e-01 0.9229759 [61,] 9.771959e-02 1.954392e-01 0.9022804 [62,] 1.220944e-01 2.441888e-01 0.8779056 [63,] 1.366604e-01 2.733207e-01 0.8633396 [64,] 1.486575e-01 2.973149e-01 0.8513425 [65,] 1.475172e-01 2.950344e-01 0.8524828 [66,] 1.435554e-01 2.871108e-01 0.8564446 [67,] 1.379602e-01 2.759204e-01 0.8620398 [68,] 1.275342e-01 2.550683e-01 0.8724658 [69,] 1.316190e-01 2.632380e-01 0.8683810 [70,] 1.190890e-01 2.381780e-01 0.8809110 [71,] 1.083178e-01 2.166355e-01 0.8916822 [72,] 9.164158e-02 1.832832e-01 0.9083584 [73,] 6.553962e-02 1.310792e-01 0.9344604 [74,] 4.315376e-02 8.630753e-02 0.9568462 [75,] 3.365026e-02 6.730052e-02 0.9663497 [76,] 4.390633e-02 8.781267e-02 0.9560937 [77,] 5.003816e-02 1.000763e-01 0.9499618 > postscript(file="/var/www/html/rcomp/tmp/1pmv21258647958.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/2f0rb1258647958.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/30f9m1258647958.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/4kacr1258647958.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/5cxv61258647958.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 = 108 Frequency = 1 1 2 3 4 5 6 -688.839541 -544.813985 -512.109541 -2.695748 77.078696 179.543141 7 8 9 10 11 12 48.563141 -98.819081 -333.712415 -318.880193 -339.041304 -640.881304 13 14 15 16 17 18 -1027.079541 -998.533985 -868.139541 -862.245748 -997.781304 -1054.686859 19 20 21 22 23 24 -1339.116859 -1203.559081 -921.602415 -928.710193 -903.921304 -889.621304 25 26 27 28 29 30 -1279.949541 -1439.983985 -1670.849541 -1764.225748 -1940.691304 -2013.616859 31 32 33 34 35 36 -2189.346859 -1943.579081 -1924.502415 -1848.850193 -2202.401304 -2127.831304 37 38 39 40 41 42 -2293.499541 -2217.133985 -1918.739541 -1585.175748 -1441.711304 -1310.726859 43 44 45 46 47 48 -1112.906859 -1078.889081 -1096.972415 -832.900193 -759.131304 -796.411304 49 50 51 52 53 54 -1043.549541 -1254.123985 -1070.929541 -845.045748 -751.511304 -546.386859 55 56 57 58 59 60 -280.136859 43.160919 94.587585 260.429807 469.398696 488.218696 61 62 63 64 65 66 88.050459 223.366015 418.200459 943.854252 1005.838696 1409.513141 67 68 69 70 71 72 1453.293141 1749.640919 1907.777585 2270.259807 2326.668696 2471.468696 73 74 75 76 77 78 2178.880459 1689.326015 1763.980459 2072.844252 2178.348696 2321.683141 79 80 81 82 83 84 2848.023141 2898.310919 3152.387585 3341.049807 3281.448696 3646.268696 85 86 87 88 89 90 3396.960459 3666.996015 3691.080459 1932.750119 1763.574563 1842.839007 91 92 93 94 95 96 2315.579007 1674.966785 1606.933452 607.105674 867.284563 498.484563 97 98 99 100 101 102 669.026326 874.901881 167.506326 109.940119 106.854563 -828.160993 103 104 105 106 107 108 -1743.950993 -2041.233215 -2484.896548 -2549.504326 -2740.305437 -2649.695437 > postscript(file="/var/www/html/rcomp/tmp/6gj2d1258647958.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 = 108 Frequency = 1 lag(myerror, k = 1) myerror 0 -688.839541 NA 1 -544.813985 -688.839541 2 -512.109541 -544.813985 3 -2.695748 -512.109541 4 77.078696 -2.695748 5 179.543141 77.078696 6 48.563141 179.543141 7 -98.819081 48.563141 8 -333.712415 -98.819081 9 -318.880193 -333.712415 10 -339.041304 -318.880193 11 -640.881304 -339.041304 12 -1027.079541 -640.881304 13 -998.533985 -1027.079541 14 -868.139541 -998.533985 15 -862.245748 -868.139541 16 -997.781304 -862.245748 17 -1054.686859 -997.781304 18 -1339.116859 -1054.686859 19 -1203.559081 -1339.116859 20 -921.602415 -1203.559081 21 -928.710193 -921.602415 22 -903.921304 -928.710193 23 -889.621304 -903.921304 24 -1279.949541 -889.621304 25 -1439.983985 -1279.949541 26 -1670.849541 -1439.983985 27 -1764.225748 -1670.849541 28 -1940.691304 -1764.225748 29 -2013.616859 -1940.691304 30 -2189.346859 -2013.616859 31 -1943.579081 -2189.346859 32 -1924.502415 -1943.579081 33 -1848.850193 -1924.502415 34 -2202.401304 -1848.850193 35 -2127.831304 -2202.401304 36 -2293.499541 -2127.831304 37 -2217.133985 -2293.499541 38 -1918.739541 -2217.133985 39 -1585.175748 -1918.739541 40 -1441.711304 -1585.175748 41 -1310.726859 -1441.711304 42 -1112.906859 -1310.726859 43 -1078.889081 -1112.906859 44 -1096.972415 -1078.889081 45 -832.900193 -1096.972415 46 -759.131304 -832.900193 47 -796.411304 -759.131304 48 -1043.549541 -796.411304 49 -1254.123985 -1043.549541 50 -1070.929541 -1254.123985 51 -845.045748 -1070.929541 52 -751.511304 -845.045748 53 -546.386859 -751.511304 54 -280.136859 -546.386859 55 43.160919 -280.136859 56 94.587585 43.160919 57 260.429807 94.587585 58 469.398696 260.429807 59 488.218696 469.398696 60 88.050459 488.218696 61 223.366015 88.050459 62 418.200459 223.366015 63 943.854252 418.200459 64 1005.838696 943.854252 65 1409.513141 1005.838696 66 1453.293141 1409.513141 67 1749.640919 1453.293141 68 1907.777585 1749.640919 69 2270.259807 1907.777585 70 2326.668696 2270.259807 71 2471.468696 2326.668696 72 2178.880459 2471.468696 73 1689.326015 2178.880459 74 1763.980459 1689.326015 75 2072.844252 1763.980459 76 2178.348696 2072.844252 77 2321.683141 2178.348696 78 2848.023141 2321.683141 79 2898.310919 2848.023141 80 3152.387585 2898.310919 81 3341.049807 3152.387585 82 3281.448696 3341.049807 83 3646.268696 3281.448696 84 3396.960459 3646.268696 85 3666.996015 3396.960459 86 3691.080459 3666.996015 87 1932.750119 3691.080459 88 1763.574563 1932.750119 89 1842.839007 1763.574563 90 2315.579007 1842.839007 91 1674.966785 2315.579007 92 1606.933452 1674.966785 93 607.105674 1606.933452 94 867.284563 607.105674 95 498.484563 867.284563 96 669.026326 498.484563 97 874.901881 669.026326 98 167.506326 874.901881 99 109.940119 167.506326 100 106.854563 109.940119 101 -828.160993 106.854563 102 -1743.950993 -828.160993 103 -2041.233215 -1743.950993 104 -2484.896548 -2041.233215 105 -2549.504326 -2484.896548 106 -2740.305437 -2549.504326 107 -2649.695437 -2740.305437 108 NA -2649.695437 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -544.813985 -688.839541 [2,] -512.109541 -544.813985 [3,] -2.695748 -512.109541 [4,] 77.078696 -2.695748 [5,] 179.543141 77.078696 [6,] 48.563141 179.543141 [7,] -98.819081 48.563141 [8,] -333.712415 -98.819081 [9,] -318.880193 -333.712415 [10,] -339.041304 -318.880193 [11,] -640.881304 -339.041304 [12,] -1027.079541 -640.881304 [13,] -998.533985 -1027.079541 [14,] -868.139541 -998.533985 [15,] -862.245748 -868.139541 [16,] -997.781304 -862.245748 [17,] -1054.686859 -997.781304 [18,] -1339.116859 -1054.686859 [19,] -1203.559081 -1339.116859 [20,] -921.602415 -1203.559081 [21,] -928.710193 -921.602415 [22,] -903.921304 -928.710193 [23,] -889.621304 -903.921304 [24,] -1279.949541 -889.621304 [25,] -1439.983985 -1279.949541 [26,] -1670.849541 -1439.983985 [27,] -1764.225748 -1670.849541 [28,] -1940.691304 -1764.225748 [29,] -2013.616859 -1940.691304 [30,] -2189.346859 -2013.616859 [31,] -1943.579081 -2189.346859 [32,] -1924.502415 -1943.579081 [33,] -1848.850193 -1924.502415 [34,] -2202.401304 -1848.850193 [35,] -2127.831304 -2202.401304 [36,] -2293.499541 -2127.831304 [37,] -2217.133985 -2293.499541 [38,] -1918.739541 -2217.133985 [39,] -1585.175748 -1918.739541 [40,] -1441.711304 -1585.175748 [41,] -1310.726859 -1441.711304 [42,] -1112.906859 -1310.726859 [43,] -1078.889081 -1112.906859 [44,] -1096.972415 -1078.889081 [45,] -832.900193 -1096.972415 [46,] -759.131304 -832.900193 [47,] -796.411304 -759.131304 [48,] -1043.549541 -796.411304 [49,] -1254.123985 -1043.549541 [50,] -1070.929541 -1254.123985 [51,] -845.045748 -1070.929541 [52,] -751.511304 -845.045748 [53,] -546.386859 -751.511304 [54,] -280.136859 -546.386859 [55,] 43.160919 -280.136859 [56,] 94.587585 43.160919 [57,] 260.429807 94.587585 [58,] 469.398696 260.429807 [59,] 488.218696 469.398696 [60,] 88.050459 488.218696 [61,] 223.366015 88.050459 [62,] 418.200459 223.366015 [63,] 943.854252 418.200459 [64,] 1005.838696 943.854252 [65,] 1409.513141 1005.838696 [66,] 1453.293141 1409.513141 [67,] 1749.640919 1453.293141 [68,] 1907.777585 1749.640919 [69,] 2270.259807 1907.777585 [70,] 2326.668696 2270.259807 [71,] 2471.468696 2326.668696 [72,] 2178.880459 2471.468696 [73,] 1689.326015 2178.880459 [74,] 1763.980459 1689.326015 [75,] 2072.844252 1763.980459 [76,] 2178.348696 2072.844252 [77,] 2321.683141 2178.348696 [78,] 2848.023141 2321.683141 [79,] 2898.310919 2848.023141 [80,] 3152.387585 2898.310919 [81,] 3341.049807 3152.387585 [82,] 3281.448696 3341.049807 [83,] 3646.268696 3281.448696 [84,] 3396.960459 3646.268696 [85,] 3666.996015 3396.960459 [86,] 3691.080459 3666.996015 [87,] 1932.750119 3691.080459 [88,] 1763.574563 1932.750119 [89,] 1842.839007 1763.574563 [90,] 2315.579007 1842.839007 [91,] 1674.966785 2315.579007 [92,] 1606.933452 1674.966785 [93,] 607.105674 1606.933452 [94,] 867.284563 607.105674 [95,] 498.484563 867.284563 [96,] 669.026326 498.484563 [97,] 874.901881 669.026326 [98,] 167.506326 874.901881 [99,] 109.940119 167.506326 [100,] 106.854563 109.940119 [101,] -828.160993 106.854563 [102,] -1743.950993 -828.160993 [103,] -2041.233215 -1743.950993 [104,] -2484.896548 -2041.233215 [105,] -2549.504326 -2484.896548 [106,] -2740.305437 -2549.504326 [107,] -2649.695437 -2740.305437 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -544.813985 -688.839541 2 -512.109541 -544.813985 3 -2.695748 -512.109541 4 77.078696 -2.695748 5 179.543141 77.078696 6 48.563141 179.543141 7 -98.819081 48.563141 8 -333.712415 -98.819081 9 -318.880193 -333.712415 10 -339.041304 -318.880193 11 -640.881304 -339.041304 12 -1027.079541 -640.881304 13 -998.533985 -1027.079541 14 -868.139541 -998.533985 15 -862.245748 -868.139541 16 -997.781304 -862.245748 17 -1054.686859 -997.781304 18 -1339.116859 -1054.686859 19 -1203.559081 -1339.116859 20 -921.602415 -1203.559081 21 -928.710193 -921.602415 22 -903.921304 -928.710193 23 -889.621304 -903.921304 24 -1279.949541 -889.621304 25 -1439.983985 -1279.949541 26 -1670.849541 -1439.983985 27 -1764.225748 -1670.849541 28 -1940.691304 -1764.225748 29 -2013.616859 -1940.691304 30 -2189.346859 -2013.616859 31 -1943.579081 -2189.346859 32 -1924.502415 -1943.579081 33 -1848.850193 -1924.502415 34 -2202.401304 -1848.850193 35 -2127.831304 -2202.401304 36 -2293.499541 -2127.831304 37 -2217.133985 -2293.499541 38 -1918.739541 -2217.133985 39 -1585.175748 -1918.739541 40 -1441.711304 -1585.175748 41 -1310.726859 -1441.711304 42 -1112.906859 -1310.726859 43 -1078.889081 -1112.906859 44 -1096.972415 -1078.889081 45 -832.900193 -1096.972415 46 -759.131304 -832.900193 47 -796.411304 -759.131304 48 -1043.549541 -796.411304 49 -1254.123985 -1043.549541 50 -1070.929541 -1254.123985 51 -845.045748 -1070.929541 52 -751.511304 -845.045748 53 -546.386859 -751.511304 54 -280.136859 -546.386859 55 43.160919 -280.136859 56 94.587585 43.160919 57 260.429807 94.587585 58 469.398696 260.429807 59 488.218696 469.398696 60 88.050459 488.218696 61 223.366015 88.050459 62 418.200459 223.366015 63 943.854252 418.200459 64 1005.838696 943.854252 65 1409.513141 1005.838696 66 1453.293141 1409.513141 67 1749.640919 1453.293141 68 1907.777585 1749.640919 69 2270.259807 1907.777585 70 2326.668696 2270.259807 71 2471.468696 2326.668696 72 2178.880459 2471.468696 73 1689.326015 2178.880459 74 1763.980459 1689.326015 75 2072.844252 1763.980459 76 2178.348696 2072.844252 77 2321.683141 2178.348696 78 2848.023141 2321.683141 79 2898.310919 2848.023141 80 3152.387585 2898.310919 81 3341.049807 3152.387585 82 3281.448696 3341.049807 83 3646.268696 3281.448696 84 3396.960459 3646.268696 85 3666.996015 3396.960459 86 3691.080459 3666.996015 87 1932.750119 3691.080459 88 1763.574563 1932.750119 89 1842.839007 1763.574563 90 2315.579007 1842.839007 91 1674.966785 2315.579007 92 1606.933452 1674.966785 93 607.105674 1606.933452 94 867.284563 607.105674 95 498.484563 867.284563 96 669.026326 498.484563 97 874.901881 669.026326 98 167.506326 874.901881 99 109.940119 167.506326 100 106.854563 109.940119 101 -828.160993 106.854563 102 -1743.950993 -828.160993 103 -2041.233215 -1743.950993 104 -2484.896548 -2041.233215 105 -2549.504326 -2484.896548 106 -2740.305437 -2549.504326 107 -2649.695437 -2740.305437 > 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/77no81258647958.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/8xlmn1258647958.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/9z4wg1258647958.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/10s15z1258647958.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/11t2hy1258647958.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/12sioq1258647958.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/13m78t1258647958.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/14bsgg1258647958.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/159xso1258647958.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/16f43j1258647958.tab") + } > > system("convert tmp/1pmv21258647958.ps tmp/1pmv21258647958.png") > system("convert tmp/2f0rb1258647958.ps tmp/2f0rb1258647958.png") > system("convert tmp/30f9m1258647958.ps tmp/30f9m1258647958.png") > system("convert tmp/4kacr1258647958.ps tmp/4kacr1258647958.png") > system("convert tmp/5cxv61258647958.ps tmp/5cxv61258647958.png") > system("convert tmp/6gj2d1258647958.ps tmp/6gj2d1258647958.png") > system("convert tmp/77no81258647958.ps tmp/77no81258647958.png") > system("convert tmp/8xlmn1258647958.ps tmp/8xlmn1258647958.png") > system("convert tmp/9z4wg1258647958.ps tmp/9z4wg1258647958.png") > system("convert tmp/10s15z1258647958.ps tmp/10s15z1258647958.png") > > > proc.time() user system elapsed 3.137 1.636 3.730