R version 2.15.2 (2012-10-26) -- "Trick or Treat" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i686-pc-linux-gnu (32-bit) 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(478 + ,184 + ,40 + ,74 + ,11 + ,31 + ,20 + ,494 + ,213 + ,32 + ,72 + ,11 + ,43 + ,18 + ,643 + ,347 + ,57 + ,70 + ,18 + ,16 + ,16 + ,341 + ,565 + ,31 + ,71 + ,11 + ,25 + ,19 + ,773 + ,327 + ,67 + ,72 + ,9 + ,29 + ,24 + ,603 + ,260 + ,25 + ,68 + ,8 + ,32 + ,15 + ,484 + ,325 + ,34 + ,68 + ,12 + ,24 + ,14 + ,546 + ,102 + ,33 + ,62 + ,13 + ,28 + ,11 + ,424 + ,38 + ,36 + ,69 + ,7 + ,25 + ,12 + ,548 + ,226 + ,31 + ,66 + ,9 + ,58 + ,15 + ,506 + ,137 + ,35 + ,60 + ,13 + ,21 + ,9 + ,819 + ,369 + ,30 + ,81 + ,4 + ,77 + ,36 + ,541 + ,109 + ,44 + ,66 + ,9 + ,37 + ,12 + ,491 + ,809 + ,32 + ,67 + ,11 + ,37 + ,16 + ,514 + ,29 + ,30 + ,65 + ,12 + ,35 + ,11 + ,371 + ,245 + ,16 + ,64 + ,10 + ,42 + ,14 + ,457 + ,118 + ,29 + ,64 + ,12 + ,21 + ,10 + ,437 + ,148 + ,36 + ,62 + ,7 + ,81 + ,27 + ,570 + ,387 + ,30 + ,59 + ,15 + ,31 + ,16 + ,432 + ,98 + ,23 + ,56 + ,15 + ,50 + ,15 + ,619 + ,608 + ,33 + ,46 + ,22 + ,24 + ,8 + ,357 + ,218 + ,35 + ,54 + ,14 + ,27 + ,13 + ,623 + ,254 + ,38 + ,54 + ,20 + ,22 + ,11 + ,547 + ,697 + ,44 + ,45 + ,26 + ,18 + ,8 + ,792 + ,827 + ,28 + ,57 + ,12 + ,23 + ,11 + ,799 + ,693 + ,35 + ,57 + ,9 + ,60 + ,18 + ,439 + ,448 + ,31 + ,61 + ,19 + ,14 + ,12 + ,867 + ,942 + ,39 + ,52 + ,17 + ,31 + ,10 + ,912 + ,1017 + ,27 + ,44 + ,21 + ,24 + ,9 + ,462 + ,216 + ,36 + ,43 + ,18 + ,23 + ,8 + ,859 + ,673 + ,38 + ,48 + ,19 + ,22 + ,10 + ,805 + ,989 + ,46 + ,57 + ,14 + ,25 + ,12 + ,652 + ,630 + ,29 + ,47 + ,19 + ,25 + ,9 + ,776 + ,404 + ,32 + ,50 + ,19 + ,21 + ,9 + ,919 + ,692 + ,39 + ,48 + ,16 + ,32 + ,11 + ,732 + ,1517 + ,44 + ,49 + ,13 + ,31 + ,14 + ,657 + ,879 + ,33 + ,72 + ,13 + ,13 + ,22 + ,1419 + ,631 + ,43 + ,59 + ,14 + ,21 + ,13 + ,989 + ,1375 + ,22 + ,49 + ,9 + ,46 + ,13 + ,821 + ,1139 + ,30 + ,54 + ,13 + ,27 + ,12 + ,1740 + ,3545 + ,86 + ,62 + ,22 + ,18 + ,15 + ,815 + ,706 + ,30 + ,47 + ,17 + ,39 + ,11 + ,760 + ,451 + ,32 + ,45 + ,34 + ,15 + ,10 + ,936 + ,433 + ,43 + ,48 + ,26 + ,23 + ,12 + ,863 + ,601 + ,20 + ,69 + ,23 + ,7 + ,12 + ,783 + ,1024 + ,55 + ,42 + ,23 + ,23 + ,11 + ,715 + ,457 + ,44 + ,49 + ,18 + ,30 + ,12 + ,1504 + ,1441 + ,37 + ,57 + ,15 + ,35 + ,13 + ,1324 + ,1022 + ,82 + ,72 + ,22 + ,15 + ,16 + ,940 + ,1244 + ,66 + ,67 + ,26 + ,18 + ,16 + ,478 + ,184 + ,40 + ,74 + ,11 + ,31 + ,20 + ,494 + ,213 + ,32 + ,72 + ,11 + ,43 + ,18 + ,643 + ,347 + ,57 + ,70 + ,18 + ,16 + ,16 + ,341 + ,565 + ,31 + ,71 + ,11 + ,25 + ,19 + ,773 + ,327 + ,67 + ,72 + ,9 + ,29 + ,24 + ,603 + ,260 + ,25 + ,68 + ,8 + ,32 + ,15 + ,484 + ,325 + ,34 + ,68 + ,12 + ,24 + ,14 + ,546 + ,102 + ,33 + ,62 + ,13 + ,28 + ,11 + ,424 + ,38 + ,36 + ,69 + ,7 + ,25 + ,12 + ,548 + ,226 + ,31 + ,66 + ,9 + ,58 + ,15 + ,506 + ,137 + ,35 + ,60 + ,13 + ,21 + ,9 + ,819 + ,369 + ,30 + ,81 + ,4 + ,77 + ,36 + ,541 + ,109 + ,44 + ,66 + ,9 + ,37 + ,12 + ,491 + ,809 + ,32 + ,67 + ,11 + ,37 + ,16 + ,514 + ,29 + ,30 + ,65 + ,12 + ,35 + ,11 + ,371 + ,245 + ,16 + ,64 + ,10 + ,42 + ,14 + ,457 + ,118 + ,29 + ,64 + ,12 + ,21 + ,10 + ,437 + ,148 + ,36 + ,62 + ,7 + ,81 + ,27 + ,570 + ,387 + ,30 + ,59 + ,15 + ,31 + ,16 + ,432 + ,98 + ,23 + ,56 + ,15 + ,50 + ,15 + ,619 + ,608 + ,33 + ,46 + ,22 + ,24 + ,8 + ,357 + ,218 + ,35 + ,54 + ,14 + ,27 + ,13 + ,623 + ,254 + ,38 + ,54 + ,20 + ,22 + ,11 + ,547 + ,697 + ,44 + ,45 + ,26 + ,18 + ,8 + ,792 + ,827 + ,28 + ,57 + ,12 + ,23 + ,11 + ,799 + ,693 + ,35 + ,57 + ,9 + ,60 + ,18 + ,439 + ,448 + ,31 + ,61 + ,19 + ,14 + ,12 + ,867 + ,942 + ,39 + ,52 + ,17 + ,31 + ,10 + ,912 + ,1017 + ,27 + ,44 + ,21 + ,24 + ,9 + ,462 + ,216 + ,36 + ,43 + ,18 + ,23 + ,8 + ,859 + ,673 + ,38 + ,48 + ,19 + ,22 + ,10 + ,805 + ,989 + ,46 + ,57 + ,14 + ,25 + ,12 + ,652 + ,630 + ,29 + ,47 + ,19 + ,25 + ,9 + ,776 + ,404 + ,32 + ,50 + ,19 + ,21 + ,9 + ,919 + ,692 + ,39 + ,48 + ,16 + ,32 + ,11 + ,732 + ,1517 + ,44 + ,49 + ,13 + ,31 + ,14 + ,657 + ,879 + ,33 + ,72 + ,13 + ,13 + ,22 + ,1419 + ,631 + ,43 + ,59 + ,14 + ,21 + ,13 + ,989 + ,1375 + ,22 + ,49 + ,9 + ,46 + ,13 + ,821 + ,1139 + ,30 + ,54 + ,13 + ,27 + ,12 + ,1740 + ,3545 + ,86 + ,62 + ,22 + ,18 + ,15 + ,815 + ,706 + ,30 + ,47 + ,17 + ,39 + ,11 + ,760 + ,451 + ,32 + ,45 + ,34 + ,15 + ,10 + ,936 + ,433 + ,43 + ,48 + ,26 + ,23 + ,12 + ,863 + ,601 + ,20 + ,69 + ,23 + ,7 + ,12 + ,783 + ,1024 + ,55 + ,42 + ,23 + ,23 + ,11 + ,715 + ,457 + ,44 + ,49 + ,18 + ,30 + ,12 + ,1504 + ,1441 + ,37 + ,57 + ,15 + ,35 + ,13 + ,1324 + ,1022 + ,82 + ,72 + ,22 + ,15 + ,16 + ,940 + ,1244 + ,66 + ,67 + ,26 + ,18 + ,16) + ,dim=c(7 + ,100) + ,dimnames=list(c('Yt' + ,'X1' + ,'X2' + ,'X3' + ,'X4' + ,'X5' + ,'X6') + ,1:100)) > y <- array(NA,dim=c(7,100),dimnames=list(c('Yt','X1','X2','X3','X4','X5','X6'),1:100)) > 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 = '1' > library(lattice) > library(lmtest) Loading required package: zoo Attaching package: 'zoo' The following object(s) are masked from 'package:base': as.Date, 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 Yt X1 X2 X3 X4 X5 X6 1 478 184 40 74 11 31 20 2 494 213 32 72 11 43 18 3 643 347 57 70 18 16 16 4 341 565 31 71 11 25 19 5 773 327 67 72 9 29 24 6 603 260 25 68 8 32 15 7 484 325 34 68 12 24 14 8 546 102 33 62 13 28 11 9 424 38 36 69 7 25 12 10 548 226 31 66 9 58 15 11 506 137 35 60 13 21 9 12 819 369 30 81 4 77 36 13 541 109 44 66 9 37 12 14 491 809 32 67 11 37 16 15 514 29 30 65 12 35 11 16 371 245 16 64 10 42 14 17 457 118 29 64 12 21 10 18 437 148 36 62 7 81 27 19 570 387 30 59 15 31 16 20 432 98 23 56 15 50 15 21 619 608 33 46 22 24 8 22 357 218 35 54 14 27 13 23 623 254 38 54 20 22 11 24 547 697 44 45 26 18 8 25 792 827 28 57 12 23 11 26 799 693 35 57 9 60 18 27 439 448 31 61 19 14 12 28 867 942 39 52 17 31 10 29 912 1017 27 44 21 24 9 30 462 216 36 43 18 23 8 31 859 673 38 48 19 22 10 32 805 989 46 57 14 25 12 33 652 630 29 47 19 25 9 34 776 404 32 50 19 21 9 35 919 692 39 48 16 32 11 36 732 1517 44 49 13 31 14 37 657 879 33 72 13 13 22 38 1419 631 43 59 14 21 13 39 989 1375 22 49 9 46 13 40 821 1139 30 54 13 27 12 41 1740 3545 86 62 22 18 15 42 815 706 30 47 17 39 11 43 760 451 32 45 34 15 10 44 936 433 43 48 26 23 12 45 863 601 20 69 23 7 12 46 783 1024 55 42 23 23 11 47 715 457 44 49 18 30 12 48 1504 1441 37 57 15 35 13 49 1324 1022 82 72 22 15 16 50 940 1244 66 67 26 18 16 51 478 184 40 74 11 31 20 52 494 213 32 72 11 43 18 53 643 347 57 70 18 16 16 54 341 565 31 71 11 25 19 55 773 327 67 72 9 29 24 56 603 260 25 68 8 32 15 57 484 325 34 68 12 24 14 58 546 102 33 62 13 28 11 59 424 38 36 69 7 25 12 60 548 226 31 66 9 58 15 61 506 137 35 60 13 21 9 62 819 369 30 81 4 77 36 63 541 109 44 66 9 37 12 64 491 809 32 67 11 37 16 65 514 29 30 65 12 35 11 66 371 245 16 64 10 42 14 67 457 118 29 64 12 21 10 68 437 148 36 62 7 81 27 69 570 387 30 59 15 31 16 70 432 98 23 56 15 50 15 71 619 608 33 46 22 24 8 72 357 218 35 54 14 27 13 73 623 254 38 54 20 22 11 74 547 697 44 45 26 18 8 75 792 827 28 57 12 23 11 76 799 693 35 57 9 60 18 77 439 448 31 61 19 14 12 78 867 942 39 52 17 31 10 79 912 1017 27 44 21 24 9 80 462 216 36 43 18 23 8 81 859 673 38 48 19 22 10 82 805 989 46 57 14 25 12 83 652 630 29 47 19 25 9 84 776 404 32 50 19 21 9 85 919 692 39 48 16 32 11 86 732 1517 44 49 13 31 14 87 657 879 33 72 13 13 22 88 1419 631 43 59 14 21 13 89 989 1375 22 49 9 46 13 90 821 1139 30 54 13 27 12 91 1740 3545 86 62 22 18 15 92 815 706 30 47 17 39 11 93 760 451 32 45 34 15 10 94 936 433 43 48 26 23 12 95 863 601 20 69 23 7 12 96 783 1024 55 42 23 23 11 97 715 457 44 49 18 30 12 98 1504 1441 37 57 15 35 13 99 1324 1022 82 72 22 15 16 100 940 1244 66 67 26 18 16 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) X1 X2 X3 X4 X5 100.3936 0.3323 3.9982 1.8579 7.8389 2.5588 X6 -3.2312 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -291.82 -108.10 -26.78 86.98 705.89 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 100.39361 252.06198 0.398 0.6913 X1 0.33234 0.04054 8.198 1.31e-12 *** X2 3.99817 1.82402 2.192 0.0309 * X3 1.85791 3.56366 0.521 0.6034 X4 7.83886 5.27652 1.486 0.1408 X5 2.55877 2.33024 1.098 0.2750 X6 -3.23116 7.28618 -0.443 0.6585 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 187.7 on 93 degrees of freedom Multiple R-squared: 0.6132, Adjusted R-squared: 0.5882 F-statistic: 24.57 on 6 and 93 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,] 0.214035908 0.428071815 0.7859641 [2,] 0.116195106 0.232390213 0.8838049 [3,] 0.067006308 0.134012615 0.9329937 [4,] 0.029054939 0.058109878 0.9709451 [5,] 0.014154047 0.028308094 0.9858460 [6,] 0.005608341 0.011216681 0.9943917 [7,] 0.004576302 0.009152604 0.9954237 [8,] 0.001880867 0.003761734 0.9981191 [9,] 0.015644581 0.031289162 0.9843554 [10,] 0.012981905 0.025963809 0.9870181 [11,] 0.006731045 0.013462091 0.9932690 [12,] 0.006672548 0.013345096 0.9933275 [13,] 0.006131808 0.012263616 0.9938682 [14,] 0.003919024 0.007838048 0.9960810 [15,] 0.002490549 0.004981098 0.9975095 [16,] 0.008558524 0.017117047 0.9914415 [17,] 0.006496737 0.012993474 0.9935033 [18,] 0.004275124 0.008550248 0.9957249 [19,] 0.003447295 0.006894591 0.9965527 [20,] 0.004591563 0.009183125 0.9954084 [21,] 0.002863034 0.005726068 0.9971370 [22,] 0.003140392 0.006280785 0.9968596 [23,] 0.001895601 0.003791203 0.9981044 [24,] 0.001062665 0.002125330 0.9989373 [25,] 0.001598189 0.003196378 0.9984018 [26,] 0.001589408 0.003178817 0.9984106 [27,] 0.004766246 0.009532491 0.9952338 [28,] 0.002954170 0.005908340 0.9970458 [29,] 0.369391370 0.738782739 0.6306086 [30,] 0.323771383 0.647542765 0.6762286 [31,] 0.268696317 0.537392635 0.7313037 [32,] 0.236204728 0.472409457 0.7637953 [33,] 0.199458770 0.398917540 0.8005412 [34,] 0.177411254 0.354822508 0.8225887 [35,] 0.189430974 0.378861948 0.8105690 [36,] 0.214653977 0.429307953 0.7853460 [37,] 0.199289055 0.398578111 0.8007109 [38,] 0.159130876 0.318261753 0.8408691 [39,] 0.466387960 0.932775921 0.5336120 [40,] 0.518493828 0.963012345 0.4815062 [41,] 0.500000000 1.000000000 0.5000000 [42,] 0.449426204 0.898852408 0.5505738 [43,] 0.398691977 0.797383953 0.6013080 [44,] 0.346708077 0.693416154 0.6532919 [45,] 0.424727783 0.849455565 0.5752722 [46,] 0.376278012 0.752556025 0.6237220 [47,] 0.330685179 0.661370357 0.6693148 [48,] 0.292510261 0.585020522 0.7074897 [49,] 0.242065058 0.484130116 0.7579349 [50,] 0.200189651 0.400379302 0.7998103 [51,] 0.161936233 0.323872465 0.8380638 [52,] 0.127600142 0.255200284 0.8723999 [53,] 0.138257111 0.276514222 0.8617429 [54,] 0.109363481 0.218726961 0.8906365 [55,] 0.140133616 0.280267232 0.8598664 [56,] 0.113309740 0.226619480 0.8866903 [57,] 0.117042161 0.234084323 0.8829578 [58,] 0.119240978 0.238481955 0.8807590 [59,] 0.100671785 0.201343570 0.8993282 [60,] 0.076259610 0.152519220 0.9237404 [61,] 0.079603987 0.159207974 0.9203960 [62,] 0.067087434 0.134174867 0.9329126 [63,] 0.079136719 0.158273437 0.9208633 [64,] 0.060284230 0.120568459 0.9397158 [65,] 0.080333465 0.160666929 0.9196665 [66,] 0.059810069 0.119620138 0.9401899 [67,] 0.054213423 0.108426845 0.9457866 [68,] 0.090258019 0.180516037 0.9097420 [69,] 0.078751239 0.157502478 0.9212488 [70,] 0.061297674 0.122595348 0.9387023 [71,] 0.057701898 0.115403796 0.9422981 [72,] 0.040506733 0.081013465 0.9594933 [73,] 0.045405777 0.090811553 0.9545942 [74,] 0.042450731 0.084901463 0.9575493 [75,] 0.033841086 0.067682172 0.9661589 [76,] 0.021466826 0.042933651 0.9785332 [77,] 0.028014019 0.056028038 0.9719860 [78,] 0.014962531 0.029925061 0.9850375 [79,] 0.255340278 0.510680557 0.7446597 [80,] 0.160532498 0.321064995 0.8394675 [81,] 0.091683720 0.183367440 0.9083163 > postscript(file="/var/wessaorg/rcomp/tmp/1dapa1353455226.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index') > points(x[,1]-mysum$resid) > grid() > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/2td581353455226.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index') > grid() > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/3k1411353455226.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals') > grid() > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/4kyp01353455226.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals') > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/5u0i91353455226.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > qqnorm(mysum$resid, main='Residual Normal Q-Q Plot') > qqline(mysum$resid) > grid() > dev.off() null device 1 > (myerror <- as.ts(mysum$resid)) Time Series: Start = 1 End = 100 Frequency = 1 1 2 3 4 5 6 -81.882079 -76.986177 -61.005364 -291.823520 95.078842 93.782435 7 8 9 10 11 12 -96.919451 26.569808 -41.219933 -44.558207 -17.893430 213.479251 13 14 15 16 17 18 -20.610430 -259.878862 15.178630 -134.313880 -32.951621 -146.595211 19 20 21 22 23 24 -33.776048 -94.017631 -108.874004 -192.932039 8.407685 -268.576777 25 26 27 28 29 30 81.539001 56.545123 -193.543584 20.734336 86.975268 -108.104677 31 32 33 34 35 36 120.913875 -48.830736 -44.861649 146.913213 171.761359 -275.496179 37 38 39 40 41 42 -65.310438 705.890659 106.398167 -10.415366 -167.605302 75.199880 43 44 45 46 47 48 25.583830 206.715290 195.281491 -159.240760 16.682628 505.741350 49 50 51 52 53 54 263.200669 -160.349435 -81.882079 -76.986177 -61.005364 -291.823520 55 56 57 58 59 60 95.078842 93.782435 -96.919451 26.569808 -41.219933 -44.558207 61 62 63 64 65 66 -17.893430 213.479251 -20.610430 -259.878862 15.178630 -134.313880 67 68 69 70 71 72 -32.951621 -146.595211 -33.776048 -94.017631 -108.874004 -192.932039 73 74 75 76 77 78 8.407685 -268.576777 81.539001 56.545123 -193.543584 20.734336 79 80 81 82 83 84 86.975268 -108.104677 120.913875 -48.830736 -44.861649 146.913213 85 86 87 88 89 90 171.761359 -275.496179 -65.310438 705.890659 106.398167 -10.415366 91 92 93 94 95 96 -167.605302 75.199880 25.583830 206.715290 195.281491 -159.240760 97 98 99 100 16.682628 505.741350 263.200669 -160.349435 > postscript(file="/var/wessaorg/rcomp/tmp/6natt1353455226.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > dum <- cbind(lag(myerror,k=1),myerror) > dum Time Series: Start = 0 End = 100 Frequency = 1 lag(myerror, k = 1) myerror 0 -81.882079 NA 1 -76.986177 -81.882079 2 -61.005364 -76.986177 3 -291.823520 -61.005364 4 95.078842 -291.823520 5 93.782435 95.078842 6 -96.919451 93.782435 7 26.569808 -96.919451 8 -41.219933 26.569808 9 -44.558207 -41.219933 10 -17.893430 -44.558207 11 213.479251 -17.893430 12 -20.610430 213.479251 13 -259.878862 -20.610430 14 15.178630 -259.878862 15 -134.313880 15.178630 16 -32.951621 -134.313880 17 -146.595211 -32.951621 18 -33.776048 -146.595211 19 -94.017631 -33.776048 20 -108.874004 -94.017631 21 -192.932039 -108.874004 22 8.407685 -192.932039 23 -268.576777 8.407685 24 81.539001 -268.576777 25 56.545123 81.539001 26 -193.543584 56.545123 27 20.734336 -193.543584 28 86.975268 20.734336 29 -108.104677 86.975268 30 120.913875 -108.104677 31 -48.830736 120.913875 32 -44.861649 -48.830736 33 146.913213 -44.861649 34 171.761359 146.913213 35 -275.496179 171.761359 36 -65.310438 -275.496179 37 705.890659 -65.310438 38 106.398167 705.890659 39 -10.415366 106.398167 40 -167.605302 -10.415366 41 75.199880 -167.605302 42 25.583830 75.199880 43 206.715290 25.583830 44 195.281491 206.715290 45 -159.240760 195.281491 46 16.682628 -159.240760 47 505.741350 16.682628 48 263.200669 505.741350 49 -160.349435 263.200669 50 -81.882079 -160.349435 51 -76.986177 -81.882079 52 -61.005364 -76.986177 53 -291.823520 -61.005364 54 95.078842 -291.823520 55 93.782435 95.078842 56 -96.919451 93.782435 57 26.569808 -96.919451 58 -41.219933 26.569808 59 -44.558207 -41.219933 60 -17.893430 -44.558207 61 213.479251 -17.893430 62 -20.610430 213.479251 63 -259.878862 -20.610430 64 15.178630 -259.878862 65 -134.313880 15.178630 66 -32.951621 -134.313880 67 -146.595211 -32.951621 68 -33.776048 -146.595211 69 -94.017631 -33.776048 70 -108.874004 -94.017631 71 -192.932039 -108.874004 72 8.407685 -192.932039 73 -268.576777 8.407685 74 81.539001 -268.576777 75 56.545123 81.539001 76 -193.543584 56.545123 77 20.734336 -193.543584 78 86.975268 20.734336 79 -108.104677 86.975268 80 120.913875 -108.104677 81 -48.830736 120.913875 82 -44.861649 -48.830736 83 146.913213 -44.861649 84 171.761359 146.913213 85 -275.496179 171.761359 86 -65.310438 -275.496179 87 705.890659 -65.310438 88 106.398167 705.890659 89 -10.415366 106.398167 90 -167.605302 -10.415366 91 75.199880 -167.605302 92 25.583830 75.199880 93 206.715290 25.583830 94 195.281491 206.715290 95 -159.240760 195.281491 96 16.682628 -159.240760 97 505.741350 16.682628 98 263.200669 505.741350 99 -160.349435 263.200669 100 NA -160.349435 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -76.986177 -81.882079 [2,] -61.005364 -76.986177 [3,] -291.823520 -61.005364 [4,] 95.078842 -291.823520 [5,] 93.782435 95.078842 [6,] -96.919451 93.782435 [7,] 26.569808 -96.919451 [8,] -41.219933 26.569808 [9,] -44.558207 -41.219933 [10,] -17.893430 -44.558207 [11,] 213.479251 -17.893430 [12,] -20.610430 213.479251 [13,] -259.878862 -20.610430 [14,] 15.178630 -259.878862 [15,] -134.313880 15.178630 [16,] -32.951621 -134.313880 [17,] -146.595211 -32.951621 [18,] -33.776048 -146.595211 [19,] -94.017631 -33.776048 [20,] -108.874004 -94.017631 [21,] -192.932039 -108.874004 [22,] 8.407685 -192.932039 [23,] -268.576777 8.407685 [24,] 81.539001 -268.576777 [25,] 56.545123 81.539001 [26,] -193.543584 56.545123 [27,] 20.734336 -193.543584 [28,] 86.975268 20.734336 [29,] -108.104677 86.975268 [30,] 120.913875 -108.104677 [31,] -48.830736 120.913875 [32,] -44.861649 -48.830736 [33,] 146.913213 -44.861649 [34,] 171.761359 146.913213 [35,] -275.496179 171.761359 [36,] -65.310438 -275.496179 [37,] 705.890659 -65.310438 [38,] 106.398167 705.890659 [39,] -10.415366 106.398167 [40,] -167.605302 -10.415366 [41,] 75.199880 -167.605302 [42,] 25.583830 75.199880 [43,] 206.715290 25.583830 [44,] 195.281491 206.715290 [45,] -159.240760 195.281491 [46,] 16.682628 -159.240760 [47,] 505.741350 16.682628 [48,] 263.200669 505.741350 [49,] -160.349435 263.200669 [50,] -81.882079 -160.349435 [51,] -76.986177 -81.882079 [52,] -61.005364 -76.986177 [53,] -291.823520 -61.005364 [54,] 95.078842 -291.823520 [55,] 93.782435 95.078842 [56,] -96.919451 93.782435 [57,] 26.569808 -96.919451 [58,] -41.219933 26.569808 [59,] -44.558207 -41.219933 [60,] -17.893430 -44.558207 [61,] 213.479251 -17.893430 [62,] -20.610430 213.479251 [63,] -259.878862 -20.610430 [64,] 15.178630 -259.878862 [65,] -134.313880 15.178630 [66,] -32.951621 -134.313880 [67,] -146.595211 -32.951621 [68,] -33.776048 -146.595211 [69,] -94.017631 -33.776048 [70,] -108.874004 -94.017631 [71,] -192.932039 -108.874004 [72,] 8.407685 -192.932039 [73,] -268.576777 8.407685 [74,] 81.539001 -268.576777 [75,] 56.545123 81.539001 [76,] -193.543584 56.545123 [77,] 20.734336 -193.543584 [78,] 86.975268 20.734336 [79,] -108.104677 86.975268 [80,] 120.913875 -108.104677 [81,] -48.830736 120.913875 [82,] -44.861649 -48.830736 [83,] 146.913213 -44.861649 [84,] 171.761359 146.913213 [85,] -275.496179 171.761359 [86,] -65.310438 -275.496179 [87,] 705.890659 -65.310438 [88,] 106.398167 705.890659 [89,] -10.415366 106.398167 [90,] -167.605302 -10.415366 [91,] 75.199880 -167.605302 [92,] 25.583830 75.199880 [93,] 206.715290 25.583830 [94,] 195.281491 206.715290 [95,] -159.240760 195.281491 [96,] 16.682628 -159.240760 [97,] 505.741350 16.682628 [98,] 263.200669 505.741350 [99,] -160.349435 263.200669 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -76.986177 -81.882079 2 -61.005364 -76.986177 3 -291.823520 -61.005364 4 95.078842 -291.823520 5 93.782435 95.078842 6 -96.919451 93.782435 7 26.569808 -96.919451 8 -41.219933 26.569808 9 -44.558207 -41.219933 10 -17.893430 -44.558207 11 213.479251 -17.893430 12 -20.610430 213.479251 13 -259.878862 -20.610430 14 15.178630 -259.878862 15 -134.313880 15.178630 16 -32.951621 -134.313880 17 -146.595211 -32.951621 18 -33.776048 -146.595211 19 -94.017631 -33.776048 20 -108.874004 -94.017631 21 -192.932039 -108.874004 22 8.407685 -192.932039 23 -268.576777 8.407685 24 81.539001 -268.576777 25 56.545123 81.539001 26 -193.543584 56.545123 27 20.734336 -193.543584 28 86.975268 20.734336 29 -108.104677 86.975268 30 120.913875 -108.104677 31 -48.830736 120.913875 32 -44.861649 -48.830736 33 146.913213 -44.861649 34 171.761359 146.913213 35 -275.496179 171.761359 36 -65.310438 -275.496179 37 705.890659 -65.310438 38 106.398167 705.890659 39 -10.415366 106.398167 40 -167.605302 -10.415366 41 75.199880 -167.605302 42 25.583830 75.199880 43 206.715290 25.583830 44 195.281491 206.715290 45 -159.240760 195.281491 46 16.682628 -159.240760 47 505.741350 16.682628 48 263.200669 505.741350 49 -160.349435 263.200669 50 -81.882079 -160.349435 51 -76.986177 -81.882079 52 -61.005364 -76.986177 53 -291.823520 -61.005364 54 95.078842 -291.823520 55 93.782435 95.078842 56 -96.919451 93.782435 57 26.569808 -96.919451 58 -41.219933 26.569808 59 -44.558207 -41.219933 60 -17.893430 -44.558207 61 213.479251 -17.893430 62 -20.610430 213.479251 63 -259.878862 -20.610430 64 15.178630 -259.878862 65 -134.313880 15.178630 66 -32.951621 -134.313880 67 -146.595211 -32.951621 68 -33.776048 -146.595211 69 -94.017631 -33.776048 70 -108.874004 -94.017631 71 -192.932039 -108.874004 72 8.407685 -192.932039 73 -268.576777 8.407685 74 81.539001 -268.576777 75 56.545123 81.539001 76 -193.543584 56.545123 77 20.734336 -193.543584 78 86.975268 20.734336 79 -108.104677 86.975268 80 120.913875 -108.104677 81 -48.830736 120.913875 82 -44.861649 -48.830736 83 146.913213 -44.861649 84 171.761359 146.913213 85 -275.496179 171.761359 86 -65.310438 -275.496179 87 705.890659 -65.310438 88 106.398167 705.890659 89 -10.415366 106.398167 90 -167.605302 -10.415366 91 75.199880 -167.605302 92 25.583830 75.199880 93 206.715290 25.583830 94 195.281491 206.715290 95 -159.240760 195.281491 96 16.682628 -159.240760 97 505.741350 16.682628 98 263.200669 505.741350 99 -160.349435 263.200669 > 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/wessaorg/rcomp/tmp/7uveg1353455226.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function') > grid() > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/83xn41353455226.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function') > grid() > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/9gjup1353455226.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0)) > plot(mylm, las = 1, sub='Residual Diagnostics') > par(opar) > dev.off() null device 1 > if (n > n25) { + postscript(file="/var/wessaorg/rcomp/tmp/10xgvj1353455226.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) + plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint') + grid() + dev.off() + } null device 1 > > #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11a2tc1353455226.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/wessaorg/rcomp/tmp/1250is1353455226.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/wessaorg/rcomp/tmp/13q2x11353455226.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/wessaorg/rcomp/tmp/14h8go1353455226.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/wessaorg/rcomp/tmp/153e5u1353455226.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/wessaorg/rcomp/tmp/164eyv1353455226.tab") + } > > try(system("convert tmp/1dapa1353455226.ps tmp/1dapa1353455226.png",intern=TRUE)) character(0) > try(system("convert tmp/2td581353455226.ps tmp/2td581353455226.png",intern=TRUE)) character(0) > try(system("convert tmp/3k1411353455226.ps tmp/3k1411353455226.png",intern=TRUE)) character(0) > try(system("convert tmp/4kyp01353455226.ps tmp/4kyp01353455226.png",intern=TRUE)) character(0) > try(system("convert tmp/5u0i91353455226.ps tmp/5u0i91353455226.png",intern=TRUE)) character(0) > try(system("convert tmp/6natt1353455226.ps tmp/6natt1353455226.png",intern=TRUE)) character(0) > try(system("convert tmp/7uveg1353455226.ps tmp/7uveg1353455226.png",intern=TRUE)) character(0) > try(system("convert tmp/83xn41353455226.ps tmp/83xn41353455226.png",intern=TRUE)) character(0) > try(system("convert tmp/9gjup1353455226.ps tmp/9gjup1353455226.png",intern=TRUE)) character(0) > try(system("convert tmp/10xgvj1353455226.ps tmp/10xgvj1353455226.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 10.529 1.803 12.394