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(2050 + ,2650 + ,13 + ,7 + ,1 + ,1 + ,0 + ,1639 + ,2150 + ,2664 + ,6 + ,5 + ,1 + ,1 + ,0 + ,1193 + ,2150 + ,2921 + ,3 + ,6 + ,1 + ,1 + ,0 + ,1635 + ,1999 + ,2580 + ,4 + ,4 + ,1 + ,1 + ,0 + ,1732 + ,1900 + ,2580 + ,4 + ,4 + ,1 + ,0 + ,0 + ,1534 + ,1800 + ,2774 + ,2 + ,4 + ,1 + ,0 + ,0 + ,1765 + ,1560 + ,1920 + ,1 + ,5 + ,1 + ,1 + ,0 + ,1161 + ,1449 + ,1710 + ,1 + ,3 + ,1 + ,1 + ,0 + ,1010 + ,1375 + ,1837 + ,4 + ,5 + ,1 + ,0 + ,0 + ,1191 + ,1270 + ,1880 + ,8 + ,6 + ,1 + ,0 + ,0 + ,930 + ,1250 + ,2150 + ,15 + ,3 + ,1 + ,0 + ,0 + ,984 + ,1235 + ,1894 + ,14 + ,5 + ,1 + ,1 + ,0 + ,1112 + ,1170 + ,1928 + ,18 + ,8 + ,1 + ,1 + ,0 + ,600 + ,1155 + ,1767 + ,16 + ,4 + ,1 + ,0 + ,0 + ,794 + ,1110 + ,1630 + ,15 + ,3 + ,1 + ,0 + ,1 + ,867 + ,1139 + ,1680 + ,17 + ,4 + ,1 + ,0 + ,1 + ,750 + ,995 + ,1500 + ,15 + ,4 + ,1 + ,0 + ,0 + ,743 + ,900 + ,1400 + ,16 + ,2 + ,1 + ,0 + ,1 + ,731 + ,960 + ,1573 + ,17 + ,6 + ,1 + ,0 + ,0 + ,768 + ,1695 + ,2931 + ,28 + ,3 + ,1 + ,0 + ,1 + ,1142 + ,1553 + ,2200 + ,28 + ,4 + ,1 + ,0 + ,0 + ,1035 + ,1020 + ,1478 + ,53 + ,3 + ,1 + ,0 + ,1 + ,626 + ,1020 + ,1713 + ,30 + ,4 + ,1 + ,0 + ,1 + ,600 + ,850 + ,1190 + ,41 + ,1 + ,1 + ,0 + ,0 + ,600 + ,720 + ,1121 + ,46 + ,4 + ,1 + ,0 + ,0 + ,398 + ,749 + ,1733 + ,43 + ,6 + ,1 + ,0 + ,0 + ,656 + ,2150 + ,2848 + ,4 + ,6 + ,1 + ,1 + ,0 + ,1487 + ,1350 + ,2253 + ,23 + ,4 + ,1 + ,1 + ,0 + ,939 + ,1299 + ,2743 + ,25 + ,5 + ,1 + ,1 + ,1 + ,1232 + ,1250 + ,2180 + ,17 + ,4 + ,1 + ,0 + ,1 + ,1141 + ,1239 + ,1706 + ,14 + ,4 + ,1 + ,0 + ,0 + ,810 + ,1125 + ,1710 + ,16 + ,4 + ,1 + ,1 + ,0 + ,800 + ,1080 + ,2200 + ,26 + ,4 + ,1 + ,0 + ,0 + ,1076 + ,1050 + ,1680 + ,13 + ,4 + ,1 + ,0 + ,0 + ,875 + ,1049 + ,1900 + ,34 + ,3 + ,1 + ,0 + ,0 + ,690 + ,934 + ,1543 + ,20 + ,3 + ,1 + ,0 + ,0 + ,820 + ,875 + ,1173 + ,6 + ,4 + ,1 + ,0 + ,0 + ,456 + ,805 + ,1258 + ,7 + ,4 + ,1 + ,0 + ,1 + ,821 + ,759 + ,997 + ,4 + ,4 + ,1 + ,0 + ,0 + ,461 + ,729 + ,1007 + ,19 + ,6 + ,1 + ,0 + ,0 + ,513 + ,710 + ,1083 + ,22 + ,4 + ,1 + ,0 + ,0 + ,504 + ,975 + ,1500 + ,7 + ,3 + ,0 + ,1 + ,1 + ,700 + ,939 + ,1428 + ,40 + ,2 + ,0 + ,0 + ,0 + ,701 + ,2100 + ,2116 + ,25 + ,3 + ,0 + ,1 + ,0 + ,1209 + ,580 + ,1051 + ,15 + ,2 + ,0 + ,0 + ,0 + ,426 + ,1844 + ,2250 + ,40 + ,6 + ,0 + ,1 + ,0 + ,915 + ,699 + ,1400 + ,45 + ,1 + ,0 + ,1 + ,1 + ,481 + ,1160 + ,1720 + ,5 + ,4 + ,0 + ,0 + ,0 + ,867 + ,1109 + ,1740 + ,4 + ,3 + ,0 + ,0 + ,0 + ,816 + ,1129 + ,1700 + ,6 + ,4 + ,0 + ,0 + ,0 + ,725 + ,1050 + ,1620 + ,6 + ,4 + ,0 + ,0 + ,0 + ,800 + ,1045 + ,1630 + ,6 + ,4 + ,0 + ,0 + ,0 + ,750 + ,1050 + ,1920 + ,8 + ,4 + ,0 + ,0 + ,0 + ,944 + ,1020 + ,1606 + ,5 + ,4 + ,0 + ,0 + ,0 + ,811 + ,1000 + ,1535 + ,7 + ,5 + ,0 + ,0 + ,1 + ,668 + ,1030 + ,1540 + ,6 + ,2 + ,0 + ,0 + ,1 + ,826 + ,975 + ,1739 + ,13 + ,3 + ,0 + ,0 + ,0 + ,880 + ,940 + ,1305 + ,5 + ,3 + ,0 + ,0 + ,0 + ,647 + ,920 + ,1415 + ,7 + ,4 + ,0 + ,0 + ,0 + ,866 + ,945 + ,1580 + ,9 + ,3 + ,0 + ,0 + ,0 + ,810 + ,874 + ,1236 + ,3 + ,4 + ,0 + ,0 + ,0 + ,707 + ,872 + ,1229 + ,6 + ,3 + ,0 + ,0 + ,0 + ,721 + ,870 + ,1273 + ,4 + ,4 + ,0 + ,0 + ,0 + ,638 + ,869 + ,1165 + ,7 + ,4 + ,0 + ,0 + ,0 + ,694 + ,766 + ,1200 + ,7 + ,4 + ,0 + ,0 + ,1 + ,634 + ,739 + ,970 + ,4 + ,4 + ,0 + ,0 + ,1 + ,541) + ,dim=c(8 + ,66) + ,dimnames=list(c('PRICE' + ,'SQFT' + ,'AGE' + ,'FEATS' + ,'NE' + ,'CUST' + ,'COR' + ,'TAX') + ,1:66)) > y <- array(NA,dim=c(8,66),dimnames=list(c('PRICE','SQFT','AGE','FEATS','NE','CUST','COR','TAX'),1:66)) > 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' > par3 <- 'No Linear Trend' > par2 <- 'Do not include Seasonal 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, 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 PRICE SQFT AGE FEATS NE CUST COR TAX 1 2050 2650 13 7 1 1 0 1639 2 2150 2664 6 5 1 1 0 1193 3 2150 2921 3 6 1 1 0 1635 4 1999 2580 4 4 1 1 0 1732 5 1900 2580 4 4 1 0 0 1534 6 1800 2774 2 4 1 0 0 1765 7 1560 1920 1 5 1 1 0 1161 8 1449 1710 1 3 1 1 0 1010 9 1375 1837 4 5 1 0 0 1191 10 1270 1880 8 6 1 0 0 930 11 1250 2150 15 3 1 0 0 984 12 1235 1894 14 5 1 1 0 1112 13 1170 1928 18 8 1 1 0 600 14 1155 1767 16 4 1 0 0 794 15 1110 1630 15 3 1 0 1 867 16 1139 1680 17 4 1 0 1 750 17 995 1500 15 4 1 0 0 743 18 900 1400 16 2 1 0 1 731 19 960 1573 17 6 1 0 0 768 20 1695 2931 28 3 1 0 1 1142 21 1553 2200 28 4 1 0 0 1035 22 1020 1478 53 3 1 0 1 626 23 1020 1713 30 4 1 0 1 600 24 850 1190 41 1 1 0 0 600 25 720 1121 46 4 1 0 0 398 26 749 1733 43 6 1 0 0 656 27 2150 2848 4 6 1 1 0 1487 28 1350 2253 23 4 1 1 0 939 29 1299 2743 25 5 1 1 1 1232 30 1250 2180 17 4 1 0 1 1141 31 1239 1706 14 4 1 0 0 810 32 1125 1710 16 4 1 1 0 800 33 1080 2200 26 4 1 0 0 1076 34 1050 1680 13 4 1 0 0 875 35 1049 1900 34 3 1 0 0 690 36 934 1543 20 3 1 0 0 820 37 875 1173 6 4 1 0 0 456 38 805 1258 7 4 1 0 1 821 39 759 997 4 4 1 0 0 461 40 729 1007 19 6 1 0 0 513 41 710 1083 22 4 1 0 0 504 42 975 1500 7 3 0 1 1 700 43 939 1428 40 2 0 0 0 701 44 2100 2116 25 3 0 1 0 1209 45 580 1051 15 2 0 0 0 426 46 1844 2250 40 6 0 1 0 915 47 699 1400 45 1 0 1 1 481 48 1160 1720 5 4 0 0 0 867 49 1109 1740 4 3 0 0 0 816 50 1129 1700 6 4 0 0 0 725 51 1050 1620 6 4 0 0 0 800 52 1045 1630 6 4 0 0 0 750 53 1050 1920 8 4 0 0 0 944 54 1020 1606 5 4 0 0 0 811 55 1000 1535 7 5 0 0 1 668 56 1030 1540 6 2 0 0 1 826 57 975 1739 13 3 0 0 0 880 58 940 1305 5 3 0 0 0 647 59 920 1415 7 4 0 0 0 866 60 945 1580 9 3 0 0 0 810 61 874 1236 3 4 0 0 0 707 62 872 1229 6 3 0 0 0 721 63 870 1273 4 4 0 0 0 638 64 869 1165 7 4 0 0 0 694 65 766 1200 7 4 0 0 1 634 66 739 970 4 4 0 0 1 541 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) SQFT AGE FEATS NE CUST 92.7448 0.3522 -0.5651 4.3896 -17.3853 174.9411 COR TAX -73.5823 0.4989 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -466.28 -82.29 6.75 78.70 484.84 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 92.74480 101.60704 0.913 0.365137 SQFT 0.35222 0.09575 3.679 0.000515 *** AGE -0.56508 2.00253 -0.282 0.778807 FEATS 4.38961 18.55499 0.237 0.813822 NE -17.38534 47.27462 -0.368 0.714397 CUST 174.94108 53.72371 3.256 0.001887 ** COR -73.58234 49.13007 -1.498 0.139633 TAX 0.49887 0.15849 3.148 0.002598 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 158.9 on 58 degrees of freedom Multiple R-squared: 0.8623, Adjusted R-squared: 0.8456 F-statistic: 51.86 on 7 and 58 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.57974342 8.405132e-01 4.202566e-01 [2,] 0.58028806 8.394239e-01 4.197119e-01 [3,] 0.55655986 8.868803e-01 4.434401e-01 [4,] 0.50573667 9.885267e-01 4.942633e-01 [5,] 0.38593969 7.718794e-01 6.140603e-01 [6,] 0.29615377 5.923075e-01 7.038462e-01 [7,] 0.23416121 4.683224e-01 7.658388e-01 [8,] 0.16583580 3.316716e-01 8.341642e-01 [9,] 0.11001897 2.200379e-01 8.899810e-01 [10,] 0.08390088 1.678018e-01 9.160991e-01 [11,] 0.23191535 4.638307e-01 7.680847e-01 [12,] 0.28398606 5.679721e-01 7.160139e-01 [13,] 0.27820611 5.564122e-01 7.217939e-01 [14,] 0.22499602 4.499920e-01 7.750040e-01 [15,] 0.17446007 3.489201e-01 8.255399e-01 [16,] 0.30146848 6.029370e-01 6.985315e-01 [17,] 0.27707770 5.541554e-01 7.229223e-01 [18,] 0.33135544 6.627109e-01 6.686446e-01 [19,] 0.84977285 3.004543e-01 1.502272e-01 [20,] 0.81323218 3.735356e-01 1.867678e-01 [21,] 0.81582271 3.683546e-01 1.841773e-01 [22,] 0.84999152 3.000170e-01 1.500085e-01 [23,] 0.94475789 1.104842e-01 5.524211e-02 [24,] 0.93117508 1.376498e-01 6.882492e-02 [25,] 0.90098142 1.980372e-01 9.901858e-02 [26,] 0.88285411 2.342918e-01 1.171459e-01 [27,] 0.88176121 2.364776e-01 1.182388e-01 [28,] 0.87001329 2.599734e-01 1.299867e-01 [29,] 0.85810293 2.837941e-01 1.418971e-01 [30,] 0.81411930 3.717614e-01 1.858807e-01 [31,] 0.75075537 4.984893e-01 2.492446e-01 [32,] 0.78859373 4.228125e-01 2.114063e-01 [33,] 0.97279824 5.440351e-02 2.720176e-02 [34,] 0.99335527 1.328946e-02 6.644731e-03 [35,] 0.99631950 7.361007e-03 3.680503e-03 [36,] 0.99998571 2.857277e-05 1.428639e-05 [37,] 0.99995799 8.402403e-05 4.201201e-05 [38,] 0.99995424 9.151967e-05 4.575984e-05 [39,] 0.99983514 3.297187e-04 1.648594e-04 [40,] 0.99964560 7.087969e-04 3.543985e-04 [41,] 0.99896882 2.062367e-03 1.031183e-03 [42,] 0.99673837 6.523252e-03 3.261626e-03 [43,] 0.99420576 1.158848e-02 5.794242e-03 [44,] 0.97997950 4.004101e-02 2.002050e-02 [45,] 0.98180810 3.638380e-02 1.819190e-02 > postscript(file="/var/wessaorg/rcomp/tmp/1ww3f1355130893.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/2le431355130893.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/30gsp1355130893.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/44yyx1355130893.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/5o8fk1355130893.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 = 66 Frequency = 1 1 2 3 4 5 6 25.291615 347.680244 30.574711 -39.364930 135.352420 -149.347086 7 8 9 10 11 12 32.869134 79.943580 38.773480 46.703895 -78.209653 -251.182495 13 14 15 16 17 18 -83.644945 52.650763 96.894031 163.391466 11.570355 40.705269 19 20 21 22 23 24 -69.262386 93.814747 184.693512 202.132007 114.944830 74.957406 25 26 27 28 29 30 59.688812 -266.051754 130.684501 -166.849024 -466.282044 -96.775902 31 32 33 34 35 36 149.023998 -135.207091 -309.890324 -63.809960 -33.750716 -95.773014 37 38 39 40 41 42 144.835727 -63.042978 87.201644 27.435230 -3.369076 -105.853354 43 44 45 46 47 48 7.403536 484.844409 -95.747921 323.622350 -207.126661 14.186275 49 50 51 52 53 54 -14.591194 61.635272 -26.602514 -10.181195 -202.975145 -57.724109 55 56 57 58 59 60 88.944711 50.965889 -175.080924 54.497911 -116.758093 -116.417628 61 62 63 64 65 66 -22.650995 -23.084796 -4.695959 6.102145 -5.710951 92.998940 > postscript(file="/var/wessaorg/rcomp/tmp/6w7jm1355130893.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 = 66 Frequency = 1 lag(myerror, k = 1) myerror 0 25.291615 NA 1 347.680244 25.291615 2 30.574711 347.680244 3 -39.364930 30.574711 4 135.352420 -39.364930 5 -149.347086 135.352420 6 32.869134 -149.347086 7 79.943580 32.869134 8 38.773480 79.943580 9 46.703895 38.773480 10 -78.209653 46.703895 11 -251.182495 -78.209653 12 -83.644945 -251.182495 13 52.650763 -83.644945 14 96.894031 52.650763 15 163.391466 96.894031 16 11.570355 163.391466 17 40.705269 11.570355 18 -69.262386 40.705269 19 93.814747 -69.262386 20 184.693512 93.814747 21 202.132007 184.693512 22 114.944830 202.132007 23 74.957406 114.944830 24 59.688812 74.957406 25 -266.051754 59.688812 26 130.684501 -266.051754 27 -166.849024 130.684501 28 -466.282044 -166.849024 29 -96.775902 -466.282044 30 149.023998 -96.775902 31 -135.207091 149.023998 32 -309.890324 -135.207091 33 -63.809960 -309.890324 34 -33.750716 -63.809960 35 -95.773014 -33.750716 36 144.835727 -95.773014 37 -63.042978 144.835727 38 87.201644 -63.042978 39 27.435230 87.201644 40 -3.369076 27.435230 41 -105.853354 -3.369076 42 7.403536 -105.853354 43 484.844409 7.403536 44 -95.747921 484.844409 45 323.622350 -95.747921 46 -207.126661 323.622350 47 14.186275 -207.126661 48 -14.591194 14.186275 49 61.635272 -14.591194 50 -26.602514 61.635272 51 -10.181195 -26.602514 52 -202.975145 -10.181195 53 -57.724109 -202.975145 54 88.944711 -57.724109 55 50.965889 88.944711 56 -175.080924 50.965889 57 54.497911 -175.080924 58 -116.758093 54.497911 59 -116.417628 -116.758093 60 -22.650995 -116.417628 61 -23.084796 -22.650995 62 -4.695959 -23.084796 63 6.102145 -4.695959 64 -5.710951 6.102145 65 92.998940 -5.710951 66 NA 92.998940 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] 347.680244 25.291615 [2,] 30.574711 347.680244 [3,] -39.364930 30.574711 [4,] 135.352420 -39.364930 [5,] -149.347086 135.352420 [6,] 32.869134 -149.347086 [7,] 79.943580 32.869134 [8,] 38.773480 79.943580 [9,] 46.703895 38.773480 [10,] -78.209653 46.703895 [11,] -251.182495 -78.209653 [12,] -83.644945 -251.182495 [13,] 52.650763 -83.644945 [14,] 96.894031 52.650763 [15,] 163.391466 96.894031 [16,] 11.570355 163.391466 [17,] 40.705269 11.570355 [18,] -69.262386 40.705269 [19,] 93.814747 -69.262386 [20,] 184.693512 93.814747 [21,] 202.132007 184.693512 [22,] 114.944830 202.132007 [23,] 74.957406 114.944830 [24,] 59.688812 74.957406 [25,] -266.051754 59.688812 [26,] 130.684501 -266.051754 [27,] -166.849024 130.684501 [28,] -466.282044 -166.849024 [29,] -96.775902 -466.282044 [30,] 149.023998 -96.775902 [31,] -135.207091 149.023998 [32,] -309.890324 -135.207091 [33,] -63.809960 -309.890324 [34,] -33.750716 -63.809960 [35,] -95.773014 -33.750716 [36,] 144.835727 -95.773014 [37,] -63.042978 144.835727 [38,] 87.201644 -63.042978 [39,] 27.435230 87.201644 [40,] -3.369076 27.435230 [41,] -105.853354 -3.369076 [42,] 7.403536 -105.853354 [43,] 484.844409 7.403536 [44,] -95.747921 484.844409 [45,] 323.622350 -95.747921 [46,] -207.126661 323.622350 [47,] 14.186275 -207.126661 [48,] -14.591194 14.186275 [49,] 61.635272 -14.591194 [50,] -26.602514 61.635272 [51,] -10.181195 -26.602514 [52,] -202.975145 -10.181195 [53,] -57.724109 -202.975145 [54,] 88.944711 -57.724109 [55,] 50.965889 88.944711 [56,] -175.080924 50.965889 [57,] 54.497911 -175.080924 [58,] -116.758093 54.497911 [59,] -116.417628 -116.758093 [60,] -22.650995 -116.417628 [61,] -23.084796 -22.650995 [62,] -4.695959 -23.084796 [63,] 6.102145 -4.695959 [64,] -5.710951 6.102145 [65,] 92.998940 -5.710951 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 347.680244 25.291615 2 30.574711 347.680244 3 -39.364930 30.574711 4 135.352420 -39.364930 5 -149.347086 135.352420 6 32.869134 -149.347086 7 79.943580 32.869134 8 38.773480 79.943580 9 46.703895 38.773480 10 -78.209653 46.703895 11 -251.182495 -78.209653 12 -83.644945 -251.182495 13 52.650763 -83.644945 14 96.894031 52.650763 15 163.391466 96.894031 16 11.570355 163.391466 17 40.705269 11.570355 18 -69.262386 40.705269 19 93.814747 -69.262386 20 184.693512 93.814747 21 202.132007 184.693512 22 114.944830 202.132007 23 74.957406 114.944830 24 59.688812 74.957406 25 -266.051754 59.688812 26 130.684501 -266.051754 27 -166.849024 130.684501 28 -466.282044 -166.849024 29 -96.775902 -466.282044 30 149.023998 -96.775902 31 -135.207091 149.023998 32 -309.890324 -135.207091 33 -63.809960 -309.890324 34 -33.750716 -63.809960 35 -95.773014 -33.750716 36 144.835727 -95.773014 37 -63.042978 144.835727 38 87.201644 -63.042978 39 27.435230 87.201644 40 -3.369076 27.435230 41 -105.853354 -3.369076 42 7.403536 -105.853354 43 484.844409 7.403536 44 -95.747921 484.844409 45 323.622350 -95.747921 46 -207.126661 323.622350 47 14.186275 -207.126661 48 -14.591194 14.186275 49 61.635272 -14.591194 50 -26.602514 61.635272 51 -10.181195 -26.602514 52 -202.975145 -10.181195 53 -57.724109 -202.975145 54 88.944711 -57.724109 55 50.965889 88.944711 56 -175.080924 50.965889 57 54.497911 -175.080924 58 -116.758093 54.497911 59 -116.417628 -116.758093 60 -22.650995 -116.417628 61 -23.084796 -22.650995 62 -4.695959 -23.084796 63 6.102145 -4.695959 64 -5.710951 6.102145 65 92.998940 -5.710951 > 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/7g3261355130893.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/82nfk1355130893.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/9cvw31355130893.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/10agss1355130893.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/1181lz1355130893.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/12r3h81355130894.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/13ek5k1355130894.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/147tgd1355130894.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/15fq9c1355130894.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/16750q1355130894.tab") + } > > try(system("convert tmp/1ww3f1355130893.ps tmp/1ww3f1355130893.png",intern=TRUE)) character(0) > try(system("convert tmp/2le431355130893.ps tmp/2le431355130893.png",intern=TRUE)) character(0) > try(system("convert tmp/30gsp1355130893.ps tmp/30gsp1355130893.png",intern=TRUE)) character(0) > try(system("convert tmp/44yyx1355130893.ps tmp/44yyx1355130893.png",intern=TRUE)) character(0) > try(system("convert tmp/5o8fk1355130893.ps tmp/5o8fk1355130893.png",intern=TRUE)) character(0) > try(system("convert tmp/6w7jm1355130893.ps tmp/6w7jm1355130893.png",intern=TRUE)) character(0) > try(system("convert tmp/7g3261355130893.ps tmp/7g3261355130893.png",intern=TRUE)) character(0) > try(system("convert tmp/82nfk1355130893.ps tmp/82nfk1355130893.png",intern=TRUE)) character(0) > try(system("convert tmp/9cvw31355130893.ps tmp/9cvw31355130893.png",intern=TRUE)) character(0) > try(system("convert tmp/10agss1355130893.ps tmp/10agss1355130893.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 5.967 0.877 6.838