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(7645 + ,27276 + ,111528 + ,27383 + ,7240 + ,27310 + ,111627 + ,27474 + ,7237 + ,27346 + ,111733 + ,27610 + ,7170 + ,27385 + ,111849 + ,27460 + ,7067 + ,27422 + ,111970 + ,27079 + ,7149 + ,27459 + ,112093 + ,26726 + ,6979 + ,27498 + ,112222 + ,25414 + ,6766 + ,27541 + ,112354 + ,25914 + ,6850 + ,27584 + ,112486 + ,26305 + ,6731 + ,27627 + ,112493 + ,25660 + ,6927 + ,27640 + ,112596 + ,26115 + ,7116 + ,27666 + ,112619 + ,26458 + ,11299 + ,27675 + ,112695 + ,26722 + ,10544 + ,27704 + ,112737 + ,26134 + ,10083 + ,27709 + ,112803 + ,26203 + ,9501 + ,27746 + ,112852 + ,26568 + ,9450 + ,27780 + ,112912 + ,26718 + ,8950 + ,27816 + ,113029 + ,26560 + ,8578 + ,27854 + ,113154 + ,25473 + ,8395 + ,27896 + ,113281 + ,25676 + ,7631 + ,27939 + ,113414 + ,25859 + ,7816 + ,27982 + ,113546 + ,25632 + ,7491 + ,28021 + ,113573 + ,25824 + ,7678 + ,28052 + ,113660 + ,26097 + ,15124 + ,28059 + ,113666 + ,26179 + ,15227 + ,28085 + ,113758 + ,25788 + ,15421 + ,28118 + ,113769 + ,26277 + ,15012 + ,28153 + ,113857 + ,26342 + ,14861 + ,28184 + ,113953 + ,26592 + ,14646 + ,28217 + ,114060 + ,26679 + ,14727 + ,28252 + ,114173 + ,25694 + ,14505 + ,28290 + ,114288 + ,26024 + ,13796 + ,28330 + ,114411 + ,26038 + ,13389 + ,28369 + ,114530 + ,25736 + ,12860 + ,28404 + ,114632 + ,25930 + ,12049 + ,28437 + ,114648 + ,26265 + ,14393 + ,28526 + ,114728 + ,26042 + ,15104 + ,28559 + ,114735 + ,24925 + ,14636 + ,28591 + ,114821 + ,25552 + ,14574 + ,28624 + ,114910 + ,26142 + ,14735 + ,28653 + ,115001 + ,26474 + ,14609 + ,28685 + ,115102 + ,26630 + ,14517 + ,28718 + ,115207 + ,25460 + ,14876 + ,28755 + ,115317 + ,25472 + ,15221 + ,28794 + ,115433 + ,25325 + ,15128 + ,28831 + ,115542 + ,25121 + ,15039 + ,28865 + ,115640 + ,25313 + ,14953 + ,28896 + ,115731 + ,25535 + ,13097 + ,28947 + ,115828 + ,25255 + ,13323 + ,28976 + ,115907 + ,24856 + ,13759 + ,29005 + ,115988 + ,25289 + ,13897 + ,29035 + ,116067 + ,25411 + ,13920 + ,29063 + ,116156 + ,25372 + ,13908 + ,29093 + ,116250 + ,25322 + ,14024 + ,29123 + ,116347 + ,24961 + ,13892 + ,29158 + ,116453 + ,24973 + ,13792 + ,29193 + ,116559 + ,25241 + ,13628 + ,29228 + ,116664 + ,24832 + ,13751 + ,29259 + ,116755 + ,24927 + ,13919 + ,29286 + ,116808 + ,25024 + ,12258 + ,29727 + ,116832 + ,25131 + ,12088 + ,29760 + ,116896 + ,24657 + ,12544 + ,29792 + ,116986 + ,24835 + ,12794 + ,29824 + ,117081 + ,25126 + ,12749 + ,29854 + ,117177 + ,25481 + ,12720 + ,29885 + ,117277 + ,25321 + ,12500 + ,29918 + ,117381 + ,24787 + ,12673 + ,29954 + ,117492 + ,24626 + ,12806 + ,29991 + ,117600 + ,24850 + ,12758 + ,30027 + ,117710 + ,24566) + ,dim=c(4 + ,70) + ,dimnames=list(c('unemployment' + ,'black' + ,'males' + ,'highschool ') + ,1:70)) > y <- array(NA,dim=c(4,70),dimnames=list(c('unemployment','black','males','highschool '),1:70)) > 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 unemployment black males highschool\r 1 7645 27276 111528 27383 2 7240 27310 111627 27474 3 7237 27346 111733 27610 4 7170 27385 111849 27460 5 7067 27422 111970 27079 6 7149 27459 112093 26726 7 6979 27498 112222 25414 8 6766 27541 112354 25914 9 6850 27584 112486 26305 10 6731 27627 112493 25660 11 6927 27640 112596 26115 12 7116 27666 112619 26458 13 11299 27675 112695 26722 14 10544 27704 112737 26134 15 10083 27709 112803 26203 16 9501 27746 112852 26568 17 9450 27780 112912 26718 18 8950 27816 113029 26560 19 8578 27854 113154 25473 20 8395 27896 113281 25676 21 7631 27939 113414 25859 22 7816 27982 113546 25632 23 7491 28021 113573 25824 24 7678 28052 113660 26097 25 15124 28059 113666 26179 26 15227 28085 113758 25788 27 15421 28118 113769 26277 28 15012 28153 113857 26342 29 14861 28184 113953 26592 30 14646 28217 114060 26679 31 14727 28252 114173 25694 32 14505 28290 114288 26024 33 13796 28330 114411 26038 34 13389 28369 114530 25736 35 12860 28404 114632 25930 36 12049 28437 114648 26265 37 14393 28526 114728 26042 38 15104 28559 114735 24925 39 14636 28591 114821 25552 40 14574 28624 114910 26142 41 14735 28653 115001 26474 42 14609 28685 115102 26630 43 14517 28718 115207 25460 44 14876 28755 115317 25472 45 15221 28794 115433 25325 46 15128 28831 115542 25121 47 15039 28865 115640 25313 48 14953 28896 115731 25535 49 13097 28947 115828 25255 50 13323 28976 115907 24856 51 13759 29005 115988 25289 52 13897 29035 116067 25411 53 13920 29063 116156 25372 54 13908 29093 116250 25322 55 14024 29123 116347 24961 56 13892 29158 116453 24973 57 13792 29193 116559 25241 58 13628 29228 116664 24832 59 13751 29259 116755 24927 60 13919 29286 116808 25024 61 12258 29727 116832 25131 62 12088 29760 116896 24657 63 12544 29792 116986 24835 64 12794 29824 117081 25126 65 12749 29854 117177 25481 66 12720 29885 117277 25321 67 12500 29918 117381 24787 68 12673 29954 117492 24626 69 12806 29991 117600 24850 70 12758 30027 117710 24566 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) black males `highschool\\r` -4.490e+05 -1.009e+01 6.162e+00 1.642e+00 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -3387.3 -1321.7 -418.8 1327.9 4456.6 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -4.490e+05 6.260e+04 -7.172 8.01e-10 *** black -1.009e+01 2.046e+00 -4.933 5.77e-06 *** males 6.162e+00 9.529e-01 6.466 1.42e-08 *** `highschool\\r` 1.642e+00 5.448e-01 3.013 0.00367 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1896 on 66 degrees of freedom Multiple R-squared: 0.6252, Adjusted R-squared: 0.6082 F-statistic: 36.7 on 3 and 66 DF, p-value: 4.493e-14 > 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,] 3.221244e-04 6.442487e-04 9.996779e-01 [2,] 3.607186e-05 7.214372e-05 9.999639e-01 [3,] 2.445644e-06 4.891287e-06 9.999976e-01 [4,] 3.470107e-07 6.940213e-07 9.999997e-01 [5,] 9.689314e-08 1.937863e-07 9.999999e-01 [6,] 4.079102e-08 8.158205e-08 1.000000e+00 [7,] 4.978421e-02 9.956841e-02 9.502158e-01 [8,] 8.461928e-02 1.692386e-01 9.153807e-01 [9,] 6.406989e-02 1.281398e-01 9.359301e-01 [10,] 4.106998e-02 8.213997e-02 9.589300e-01 [11,] 2.913015e-02 5.826031e-02 9.708698e-01 [12,] 2.695256e-02 5.390512e-02 9.730474e-01 [13,] 1.687856e-02 3.375712e-02 9.831214e-01 [14,] 1.552144e-02 3.104289e-02 9.844786e-01 [15,] 4.756650e-02 9.513301e-02 9.524335e-01 [16,] 1.145525e-01 2.291050e-01 8.854475e-01 [17,] 4.752999e-01 9.505998e-01 5.247001e-01 [18,] 9.993189e-01 1.362282e-03 6.811412e-04 [19,] 9.999994e-01 1.245704e-06 6.228518e-07 [20,] 1.000000e+00 3.410970e-08 1.705485e-08 [21,] 1.000000e+00 1.769866e-08 8.849328e-09 [22,] 1.000000e+00 3.520721e-08 1.760360e-08 [23,] 1.000000e+00 8.692367e-08 4.346184e-08 [24,] 9.999999e-01 2.294003e-07 1.147001e-07 [25,] 9.999998e-01 4.880452e-07 2.440226e-07 [26,] 9.999994e-01 1.260031e-06 6.300156e-07 [27,] 9.999988e-01 2.380729e-06 1.190364e-06 [28,] 9.999989e-01 2.219691e-06 1.109846e-06 [29,] 9.999998e-01 3.716724e-07 1.858362e-07 [30,] 1.000000e+00 5.646232e-12 2.823116e-12 [31,] 1.000000e+00 5.213006e-12 2.606503e-12 [32,] 1.000000e+00 1.641768e-11 8.208838e-12 [33,] 1.000000e+00 5.538983e-11 2.769492e-11 [34,] 1.000000e+00 1.376863e-10 6.884315e-11 [35,] 1.000000e+00 3.730016e-10 1.865008e-10 [36,] 1.000000e+00 7.997489e-10 3.998744e-10 [37,] 1.000000e+00 2.960320e-09 1.480160e-09 [38,] 1.000000e+00 9.494406e-09 4.747203e-09 [39,] 1.000000e+00 8.456122e-09 4.228061e-09 [40,] 1.000000e+00 1.795674e-09 8.978369e-10 [41,] 1.000000e+00 7.116255e-11 3.558127e-11 [42,] 1.000000e+00 3.139012e-14 1.569506e-14 [43,] 1.000000e+00 5.758598e-15 2.879299e-15 [44,] 1.000000e+00 3.114019e-14 1.557009e-14 [45,] 1.000000e+00 2.653077e-13 1.326538e-13 [46,] 1.000000e+00 1.925480e-12 9.627400e-13 [47,] 1.000000e+00 1.334327e-11 6.671634e-12 [48,] 1.000000e+00 9.848228e-11 4.924114e-11 [49,] 1.000000e+00 7.866014e-11 3.933007e-11 [50,] 1.000000e+00 1.727153e-10 8.635766e-11 [51,] 1.000000e+00 2.223358e-09 1.111679e-09 [52,] 1.000000e+00 2.708328e-08 1.354164e-08 [53,] 9.999998e-01 3.040531e-07 1.520265e-07 [54,] 9.999996e-01 7.233939e-07 3.616970e-07 [55,] 9.999995e-01 1.068106e-06 5.340531e-07 [56,] 9.999965e-01 6.921798e-06 3.460899e-06 [57,] 9.999153e-01 1.694371e-04 8.471854e-05 > postscript(file="/var/fisher/rcomp/tmp/1t9pc1352137113.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/fisher/rcomp/tmp/2dvr91352137113.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/fisher/rcomp/tmp/3giz31352137113.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/fisher/rcomp/tmp/4fmnc1352137113.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/fisher/rcomp/tmp/5qeaj1352137113.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 = 70 Frequency = 1 1 2 3 4 5 6 7 -227.8248 -1049.0208 -1565.0344 -1706.8722 -1556.4871 -1279.3915 303.2590 8 9 10 11 12 13 14 -1109.8581 -2047.0359 -716.2570 -1770.6364 -2023.9906 1348.1733 1592.4001 15 16 17 18 19 20 21 661.9249 -447.7137 -771.4601 -1369.6091 -343.7758 -1218.6103 -2668.4878 22 23 24 25 26 27 28 -2490.1296 -2903.0178 -3387.3351 3957.7392 4398.1940 4054.7597 3350.1194 29 30 31 32 33 34 35 2510.1045 1826.0872 3181.1306 2092.3711 1006.2658 755.4728 -367.2035 36 37 38 39 40 41 42 -1493.6332 1621.9003 4456.5933 2752.3898 1506.5362 854.5266 173.1119 43 44 45 46 47 48 49 1687.9648 1722.9625 1988.1997 1931.9543 1267.1139 568.8735 -910.3478 50 51 52 53 54 55 56 -223.3749 -704.5735 -950.8032 -1129.5335 -1335.8260 -922.0522 -1373.5961 57 58 59 60 61 62 63 -2213.4009 -1999.6508 -2280.4024 -2325.6691 141.3341 688.2326 620.4805 64 65 66 67 68 69 70 130.4141 -786.0677 -855.6553 -506.7264 -389.9802 -916.6877 -812.8575 > postscript(file="/var/fisher/rcomp/tmp/659rd1352137113.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 = 70 Frequency = 1 lag(myerror, k = 1) myerror 0 -227.8248 NA 1 -1049.0208 -227.8248 2 -1565.0344 -1049.0208 3 -1706.8722 -1565.0344 4 -1556.4871 -1706.8722 5 -1279.3915 -1556.4871 6 303.2590 -1279.3915 7 -1109.8581 303.2590 8 -2047.0359 -1109.8581 9 -716.2570 -2047.0359 10 -1770.6364 -716.2570 11 -2023.9906 -1770.6364 12 1348.1733 -2023.9906 13 1592.4001 1348.1733 14 661.9249 1592.4001 15 -447.7137 661.9249 16 -771.4601 -447.7137 17 -1369.6091 -771.4601 18 -343.7758 -1369.6091 19 -1218.6103 -343.7758 20 -2668.4878 -1218.6103 21 -2490.1296 -2668.4878 22 -2903.0178 -2490.1296 23 -3387.3351 -2903.0178 24 3957.7392 -3387.3351 25 4398.1940 3957.7392 26 4054.7597 4398.1940 27 3350.1194 4054.7597 28 2510.1045 3350.1194 29 1826.0872 2510.1045 30 3181.1306 1826.0872 31 2092.3711 3181.1306 32 1006.2658 2092.3711 33 755.4728 1006.2658 34 -367.2035 755.4728 35 -1493.6332 -367.2035 36 1621.9003 -1493.6332 37 4456.5933 1621.9003 38 2752.3898 4456.5933 39 1506.5362 2752.3898 40 854.5266 1506.5362 41 173.1119 854.5266 42 1687.9648 173.1119 43 1722.9625 1687.9648 44 1988.1997 1722.9625 45 1931.9543 1988.1997 46 1267.1139 1931.9543 47 568.8735 1267.1139 48 -910.3478 568.8735 49 -223.3749 -910.3478 50 -704.5735 -223.3749 51 -950.8032 -704.5735 52 -1129.5335 -950.8032 53 -1335.8260 -1129.5335 54 -922.0522 -1335.8260 55 -1373.5961 -922.0522 56 -2213.4009 -1373.5961 57 -1999.6508 -2213.4009 58 -2280.4024 -1999.6508 59 -2325.6691 -2280.4024 60 141.3341 -2325.6691 61 688.2326 141.3341 62 620.4805 688.2326 63 130.4141 620.4805 64 -786.0677 130.4141 65 -855.6553 -786.0677 66 -506.7264 -855.6553 67 -389.9802 -506.7264 68 -916.6877 -389.9802 69 -812.8575 -916.6877 70 NA -812.8575 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -1049.0208 -227.8248 [2,] -1565.0344 -1049.0208 [3,] -1706.8722 -1565.0344 [4,] -1556.4871 -1706.8722 [5,] -1279.3915 -1556.4871 [6,] 303.2590 -1279.3915 [7,] -1109.8581 303.2590 [8,] -2047.0359 -1109.8581 [9,] -716.2570 -2047.0359 [10,] -1770.6364 -716.2570 [11,] -2023.9906 -1770.6364 [12,] 1348.1733 -2023.9906 [13,] 1592.4001 1348.1733 [14,] 661.9249 1592.4001 [15,] -447.7137 661.9249 [16,] -771.4601 -447.7137 [17,] -1369.6091 -771.4601 [18,] -343.7758 -1369.6091 [19,] -1218.6103 -343.7758 [20,] -2668.4878 -1218.6103 [21,] -2490.1296 -2668.4878 [22,] -2903.0178 -2490.1296 [23,] -3387.3351 -2903.0178 [24,] 3957.7392 -3387.3351 [25,] 4398.1940 3957.7392 [26,] 4054.7597 4398.1940 [27,] 3350.1194 4054.7597 [28,] 2510.1045 3350.1194 [29,] 1826.0872 2510.1045 [30,] 3181.1306 1826.0872 [31,] 2092.3711 3181.1306 [32,] 1006.2658 2092.3711 [33,] 755.4728 1006.2658 [34,] -367.2035 755.4728 [35,] -1493.6332 -367.2035 [36,] 1621.9003 -1493.6332 [37,] 4456.5933 1621.9003 [38,] 2752.3898 4456.5933 [39,] 1506.5362 2752.3898 [40,] 854.5266 1506.5362 [41,] 173.1119 854.5266 [42,] 1687.9648 173.1119 [43,] 1722.9625 1687.9648 [44,] 1988.1997 1722.9625 [45,] 1931.9543 1988.1997 [46,] 1267.1139 1931.9543 [47,] 568.8735 1267.1139 [48,] -910.3478 568.8735 [49,] -223.3749 -910.3478 [50,] -704.5735 -223.3749 [51,] -950.8032 -704.5735 [52,] -1129.5335 -950.8032 [53,] -1335.8260 -1129.5335 [54,] -922.0522 -1335.8260 [55,] -1373.5961 -922.0522 [56,] -2213.4009 -1373.5961 [57,] -1999.6508 -2213.4009 [58,] -2280.4024 -1999.6508 [59,] -2325.6691 -2280.4024 [60,] 141.3341 -2325.6691 [61,] 688.2326 141.3341 [62,] 620.4805 688.2326 [63,] 130.4141 620.4805 [64,] -786.0677 130.4141 [65,] -855.6553 -786.0677 [66,] -506.7264 -855.6553 [67,] -389.9802 -506.7264 [68,] -916.6877 -389.9802 [69,] -812.8575 -916.6877 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -1049.0208 -227.8248 2 -1565.0344 -1049.0208 3 -1706.8722 -1565.0344 4 -1556.4871 -1706.8722 5 -1279.3915 -1556.4871 6 303.2590 -1279.3915 7 -1109.8581 303.2590 8 -2047.0359 -1109.8581 9 -716.2570 -2047.0359 10 -1770.6364 -716.2570 11 -2023.9906 -1770.6364 12 1348.1733 -2023.9906 13 1592.4001 1348.1733 14 661.9249 1592.4001 15 -447.7137 661.9249 16 -771.4601 -447.7137 17 -1369.6091 -771.4601 18 -343.7758 -1369.6091 19 -1218.6103 -343.7758 20 -2668.4878 -1218.6103 21 -2490.1296 -2668.4878 22 -2903.0178 -2490.1296 23 -3387.3351 -2903.0178 24 3957.7392 -3387.3351 25 4398.1940 3957.7392 26 4054.7597 4398.1940 27 3350.1194 4054.7597 28 2510.1045 3350.1194 29 1826.0872 2510.1045 30 3181.1306 1826.0872 31 2092.3711 3181.1306 32 1006.2658 2092.3711 33 755.4728 1006.2658 34 -367.2035 755.4728 35 -1493.6332 -367.2035 36 1621.9003 -1493.6332 37 4456.5933 1621.9003 38 2752.3898 4456.5933 39 1506.5362 2752.3898 40 854.5266 1506.5362 41 173.1119 854.5266 42 1687.9648 173.1119 43 1722.9625 1687.9648 44 1988.1997 1722.9625 45 1931.9543 1988.1997 46 1267.1139 1931.9543 47 568.8735 1267.1139 48 -910.3478 568.8735 49 -223.3749 -910.3478 50 -704.5735 -223.3749 51 -950.8032 -704.5735 52 -1129.5335 -950.8032 53 -1335.8260 -1129.5335 54 -922.0522 -1335.8260 55 -1373.5961 -922.0522 56 -2213.4009 -1373.5961 57 -1999.6508 -2213.4009 58 -2280.4024 -1999.6508 59 -2325.6691 -2280.4024 60 141.3341 -2325.6691 61 688.2326 141.3341 62 620.4805 688.2326 63 130.4141 620.4805 64 -786.0677 130.4141 65 -855.6553 -786.0677 66 -506.7264 -855.6553 67 -389.9802 -506.7264 68 -916.6877 -389.9802 69 -812.8575 -916.6877 > 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/fisher/rcomp/tmp/7x3j91352137113.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/fisher/rcomp/tmp/8zvdp1352137113.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/fisher/rcomp/tmp/9o8lv1352137113.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/fisher/rcomp/tmp/102vfm1352137113.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/fisher/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/fisher/rcomp/tmp/11htyc1352137113.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/fisher/rcomp/tmp/12x9x31352137113.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/fisher/rcomp/tmp/13i0ze1352137113.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/fisher/rcomp/tmp/14iuxf1352137113.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/fisher/rcomp/tmp/1531tm1352137113.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/fisher/rcomp/tmp/16r61v1352137113.tab") + } > > try(system("convert tmp/1t9pc1352137113.ps tmp/1t9pc1352137113.png",intern=TRUE)) character(0) > try(system("convert tmp/2dvr91352137113.ps tmp/2dvr91352137113.png",intern=TRUE)) character(0) > try(system("convert tmp/3giz31352137113.ps tmp/3giz31352137113.png",intern=TRUE)) character(0) > try(system("convert tmp/4fmnc1352137113.ps tmp/4fmnc1352137113.png",intern=TRUE)) character(0) > try(system("convert tmp/5qeaj1352137113.ps tmp/5qeaj1352137113.png",intern=TRUE)) character(0) > try(system("convert tmp/659rd1352137113.ps tmp/659rd1352137113.png",intern=TRUE)) character(0) > try(system("convert tmp/7x3j91352137113.ps tmp/7x3j91352137113.png",intern=TRUE)) character(0) > try(system("convert tmp/8zvdp1352137113.ps tmp/8zvdp1352137113.png",intern=TRUE)) character(0) > try(system("convert tmp/9o8lv1352137113.ps tmp/9o8lv1352137113.png",intern=TRUE)) character(0) > try(system("convert tmp/102vfm1352137113.ps tmp/102vfm1352137113.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 6.117 1.086 7.201