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(2440.25 + ,0 + ,2350.44 + ,2408.64 + ,0 + ,2440.25 + ,2472.81 + ,0 + ,2408.64 + ,2407.6 + ,0 + ,2472.81 + ,2454.62 + ,0 + ,2407.6 + ,2448.05 + ,0 + ,2454.62 + ,2497.84 + ,0 + ,2448.05 + ,2645.64 + ,0 + ,2497.84 + ,2756.76 + ,0 + ,2645.64 + ,2849.27 + ,0 + ,2756.76 + ,2921.44 + ,0 + ,2849.27 + ,2981.85 + ,0 + ,2921.44 + ,3080.58 + ,0 + ,2981.85 + ,3106.22 + ,0 + ,3080.58 + ,3119.31 + ,0 + ,3106.22 + ,3061.26 + ,0 + ,3119.31 + ,3097.31 + ,0 + ,3061.26 + ,3161.69 + ,0 + ,3097.31 + ,3257.16 + ,0 + ,3161.69 + ,3277.01 + ,0 + ,3257.16 + ,3295.32 + ,0 + ,3277.01 + ,3363.99 + ,0 + ,3295.32 + ,3494.17 + ,0 + ,3363.99 + ,3667.03 + ,0 + ,3494.17 + ,3813.06 + ,0 + ,3667.03 + ,3917.96 + ,0 + ,3813.06 + ,3895.51 + ,0 + ,3917.96 + ,3801.06 + ,0 + ,3895.51 + ,3570.12 + ,0 + ,3801.06 + ,3701.61 + ,0 + ,3570.12 + ,3862.27 + ,0 + ,3701.61 + ,3970.1 + ,0 + ,3862.27 + ,4138.52 + ,0 + ,3970.1 + ,4199.75 + ,0 + ,4138.52 + ,4290.89 + ,0 + ,4199.75 + ,4443.91 + ,0 + ,4290.89 + ,4502.64 + ,0 + ,4443.91 + ,4356.98 + ,0 + ,4502.64 + ,4591.27 + ,0 + ,4356.98 + ,4696.96 + ,0 + ,4591.27 + ,4621.4 + ,0 + ,4696.96 + ,4562.84 + ,0 + ,4621.4 + ,4202.52 + ,0 + ,4562.84 + ,4296.49 + ,0 + ,4202.52 + ,4435.23 + ,0 + ,4296.49 + ,4105.18 + ,0 + ,4435.23 + ,4116.68 + ,0 + ,4105.18 + ,3844.49 + ,0 + ,4116.68 + ,3720.98 + ,0 + ,3844.49 + ,3674.4 + ,0 + ,3720.98 + ,3857.62 + ,0 + ,3674.4 + ,3801.06 + ,0 + ,3857.62 + ,3504.37 + ,0 + ,3801.06 + ,3032.6 + ,0 + ,3504.37 + ,3047.03 + ,0 + ,3032.6 + ,2962.34 + ,1 + ,3047.03 + ,2197.82 + ,1 + ,2962.34 + ,2014.45 + ,1 + ,2197.82 + ,1862.83 + ,1 + ,2014.45 + ,1905.41 + ,1 + ,1862.83 + ,1810.99 + ,1 + ,1905.41 + ,1670.07 + ,1 + ,1810.99 + ,1864.44 + ,1 + ,1670.07 + ,2052.02 + ,1 + ,1864.44 + ,2029.6 + ,1 + ,2052.02 + ,2070.83 + ,1 + ,2029.6 + ,2293.41 + ,1 + ,2070.83 + ,2443.27 + ,1 + ,2293.41 + ,2513.17 + ,1 + ,2443.27 + ,2466.92 + ,1 + ,2513.17) + ,dim=c(3 + ,70) + ,dimnames=list(c('Y' + ,'X' + ,'Y1') + ,1:70)) > y <- array(NA,dim=c(3,70),dimnames=list(c('Y','X','Y1'),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 = '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.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 Y1 t 1 2440.25 0 2350.44 1 2 2408.64 0 2440.25 2 3 2472.81 0 2408.64 3 4 2407.60 0 2472.81 4 5 2454.62 0 2407.60 5 6 2448.05 0 2454.62 6 7 2497.84 0 2448.05 7 8 2645.64 0 2497.84 8 9 2756.76 0 2645.64 9 10 2849.27 0 2756.76 10 11 2921.44 0 2849.27 11 12 2981.85 0 2921.44 12 13 3080.58 0 2981.85 13 14 3106.22 0 3080.58 14 15 3119.31 0 3106.22 15 16 3061.26 0 3119.31 16 17 3097.31 0 3061.26 17 18 3161.69 0 3097.31 18 19 3257.16 0 3161.69 19 20 3277.01 0 3257.16 20 21 3295.32 0 3277.01 21 22 3363.99 0 3295.32 22 23 3494.17 0 3363.99 23 24 3667.03 0 3494.17 24 25 3813.06 0 3667.03 25 26 3917.96 0 3813.06 26 27 3895.51 0 3917.96 27 28 3801.06 0 3895.51 28 29 3570.12 0 3801.06 29 30 3701.61 0 3570.12 30 31 3862.27 0 3701.61 31 32 3970.10 0 3862.27 32 33 4138.52 0 3970.10 33 34 4199.75 0 4138.52 34 35 4290.89 0 4199.75 35 36 4443.91 0 4290.89 36 37 4502.64 0 4443.91 37 38 4356.98 0 4502.64 38 39 4591.27 0 4356.98 39 40 4696.96 0 4591.27 40 41 4621.40 0 4696.96 41 42 4562.84 0 4621.40 42 43 4202.52 0 4562.84 43 44 4296.49 0 4202.52 44 45 4435.23 0 4296.49 45 46 4105.18 0 4435.23 46 47 4116.68 0 4105.18 47 48 3844.49 0 4116.68 48 49 3720.98 0 3844.49 49 50 3674.40 0 3720.98 50 51 3857.62 0 3674.40 51 52 3801.06 0 3857.62 52 53 3504.37 0 3801.06 53 54 3032.60 0 3504.37 54 55 3047.03 0 3032.60 55 56 2962.34 1 3047.03 56 57 2197.82 1 2962.34 57 58 2014.45 1 2197.82 58 59 1862.83 1 2014.45 59 60 1905.41 1 1862.83 60 61 1810.99 1 1905.41 61 62 1670.07 1 1810.99 62 63 1864.44 1 1670.07 63 64 2052.02 1 1864.44 64 65 2029.60 1 2052.02 65 66 2070.83 1 2029.60 66 67 2293.41 1 2070.83 67 68 2443.27 1 2293.41 68 69 2513.17 1 2443.27 69 70 2466.92 1 2513.17 70 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) X Y1 t 227.0398 -97.2525 0.9455 -0.7861 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -688.066 -65.536 5.529 113.839 275.344 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 227.0399 132.0454 1.719 0.0902 . X -97.2525 139.7743 -0.696 0.4890 Y1 0.9455 0.0480 19.698 <2e-16 *** t -0.7861 2.1460 -0.366 0.7153 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 167.1 on 66 degrees of freedom Multiple R-squared: 0.9628, Adjusted R-squared: 0.9611 F-statistic: 568.7 on 3 and 66 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,] 4.369972e-03 8.739944e-03 0.9956300 [2,] 4.470491e-02 8.940982e-02 0.9552951 [3,] 3.582796e-02 7.165593e-02 0.9641720 [4,] 1.383104e-02 2.766208e-02 0.9861690 [5,] 4.843815e-03 9.687630e-03 0.9951562 [6,] 1.626840e-03 3.253679e-03 0.9983732 [7,] 5.191990e-04 1.038398e-03 0.9994808 [8,] 2.096952e-04 4.193904e-04 0.9997903 [9,] 9.198986e-05 1.839797e-04 0.9999080 [10,] 1.534390e-04 3.068781e-04 0.9998466 [11,] 5.642143e-05 1.128429e-04 0.9999436 [12,] 1.774805e-05 3.549610e-05 0.9999823 [13,] 5.867288e-06 1.173458e-05 0.9999941 [14,] 2.156874e-06 4.313748e-06 0.9999978 [15,] 8.136162e-07 1.627232e-06 0.9999992 [16,] 2.342808e-07 4.685617e-07 0.9999998 [17,] 1.201232e-07 2.402463e-07 0.9999999 [18,] 1.786012e-07 3.572024e-07 0.9999998 [19,] 1.239152e-07 2.478305e-07 0.9999999 [20,] 4.066552e-08 8.133105e-08 1.0000000 [21,] 2.697976e-08 5.395951e-08 1.0000000 [22,] 7.934359e-08 1.586872e-07 0.9999999 [23,] 1.550584e-05 3.101167e-05 0.9999845 [24,] 6.805784e-06 1.361157e-05 0.9999932 [25,] 3.463860e-06 6.927720e-06 0.9999965 [26,] 1.431727e-06 2.863453e-06 0.9999986 [27,] 1.017852e-06 2.035704e-06 0.9999990 [28,] 3.868120e-07 7.736240e-07 0.9999996 [29,] 1.652942e-07 3.305884e-07 0.9999998 [30,] 1.499335e-07 2.998670e-07 0.9999999 [31,] 6.614125e-08 1.322825e-07 0.9999999 [32,] 1.155835e-07 2.311670e-07 0.9999999 [33,] 4.732078e-07 9.464157e-07 0.9999995 [34,] 5.821634e-07 1.164327e-06 0.9999994 [35,] 4.821131e-07 9.642262e-07 0.9999995 [36,] 5.250446e-07 1.050089e-06 0.9999995 [37,] 6.838065e-05 1.367613e-04 0.9999316 [38,] 8.965965e-05 1.793193e-04 0.9999103 [39,] 3.269691e-04 6.539381e-04 0.9996730 [40,] 2.525709e-03 5.051418e-03 0.9974743 [41,] 3.692723e-03 7.385446e-03 0.9963073 [42,] 5.565967e-03 1.113193e-02 0.9944340 [43,] 3.795075e-03 7.590150e-03 0.9962049 [44,] 2.860568e-03 5.721137e-03 0.9971394 [45,] 2.210554e-02 4.421108e-02 0.9778945 [46,] 4.391943e-02 8.783886e-02 0.9560806 [47,] 4.915525e-02 9.831050e-02 0.9508447 [48,] 1.121181e-01 2.242362e-01 0.8878819 [49,] 8.845594e-02 1.769119e-01 0.9115441 [50,] 8.268760e-01 3.462480e-01 0.1731240 [51,] 8.771669e-01 2.456662e-01 0.1228331 [52,] 8.530488e-01 2.939024e-01 0.1469512 [53,] 7.919925e-01 4.160151e-01 0.2080075 [54,] 8.809853e-01 2.380295e-01 0.1190147 [55,] 8.943004e-01 2.113992e-01 0.1056996 [56,] 8.363377e-01 3.273246e-01 0.1636623 [57,] 7.256217e-01 5.487566e-01 0.2743783 > postscript(file="/var/www/html/rcomp/tmp/1qbo71260885394.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/28urf1260885394.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/3ivcj1260885394.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/40x171260885394.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/5wepf1260885394.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 = 70 Frequency = 1 1 2 3 4 5 6 -8.355451 -124.095113 -29.251613 -154.348538 -44.886085 -95.127607 7 8 9 10 11 12 -38.339540 63.169890 35.330419 23.562056 9.049532 2.008570 13 14 15 16 17 18 44.406742 -22.516820 -32.883455 -102.524008 -10.801368 20.279295 19 20 21 22 23 24 55.663814 -13.967403 -13.639567 38.504346 104.542651 155.102970 25 26 27 28 29 30 138.479156 106.093229 -14.754097 -87.191418 -228.042412 122.588510 31 32 33 34 35 36 159.710218 116.421560 183.673906 86.448132 120.480990 188.113807 37 38 39 40 41 42 102.948803 -97.454578 275.343716 160.297558 -14.406716 -0.738291 43 44 45 46 47 48 -304.903443 130.536857 181.213896 -279.229303 45.120574 -237.156626 49 50 51 52 53 54 -102.523642 -31.538274 196.509430 -32.499812 -274.925973 -465.388127 55 56 57 58 59 60 -4.111345 -4.406412 -688.065530 -147.792292 -125.249021 61.474480 61 62 63 64 65 66 -72.419002 -123.278361 205.118241 209.706624 10.714983 73.929296 67 68 69 70 258.312245 198.507945 127.500735 15.946069 > postscript(file="/var/www/html/rcomp/tmp/6i2hr1260885394.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 = 70 Frequency = 1 lag(myerror, k = 1) myerror 0 -8.355451 NA 1 -124.095113 -8.355451 2 -29.251613 -124.095113 3 -154.348538 -29.251613 4 -44.886085 -154.348538 5 -95.127607 -44.886085 6 -38.339540 -95.127607 7 63.169890 -38.339540 8 35.330419 63.169890 9 23.562056 35.330419 10 9.049532 23.562056 11 2.008570 9.049532 12 44.406742 2.008570 13 -22.516820 44.406742 14 -32.883455 -22.516820 15 -102.524008 -32.883455 16 -10.801368 -102.524008 17 20.279295 -10.801368 18 55.663814 20.279295 19 -13.967403 55.663814 20 -13.639567 -13.967403 21 38.504346 -13.639567 22 104.542651 38.504346 23 155.102970 104.542651 24 138.479156 155.102970 25 106.093229 138.479156 26 -14.754097 106.093229 27 -87.191418 -14.754097 28 -228.042412 -87.191418 29 122.588510 -228.042412 30 159.710218 122.588510 31 116.421560 159.710218 32 183.673906 116.421560 33 86.448132 183.673906 34 120.480990 86.448132 35 188.113807 120.480990 36 102.948803 188.113807 37 -97.454578 102.948803 38 275.343716 -97.454578 39 160.297558 275.343716 40 -14.406716 160.297558 41 -0.738291 -14.406716 42 -304.903443 -0.738291 43 130.536857 -304.903443 44 181.213896 130.536857 45 -279.229303 181.213896 46 45.120574 -279.229303 47 -237.156626 45.120574 48 -102.523642 -237.156626 49 -31.538274 -102.523642 50 196.509430 -31.538274 51 -32.499812 196.509430 52 -274.925973 -32.499812 53 -465.388127 -274.925973 54 -4.111345 -465.388127 55 -4.406412 -4.111345 56 -688.065530 -4.406412 57 -147.792292 -688.065530 58 -125.249021 -147.792292 59 61.474480 -125.249021 60 -72.419002 61.474480 61 -123.278361 -72.419002 62 205.118241 -123.278361 63 209.706624 205.118241 64 10.714983 209.706624 65 73.929296 10.714983 66 258.312245 73.929296 67 198.507945 258.312245 68 127.500735 198.507945 69 15.946069 127.500735 70 NA 15.946069 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -124.095113 -8.355451 [2,] -29.251613 -124.095113 [3,] -154.348538 -29.251613 [4,] -44.886085 -154.348538 [5,] -95.127607 -44.886085 [6,] -38.339540 -95.127607 [7,] 63.169890 -38.339540 [8,] 35.330419 63.169890 [9,] 23.562056 35.330419 [10,] 9.049532 23.562056 [11,] 2.008570 9.049532 [12,] 44.406742 2.008570 [13,] -22.516820 44.406742 [14,] -32.883455 -22.516820 [15,] -102.524008 -32.883455 [16,] -10.801368 -102.524008 [17,] 20.279295 -10.801368 [18,] 55.663814 20.279295 [19,] -13.967403 55.663814 [20,] -13.639567 -13.967403 [21,] 38.504346 -13.639567 [22,] 104.542651 38.504346 [23,] 155.102970 104.542651 [24,] 138.479156 155.102970 [25,] 106.093229 138.479156 [26,] -14.754097 106.093229 [27,] -87.191418 -14.754097 [28,] -228.042412 -87.191418 [29,] 122.588510 -228.042412 [30,] 159.710218 122.588510 [31,] 116.421560 159.710218 [32,] 183.673906 116.421560 [33,] 86.448132 183.673906 [34,] 120.480990 86.448132 [35,] 188.113807 120.480990 [36,] 102.948803 188.113807 [37,] -97.454578 102.948803 [38,] 275.343716 -97.454578 [39,] 160.297558 275.343716 [40,] -14.406716 160.297558 [41,] -0.738291 -14.406716 [42,] -304.903443 -0.738291 [43,] 130.536857 -304.903443 [44,] 181.213896 130.536857 [45,] -279.229303 181.213896 [46,] 45.120574 -279.229303 [47,] -237.156626 45.120574 [48,] -102.523642 -237.156626 [49,] -31.538274 -102.523642 [50,] 196.509430 -31.538274 [51,] -32.499812 196.509430 [52,] -274.925973 -32.499812 [53,] -465.388127 -274.925973 [54,] -4.111345 -465.388127 [55,] -4.406412 -4.111345 [56,] -688.065530 -4.406412 [57,] -147.792292 -688.065530 [58,] -125.249021 -147.792292 [59,] 61.474480 -125.249021 [60,] -72.419002 61.474480 [61,] -123.278361 -72.419002 [62,] 205.118241 -123.278361 [63,] 209.706624 205.118241 [64,] 10.714983 209.706624 [65,] 73.929296 10.714983 [66,] 258.312245 73.929296 [67,] 198.507945 258.312245 [68,] 127.500735 198.507945 [69,] 15.946069 127.500735 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -124.095113 -8.355451 2 -29.251613 -124.095113 3 -154.348538 -29.251613 4 -44.886085 -154.348538 5 -95.127607 -44.886085 6 -38.339540 -95.127607 7 63.169890 -38.339540 8 35.330419 63.169890 9 23.562056 35.330419 10 9.049532 23.562056 11 2.008570 9.049532 12 44.406742 2.008570 13 -22.516820 44.406742 14 -32.883455 -22.516820 15 -102.524008 -32.883455 16 -10.801368 -102.524008 17 20.279295 -10.801368 18 55.663814 20.279295 19 -13.967403 55.663814 20 -13.639567 -13.967403 21 38.504346 -13.639567 22 104.542651 38.504346 23 155.102970 104.542651 24 138.479156 155.102970 25 106.093229 138.479156 26 -14.754097 106.093229 27 -87.191418 -14.754097 28 -228.042412 -87.191418 29 122.588510 -228.042412 30 159.710218 122.588510 31 116.421560 159.710218 32 183.673906 116.421560 33 86.448132 183.673906 34 120.480990 86.448132 35 188.113807 120.480990 36 102.948803 188.113807 37 -97.454578 102.948803 38 275.343716 -97.454578 39 160.297558 275.343716 40 -14.406716 160.297558 41 -0.738291 -14.406716 42 -304.903443 -0.738291 43 130.536857 -304.903443 44 181.213896 130.536857 45 -279.229303 181.213896 46 45.120574 -279.229303 47 -237.156626 45.120574 48 -102.523642 -237.156626 49 -31.538274 -102.523642 50 196.509430 -31.538274 51 -32.499812 196.509430 52 -274.925973 -32.499812 53 -465.388127 -274.925973 54 -4.111345 -465.388127 55 -4.406412 -4.111345 56 -688.065530 -4.406412 57 -147.792292 -688.065530 58 -125.249021 -147.792292 59 61.474480 -125.249021 60 -72.419002 61.474480 61 -123.278361 -72.419002 62 205.118241 -123.278361 63 209.706624 205.118241 64 10.714983 209.706624 65 73.929296 10.714983 66 258.312245 73.929296 67 198.507945 258.312245 68 127.500735 198.507945 69 15.946069 127.500735 > 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/727l81260885394.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/88cek1260885394.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/99i031260885394.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/1070cj1260885394.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/11166y1260885394.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/12p5rh1260885394.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/13rwlr1260885394.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/14jt3e1260885394.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/15fqfr1260885394.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/16t3qh1260885394.tab") + } > > try(system("convert tmp/1qbo71260885394.ps tmp/1qbo71260885394.png",intern=TRUE)) character(0) > try(system("convert tmp/28urf1260885394.ps tmp/28urf1260885394.png",intern=TRUE)) character(0) > try(system("convert tmp/3ivcj1260885394.ps tmp/3ivcj1260885394.png",intern=TRUE)) character(0) > try(system("convert tmp/40x171260885394.ps tmp/40x171260885394.png",intern=TRUE)) character(0) > try(system("convert tmp/5wepf1260885394.ps tmp/5wepf1260885394.png",intern=TRUE)) character(0) > try(system("convert tmp/6i2hr1260885394.ps tmp/6i2hr1260885394.png",intern=TRUE)) character(0) > try(system("convert tmp/727l81260885394.ps tmp/727l81260885394.png",intern=TRUE)) character(0) > try(system("convert tmp/88cek1260885394.ps tmp/88cek1260885394.png",intern=TRUE)) character(0) > try(system("convert tmp/99i031260885394.ps tmp/99i031260885394.png",intern=TRUE)) character(0) > try(system("convert tmp/1070cj1260885394.ps tmp/1070cj1260885394.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.532 1.539 3.641