R version 2.8.0 (2008-10-20) Copyright (C) 2008 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. Natural language support but running in an English locale 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(2649.2 + ,31077 + ,0 + ,2579.4 + ,31293 + ,0 + ,2504.6 + ,30236 + ,0 + ,2462.3 + ,30160 + ,0 + ,2467.4 + ,32436 + ,0 + ,2446.7 + ,30695 + ,0 + ,2656.3 + ,27525 + ,0 + ,2626.2 + ,26434 + ,0 + ,2482.6 + ,25739 + ,0 + ,2539.9 + ,25204 + ,0 + ,2502.7 + ,24977 + ,0 + ,2466.9 + ,24320 + ,0 + ,2513.2 + ,22680 + ,1 + ,2443.3 + ,22052 + ,1 + ,2293.4 + ,21467 + ,1 + ,2070.8 + ,21383 + ,1 + ,2029.6 + ,21777 + ,1 + ,2052 + ,21928 + ,1 + ,1864.4 + ,21814 + ,1 + ,1670.1 + ,22937 + ,1 + ,1811 + ,23595 + ,1 + ,1905.4 + ,20830 + ,1 + ,1862.8 + ,19650 + ,1 + ,2014.5 + ,19195 + ,1 + ,2197.8 + ,19644 + ,1 + ,2962.3 + ,18483 + ,0 + ,3047 + ,18079 + ,0 + ,3032.6 + ,19178 + ,0 + ,3504.4 + ,18391 + ,0 + ,3801.1 + ,18441 + ,0 + ,3857.6 + ,18584 + ,0 + ,3674.4 + ,20108 + ,0 + ,3721 + ,20148 + ,0 + ,3844.5 + ,19394 + ,0 + ,4116.7 + ,17745 + ,0 + ,4105.2 + ,17696 + ,0 + ,4435.2 + ,17032 + ,0 + ,4296.5 + ,16438 + ,0 + ,4202.5 + ,15683 + ,0 + ,4562.8 + ,15594 + ,0 + ,4621.4 + ,15713 + ,0 + ,4697 + ,15937 + ,0 + ,4591.3 + ,16171 + ,0 + ,4357 + ,15928 + ,0 + ,4502.6 + ,16348 + ,0 + ,4443.9 + ,15579 + ,0 + ,4290.9 + ,15305 + ,0 + ,4199.8 + ,15648 + ,0 + ,4138.5 + ,14954 + ,0 + ,3970.1 + ,15137 + ,0 + ,3862.3 + ,15839 + ,0 + ,3701.6 + ,16050 + ,0 + ,3570.12 + ,15168 + ,0 + ,3801.06 + ,17064 + ,0 + ,3895.51 + ,16005 + ,0 + ,3917.96 + ,14886 + ,0 + ,3813.06 + ,14931 + ,0 + ,3667.03 + ,14544 + ,0 + ,3494.17 + ,13812 + ,0 + ,3364 + ,13031 + ,0 + ,3295.3 + ,12574 + ,0 + ,3277.0 + ,11964 + ,0 + ,3257.2 + ,11451 + ,0 + ,3161.7 + ,11346 + ,0 + ,3097.3 + ,11353 + ,0 + ,3061.3 + ,10702 + ,0 + ,3119.3 + ,10646 + ,0 + ,3106.22 + ,10556 + ,0 + ,3080.58 + ,10463 + ,0 + ,2981.85 + ,10407 + ,0) + ,dim=c(3 + ,70) + ,dimnames=list(c('Bel20' + ,'Goudprijs' + ,'Crisis') + ,1:70)) > y <- array(NA,dim=c(3,70),dimnames=list(c('Bel20','Goudprijs','Crisis'),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' > 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 Bel20 Goudprijs Crisis t 1 2649.20 31077 0 1 2 2579.40 31293 0 2 3 2504.60 30236 0 3 4 2462.30 30160 0 4 5 2467.40 32436 0 5 6 2446.70 30695 0 6 7 2656.30 27525 0 7 8 2626.20 26434 0 8 9 2482.60 25739 0 9 10 2539.90 25204 0 10 11 2502.70 24977 0 11 12 2466.90 24320 0 12 13 2513.20 22680 1 13 14 2443.30 22052 1 14 15 2293.40 21467 1 15 16 2070.80 21383 1 16 17 2029.60 21777 1 17 18 2052.00 21928 1 18 19 1864.40 21814 1 19 20 1670.10 22937 1 20 21 1811.00 23595 1 21 22 1905.40 20830 1 22 23 1862.80 19650 1 23 24 2014.50 19195 1 24 25 2197.80 19644 1 25 26 2962.30 18483 0 26 27 3047.00 18079 0 27 28 3032.60 19178 0 28 29 3504.40 18391 0 29 30 3801.10 18441 0 30 31 3857.60 18584 0 31 32 3674.40 20108 0 32 33 3721.00 20148 0 33 34 3844.50 19394 0 34 35 4116.70 17745 0 35 36 4105.20 17696 0 36 37 4435.20 17032 0 37 38 4296.50 16438 0 38 39 4202.50 15683 0 39 40 4562.80 15594 0 40 41 4621.40 15713 0 41 42 4697.00 15937 0 42 43 4591.30 16171 0 43 44 4357.00 15928 0 44 45 4502.60 16348 0 45 46 4443.90 15579 0 46 47 4290.90 15305 0 47 48 4199.80 15648 0 48 49 4138.50 14954 0 49 50 3970.10 15137 0 50 51 3862.30 15839 0 51 52 3701.60 16050 0 52 53 3570.12 15168 0 53 54 3801.06 17064 0 54 55 3895.51 16005 0 55 56 3917.96 14886 0 56 57 3813.06 14931 0 57 58 3667.03 14544 0 58 59 3494.17 13812 0 59 60 3364.00 13031 0 60 61 3295.30 12574 0 61 62 3277.00 11964 0 62 63 3257.20 11451 0 63 64 3161.70 11346 0 64 65 3097.30 11353 0 65 66 3061.30 10702 0 66 67 3119.30 10646 0 67 68 3106.22 10556 0 68 69 3080.58 10463 0 69 70 2981.85 10407 0 70 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) Goudprijs Crisis t 7157.5498 -0.1459 -1480.7553 -25.7971 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -851.55 -429.04 38.88 389.39 948.08 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.158e+03 1.150e+03 6.224 3.79e-08 *** Goudprijs -1.459e-01 3.933e-02 -3.710 0.000428 *** Crisis -1.481e+03 1.970e+02 -7.517 1.93e-10 *** t -2.580e+01 1.144e+01 -2.255 0.027484 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 527 on 66 degrees of freedom Multiple R-squared: 0.6312, Adjusted R-squared: 0.6145 F-statistic: 37.66 on 3 and 66 DF, p-value: 2.654e-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,] 5.800093e-03 1.160019e-02 9.941999e-01 [2,] 7.853531e-04 1.570706e-03 9.992146e-01 [3,] 2.828283e-04 5.656565e-04 9.997172e-01 [4,] 4.781839e-05 9.563678e-05 9.999522e-01 [5,] 9.957801e-06 1.991560e-05 9.999900e-01 [6,] 3.963871e-06 7.927742e-06 9.999960e-01 [7,] 5.859381e-07 1.171876e-06 9.999994e-01 [8,] 1.035375e-07 2.070749e-07 9.999999e-01 [9,] 1.092538e-07 2.185077e-07 9.999999e-01 [10,] 1.231271e-06 2.462541e-06 9.999988e-01 [11,] 8.130403e-07 1.626081e-06 9.999992e-01 [12,] 1.873481e-07 3.746962e-07 9.999998e-01 [13,] 8.973008e-08 1.794602e-07 9.999999e-01 [14,] 3.809700e-08 7.619399e-08 1.000000e+00 [15,] 2.087936e-08 4.175872e-08 1.000000e+00 [16,] 5.432558e-09 1.086512e-08 1.000000e+00 [17,] 1.197393e-09 2.394786e-09 1.000000e+00 [18,] 7.554235e-10 1.510847e-09 1.000000e+00 [19,] 1.677307e-08 3.354613e-08 1.000000e+00 [20,] 3.571394e-06 7.142789e-06 9.999964e-01 [21,] 2.897564e-05 5.795128e-05 9.999710e-01 [22,] 3.471383e-04 6.942766e-04 9.996529e-01 [23,] 8.920675e-03 1.784135e-02 9.910793e-01 [24,] 8.828805e-02 1.765761e-01 9.117120e-01 [25,] 2.490482e-01 4.980964e-01 7.509518e-01 [26,] 4.434805e-01 8.869611e-01 5.565195e-01 [27,] 7.131806e-01 5.736388e-01 2.868194e-01 [28,] 9.339399e-01 1.321202e-01 6.606009e-02 [29,] 9.849822e-01 3.003559e-02 1.501779e-02 [30,] 9.986141e-01 2.771773e-03 1.385887e-03 [31,] 9.994746e-01 1.050815e-03 5.254074e-04 [32,] 9.998477e-01 3.045430e-04 1.522715e-04 [33,] 9.999925e-01 1.500897e-05 7.504485e-06 [34,] 9.999924e-01 1.514050e-05 7.570251e-06 [35,] 9.999881e-01 2.375225e-05 1.187613e-05 [36,] 9.999867e-01 2.668547e-05 1.334274e-05 [37,] 9.999780e-01 4.391211e-05 2.195606e-05 [38,] 9.999492e-01 1.016421e-04 5.082106e-05 [39,] 9.999323e-01 1.354033e-04 6.770165e-05 [40,] 9.999406e-01 1.187157e-04 5.935784e-05 [41,] 9.999372e-01 1.256692e-04 6.283459e-05 [42,] 9.999365e-01 1.270641e-04 6.353207e-05 [43,] 9.999801e-01 3.973950e-05 1.986975e-05 [44,] 9.999917e-01 1.652214e-05 8.261071e-06 [45,] 9.999878e-01 2.444687e-05 1.222343e-05 [46,] 9.999878e-01 2.435244e-05 1.217622e-05 [47,] 9.999949e-01 1.024577e-05 5.122885e-06 [48,] 9.999991e-01 1.795903e-06 8.979516e-07 [49,] 9.999967e-01 6.604415e-06 3.302208e-06 [50,] 9.999992e-01 1.685214e-06 8.426069e-07 [51,] 9.999993e-01 1.462253e-06 7.311267e-07 [52,] 9.999989e-01 2.268219e-06 1.134109e-06 [53,] 9.999963e-01 7.456531e-06 3.728266e-06 [54,] 9.999793e-01 4.146254e-05 2.073127e-05 [55,] 9.998681e-01 2.638210e-04 1.319105e-04 [56,] 9.993165e-01 1.366968e-03 6.834842e-04 [57,] 9.976540e-01 4.691912e-03 2.345956e-03 > postscript(file="/var/www/html/freestat/rcomp/tmp/1aj2w1291108641.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/www/html/freestat/rcomp/tmp/2aj2w1291108641.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/www/html/freestat/rcomp/tmp/3aj2w1291108641.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/www/html/freestat/rcomp/tmp/4la1z1291108641.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/www/html/freestat/rcomp/tmp/5la1z1291108641.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 51.474372 38.985168 -164.230340 -191.821336 171.136347 -77.772400 7 8 9 10 11 12 -304.867312 -468.343303 -687.544262 -682.501773 -727.023147 -832.880037 13 14 15 16 17 18 480.702099 344.976209 135.523871 -73.534298 -31.453916 38.773604 19 20 21 22 23 24 -139.661461 -144.322497 118.374573 -164.832238 -353.793021 -242.678808 25 26 27 28 29 30 31.925883 -827.918197 -776.363260 -604.625813 -221.849254 107.942715 31 32 33 34 35 36 211.103063 276.046542 354.279546 393.570691 450.984427 458.132638 37 38 39 40 41 42 717.054472 517.489065 339.134314 712.246663 814.005494 948.083462 43 44 45 46 47 48 902.320396 658.364677 891.038368 745.941065 578.762553 563.502211 49 50 51 52 53 54 426.747149 310.843359 331.259876 227.141189 -7.222423 526.134572 55 56 57 58 59 60 491.877270 376.866176 304.328662 127.633840 -126.225291 -344.543352 61 62 63 64 65 66 -454.120932 -535.620684 -604.468471 -689.490467 -727.072049 -832.253560 67 68 69 70 -756.626624 -757.040172 -770.451409 -851.554474 > postscript(file="/var/www/html/freestat/rcomp/tmp/6la1z1291108641.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 51.474372 NA 1 38.985168 51.474372 2 -164.230340 38.985168 3 -191.821336 -164.230340 4 171.136347 -191.821336 5 -77.772400 171.136347 6 -304.867312 -77.772400 7 -468.343303 -304.867312 8 -687.544262 -468.343303 9 -682.501773 -687.544262 10 -727.023147 -682.501773 11 -832.880037 -727.023147 12 480.702099 -832.880037 13 344.976209 480.702099 14 135.523871 344.976209 15 -73.534298 135.523871 16 -31.453916 -73.534298 17 38.773604 -31.453916 18 -139.661461 38.773604 19 -144.322497 -139.661461 20 118.374573 -144.322497 21 -164.832238 118.374573 22 -353.793021 -164.832238 23 -242.678808 -353.793021 24 31.925883 -242.678808 25 -827.918197 31.925883 26 -776.363260 -827.918197 27 -604.625813 -776.363260 28 -221.849254 -604.625813 29 107.942715 -221.849254 30 211.103063 107.942715 31 276.046542 211.103063 32 354.279546 276.046542 33 393.570691 354.279546 34 450.984427 393.570691 35 458.132638 450.984427 36 717.054472 458.132638 37 517.489065 717.054472 38 339.134314 517.489065 39 712.246663 339.134314 40 814.005494 712.246663 41 948.083462 814.005494 42 902.320396 948.083462 43 658.364677 902.320396 44 891.038368 658.364677 45 745.941065 891.038368 46 578.762553 745.941065 47 563.502211 578.762553 48 426.747149 563.502211 49 310.843359 426.747149 50 331.259876 310.843359 51 227.141189 331.259876 52 -7.222423 227.141189 53 526.134572 -7.222423 54 491.877270 526.134572 55 376.866176 491.877270 56 304.328662 376.866176 57 127.633840 304.328662 58 -126.225291 127.633840 59 -344.543352 -126.225291 60 -454.120932 -344.543352 61 -535.620684 -454.120932 62 -604.468471 -535.620684 63 -689.490467 -604.468471 64 -727.072049 -689.490467 65 -832.253560 -727.072049 66 -756.626624 -832.253560 67 -757.040172 -756.626624 68 -770.451409 -757.040172 69 -851.554474 -770.451409 70 NA -851.554474 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] 38.985168 51.474372 [2,] -164.230340 38.985168 [3,] -191.821336 -164.230340 [4,] 171.136347 -191.821336 [5,] -77.772400 171.136347 [6,] -304.867312 -77.772400 [7,] -468.343303 -304.867312 [8,] -687.544262 -468.343303 [9,] -682.501773 -687.544262 [10,] -727.023147 -682.501773 [11,] -832.880037 -727.023147 [12,] 480.702099 -832.880037 [13,] 344.976209 480.702099 [14,] 135.523871 344.976209 [15,] -73.534298 135.523871 [16,] -31.453916 -73.534298 [17,] 38.773604 -31.453916 [18,] -139.661461 38.773604 [19,] -144.322497 -139.661461 [20,] 118.374573 -144.322497 [21,] -164.832238 118.374573 [22,] -353.793021 -164.832238 [23,] -242.678808 -353.793021 [24,] 31.925883 -242.678808 [25,] -827.918197 31.925883 [26,] -776.363260 -827.918197 [27,] -604.625813 -776.363260 [28,] -221.849254 -604.625813 [29,] 107.942715 -221.849254 [30,] 211.103063 107.942715 [31,] 276.046542 211.103063 [32,] 354.279546 276.046542 [33,] 393.570691 354.279546 [34,] 450.984427 393.570691 [35,] 458.132638 450.984427 [36,] 717.054472 458.132638 [37,] 517.489065 717.054472 [38,] 339.134314 517.489065 [39,] 712.246663 339.134314 [40,] 814.005494 712.246663 [41,] 948.083462 814.005494 [42,] 902.320396 948.083462 [43,] 658.364677 902.320396 [44,] 891.038368 658.364677 [45,] 745.941065 891.038368 [46,] 578.762553 745.941065 [47,] 563.502211 578.762553 [48,] 426.747149 563.502211 [49,] 310.843359 426.747149 [50,] 331.259876 310.843359 [51,] 227.141189 331.259876 [52,] -7.222423 227.141189 [53,] 526.134572 -7.222423 [54,] 491.877270 526.134572 [55,] 376.866176 491.877270 [56,] 304.328662 376.866176 [57,] 127.633840 304.328662 [58,] -126.225291 127.633840 [59,] -344.543352 -126.225291 [60,] -454.120932 -344.543352 [61,] -535.620684 -454.120932 [62,] -604.468471 -535.620684 [63,] -689.490467 -604.468471 [64,] -727.072049 -689.490467 [65,] -832.253560 -727.072049 [66,] -756.626624 -832.253560 [67,] -757.040172 -756.626624 [68,] -770.451409 -757.040172 [69,] -851.554474 -770.451409 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 38.985168 51.474372 2 -164.230340 38.985168 3 -191.821336 -164.230340 4 171.136347 -191.821336 5 -77.772400 171.136347 6 -304.867312 -77.772400 7 -468.343303 -304.867312 8 -687.544262 -468.343303 9 -682.501773 -687.544262 10 -727.023147 -682.501773 11 -832.880037 -727.023147 12 480.702099 -832.880037 13 344.976209 480.702099 14 135.523871 344.976209 15 -73.534298 135.523871 16 -31.453916 -73.534298 17 38.773604 -31.453916 18 -139.661461 38.773604 19 -144.322497 -139.661461 20 118.374573 -144.322497 21 -164.832238 118.374573 22 -353.793021 -164.832238 23 -242.678808 -353.793021 24 31.925883 -242.678808 25 -827.918197 31.925883 26 -776.363260 -827.918197 27 -604.625813 -776.363260 28 -221.849254 -604.625813 29 107.942715 -221.849254 30 211.103063 107.942715 31 276.046542 211.103063 32 354.279546 276.046542 33 393.570691 354.279546 34 450.984427 393.570691 35 458.132638 450.984427 36 717.054472 458.132638 37 517.489065 717.054472 38 339.134314 517.489065 39 712.246663 339.134314 40 814.005494 712.246663 41 948.083462 814.005494 42 902.320396 948.083462 43 658.364677 902.320396 44 891.038368 658.364677 45 745.941065 891.038368 46 578.762553 745.941065 47 563.502211 578.762553 48 426.747149 563.502211 49 310.843359 426.747149 50 331.259876 310.843359 51 227.141189 331.259876 52 -7.222423 227.141189 53 526.134572 -7.222423 54 491.877270 526.134572 55 376.866176 491.877270 56 304.328662 376.866176 57 127.633840 304.328662 58 -126.225291 127.633840 59 -344.543352 -126.225291 60 -454.120932 -344.543352 61 -535.620684 -454.120932 62 -604.468471 -535.620684 63 -689.490467 -604.468471 64 -727.072049 -689.490467 65 -832.253560 -727.072049 66 -756.626624 -832.253560 67 -757.040172 -756.626624 68 -770.451409 -757.040172 69 -851.554474 -770.451409 > 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/freestat/rcomp/tmp/7e1ik1291108641.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/www/html/freestat/rcomp/tmp/8osin1291108641.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/www/html/freestat/rcomp/tmp/9osin1291108641.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/www/html/freestat/rcomp/tmp/10osin1291108641.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/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11k2xe1291108641.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/freestat/rcomp/tmp/1263wk1291108641.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/freestat/rcomp/tmp/13kdub1291108641.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/freestat/rcomp/tmp/14ndsy1291108641.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/freestat/rcomp/tmp/15g4911291108641.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/freestat/rcomp/tmp/16uepa1291108641.tab") + } > > try(system("convert tmp/1aj2w1291108641.ps tmp/1aj2w1291108641.png",intern=TRUE)) character(0) > try(system("convert tmp/2aj2w1291108641.ps tmp/2aj2w1291108641.png",intern=TRUE)) character(0) > try(system("convert tmp/3aj2w1291108641.ps tmp/3aj2w1291108641.png",intern=TRUE)) character(0) > try(system("convert tmp/4la1z1291108641.ps tmp/4la1z1291108641.png",intern=TRUE)) character(0) > try(system("convert tmp/5la1z1291108641.ps tmp/5la1z1291108641.png",intern=TRUE)) character(0) > try(system("convert tmp/6la1z1291108641.ps tmp/6la1z1291108641.png",intern=TRUE)) character(0) > try(system("convert tmp/7e1ik1291108641.ps tmp/7e1ik1291108641.png",intern=TRUE)) character(0) > try(system("convert tmp/8osin1291108641.ps tmp/8osin1291108641.png",intern=TRUE)) character(0) > try(system("convert tmp/9osin1291108641.ps tmp/9osin1291108641.png",intern=TRUE)) character(0) > try(system("convert tmp/10osin1291108641.ps tmp/10osin1291108641.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.953 2.448 4.294