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 = 'No Linear Trend' > par2 = 'Include Monthly 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 1 2649.20 31077 0 1 0 0 0 0 0 0 0 0 0 0 2 2579.40 31293 0 0 1 0 0 0 0 0 0 0 0 0 3 2504.60 30236 0 0 0 1 0 0 0 0 0 0 0 0 4 2462.30 30160 0 0 0 0 1 0 0 0 0 0 0 0 5 2467.40 32436 0 0 0 0 0 1 0 0 0 0 0 0 6 2446.70 30695 0 0 0 0 0 0 1 0 0 0 0 0 7 2656.30 27525 0 0 0 0 0 0 0 1 0 0 0 0 8 2626.20 26434 0 0 0 0 0 0 0 0 1 0 0 0 9 2482.60 25739 0 0 0 0 0 0 0 0 0 1 0 0 10 2539.90 25204 0 0 0 0 0 0 0 0 0 0 1 0 11 2502.70 24977 0 0 0 0 0 0 0 0 0 0 0 1 12 2466.90 24320 0 0 0 0 0 0 0 0 0 0 0 0 13 2513.20 22680 1 1 0 0 0 0 0 0 0 0 0 0 14 2443.30 22052 1 0 1 0 0 0 0 0 0 0 0 0 15 2293.40 21467 1 0 0 1 0 0 0 0 0 0 0 0 16 2070.80 21383 1 0 0 0 1 0 0 0 0 0 0 0 17 2029.60 21777 1 0 0 0 0 1 0 0 0 0 0 0 18 2052.00 21928 1 0 0 0 0 0 1 0 0 0 0 0 19 1864.40 21814 1 0 0 0 0 0 0 1 0 0 0 0 20 1670.10 22937 1 0 0 0 0 0 0 0 1 0 0 0 21 1811.00 23595 1 0 0 0 0 0 0 0 0 1 0 0 22 1905.40 20830 1 0 0 0 0 0 0 0 0 0 1 0 23 1862.80 19650 1 0 0 0 0 0 0 0 0 0 0 1 24 2014.50 19195 1 0 0 0 0 0 0 0 0 0 0 0 25 2197.80 19644 1 1 0 0 0 0 0 0 0 0 0 0 26 2962.30 18483 0 0 1 0 0 0 0 0 0 0 0 0 27 3047.00 18079 0 0 0 1 0 0 0 0 0 0 0 0 28 3032.60 19178 0 0 0 0 1 0 0 0 0 0 0 0 29 3504.40 18391 0 0 0 0 0 1 0 0 0 0 0 0 30 3801.10 18441 0 0 0 0 0 0 1 0 0 0 0 0 31 3857.60 18584 0 0 0 0 0 0 0 1 0 0 0 0 32 3674.40 20108 0 0 0 0 0 0 0 0 1 0 0 0 33 3721.00 20148 0 0 0 0 0 0 0 0 0 1 0 0 34 3844.50 19394 0 0 0 0 0 0 0 0 0 0 1 0 35 4116.70 17745 0 0 0 0 0 0 0 0 0 0 0 1 36 4105.20 17696 0 0 0 0 0 0 0 0 0 0 0 0 37 4435.20 17032 0 1 0 0 0 0 0 0 0 0 0 0 38 4296.50 16438 0 0 1 0 0 0 0 0 0 0 0 0 39 4202.50 15683 0 0 0 1 0 0 0 0 0 0 0 0 40 4562.80 15594 0 0 0 0 1 0 0 0 0 0 0 0 41 4621.40 15713 0 0 0 0 0 1 0 0 0 0 0 0 42 4697.00 15937 0 0 0 0 0 0 1 0 0 0 0 0 43 4591.30 16171 0 0 0 0 0 0 0 1 0 0 0 0 44 4357.00 15928 0 0 0 0 0 0 0 0 1 0 0 0 45 4502.60 16348 0 0 0 0 0 0 0 0 0 1 0 0 46 4443.90 15579 0 0 0 0 0 0 0 0 0 0 1 0 47 4290.90 15305 0 0 0 0 0 0 0 0 0 0 0 1 48 4199.80 15648 0 0 0 0 0 0 0 0 0 0 0 0 49 4138.50 14954 0 1 0 0 0 0 0 0 0 0 0 0 50 3970.10 15137 0 0 1 0 0 0 0 0 0 0 0 0 51 3862.30 15839 0 0 0 1 0 0 0 0 0 0 0 0 52 3701.60 16050 0 0 0 0 1 0 0 0 0 0 0 0 53 3570.12 15168 0 0 0 0 0 1 0 0 0 0 0 0 54 3801.06 17064 0 0 0 0 0 0 1 0 0 0 0 0 55 3895.51 16005 0 0 0 0 0 0 0 1 0 0 0 0 56 3917.96 14886 0 0 0 0 0 0 0 0 1 0 0 0 57 3813.06 14931 0 0 0 0 0 0 0 0 0 1 0 0 58 3667.03 14544 0 0 0 0 0 0 0 0 0 0 1 0 59 3494.17 13812 0 0 0 0 0 0 0 0 0 0 0 1 60 3364.00 13031 0 0 0 0 0 0 0 0 0 0 0 0 61 3295.30 12574 0 1 0 0 0 0 0 0 0 0 0 0 62 3277.00 11964 0 0 1 0 0 0 0 0 0 0 0 0 63 3257.20 11451 0 0 0 1 0 0 0 0 0 0 0 0 64 3161.70 11346 0 0 0 0 1 0 0 0 0 0 0 0 65 3097.30 11353 0 0 0 0 0 1 0 0 0 0 0 0 66 3061.30 10702 0 0 0 0 0 0 1 0 0 0 0 0 67 3119.30 10646 0 0 0 0 0 0 0 1 0 0 0 0 68 3106.22 10556 0 0 0 0 0 0 0 0 1 0 0 0 69 3080.58 10463 0 0 0 0 0 0 0 0 0 1 0 0 70 2981.85 10407 0 0 0 0 0 0 0 0 0 0 1 0 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) Goudprijs Crisis M1 M2 M3 4.594e+03 -6.184e-02 -1.262e+03 2.470e+02 5.992e+01 -2.727e+01 M4 M5 M6 M7 M8 M9 -4.661e+01 1.474e+01 1.088e+02 8.825e+01 -1.610e+01 -2.405e+00 M10 M11 -6.139e+01 4.315e+01 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -979.85 -350.85 -46.72 338.89 984.20 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.594e+03 3.476e+02 13.217 < 2e-16 *** Goudprijs -6.184e-02 1.301e-02 -4.754 1.44e-05 *** Crisis -1.262e+03 1.868e+02 -6.755 8.80e-09 *** M1 2.470e+02 3.564e+02 0.693 0.491 M2 5.992e+01 3.557e+02 0.168 0.867 M3 -2.727e+01 3.555e+02 -0.077 0.939 M4 -4.661e+01 3.555e+02 -0.131 0.896 M5 1.474e+01 3.556e+02 0.041 0.967 M6 1.088e+02 3.556e+02 0.306 0.761 M7 8.825e+01 3.553e+02 0.248 0.805 M8 -1.610e+01 3.553e+02 -0.045 0.964 M9 -2.405e+00 3.554e+02 -0.007 0.995 M10 -6.139e+01 3.553e+02 -0.173 0.863 M11 4.315e+01 3.710e+02 0.116 0.908 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 586.6 on 56 degrees of freedom Multiple R-squared: 0.6123, Adjusted R-squared: 0.5223 F-statistic: 6.802 on 13 and 56 DF, p-value: 1.320e-07 > 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,] 8.326970e-03 1.665394e-02 0.9916730296 [2,] 2.629341e-03 5.258682e-03 0.9973706591 [3,] 5.075126e-03 1.015025e-02 0.9949248740 [4,] 1.349320e-03 2.698640e-03 0.9986506802 [5,] 8.232962e-04 1.646592e-03 0.9991767038 [6,] 2.057133e-04 4.114265e-04 0.9997942867 [7,] 4.827436e-05 9.654871e-05 0.9999517256 [8,] 1.554944e-05 3.109889e-05 0.9999844506 [9,] 3.467255e-05 6.934510e-05 0.9999653275 [10,] 3.659998e-05 7.319997e-05 0.9999634000 [11,] 1.994277e-05 3.988553e-05 0.9999800572 [12,] 2.680358e-05 5.360716e-05 0.9999731964 [13,] 1.818609e-04 3.637218e-04 0.9998181391 [14,] 1.975185e-03 3.950370e-03 0.9980248152 [15,] 9.736752e-03 1.947350e-02 0.9902632478 [16,] 4.295706e-02 8.591413e-02 0.9570429352 [17,] 1.670392e-01 3.340784e-01 0.8329607882 [18,] 4.348909e-01 8.697818e-01 0.5651090918 [19,] 6.120389e-01 7.759221e-01 0.3879610642 [20,] 7.064732e-01 5.870537e-01 0.2935268472 [21,] 6.544444e-01 6.911112e-01 0.3455555769 [22,] 6.024958e-01 7.950084e-01 0.3975042103 [23,] 5.440719e-01 9.118562e-01 0.4559281000 [24,] 6.761026e-01 6.477948e-01 0.3238973791 [25,] 7.930440e-01 4.139121e-01 0.2069560467 [26,] 9.255236e-01 1.489528e-01 0.0744764065 [27,] 9.523384e-01 9.532325e-02 0.0476616270 [28,] 9.370385e-01 1.259231e-01 0.0629615264 [29,] 9.430708e-01 1.138583e-01 0.0569291588 [30,] 9.779245e-01 4.415093e-02 0.0220754674 [31,] 9.898728e-01 2.025437e-02 0.0101271838 [32,] 9.934382e-01 1.312367e-02 0.0065618364 [33,] 9.991133e-01 1.773498e-03 0.0008867489 [34,] 9.990554e-01 1.889175e-03 0.0009445874 [35,] 9.959804e-01 8.039184e-03 0.0040195922 [36,] 9.891892e-01 2.162165e-02 0.0108108255 [37,] 9.668537e-01 6.629268e-02 0.0331463396 > postscript(file="/var/www/html/freestat/rcomp/tmp/185jx1291107212.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/285jx1291107212.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/3jei01291107212.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/4jei01291107212.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/5jei01291107212.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 -270.21392 -139.55831 -192.53361 -220.18590 -135.69672 -358.14644 -323.99318 8 9 10 11 12 13 14 -317.20746 -517.47561 -434.27615 -590.05029 -623.32771 336.12926 414.49450 15 16 17 18 19 20 21 315.60624 107.15925 28.97130 -33.38292 -207.45597 -227.96325 -60.06602 22 23 24 25 26 27 28 -77.66295 -297.76769 -131.05404 -167.00770 -548.78955 -401.88533 -328.97921 29 30 31 32 33 34 35 32.80341 238.50366 324.42271 339.81203 375.19404 511.05121 576.74494 36 37 38 39 40 41 42 605.36438 647.28621 658.95390 605.45336 979.59718 984.20407 979.56395 43 44 45 46 47 48 49 908.91016 763.93361 921.81365 874.54327 600.06280 573.32232 222.08904 50 51 52 53 54 55 56 252.10404 274.89992 146.59483 -100.77706 153.31419 202.85523 260.45952 57 58 59 60 61 62 63 144.65070 33.67203 -288.98976 -424.30496 -768.28289 -637.20458 -601.54059 64 65 66 67 68 69 70 -684.18615 -809.50500 -979.85244 -904.73895 -819.03445 -864.11676 -907.32740 > postscript(file="/var/www/html/freestat/rcomp/tmp/6u6il1291107212.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 -270.21392 NA 1 -139.55831 -270.21392 2 -192.53361 -139.55831 3 -220.18590 -192.53361 4 -135.69672 -220.18590 5 -358.14644 -135.69672 6 -323.99318 -358.14644 7 -317.20746 -323.99318 8 -517.47561 -317.20746 9 -434.27615 -517.47561 10 -590.05029 -434.27615 11 -623.32771 -590.05029 12 336.12926 -623.32771 13 414.49450 336.12926 14 315.60624 414.49450 15 107.15925 315.60624 16 28.97130 107.15925 17 -33.38292 28.97130 18 -207.45597 -33.38292 19 -227.96325 -207.45597 20 -60.06602 -227.96325 21 -77.66295 -60.06602 22 -297.76769 -77.66295 23 -131.05404 -297.76769 24 -167.00770 -131.05404 25 -548.78955 -167.00770 26 -401.88533 -548.78955 27 -328.97921 -401.88533 28 32.80341 -328.97921 29 238.50366 32.80341 30 324.42271 238.50366 31 339.81203 324.42271 32 375.19404 339.81203 33 511.05121 375.19404 34 576.74494 511.05121 35 605.36438 576.74494 36 647.28621 605.36438 37 658.95390 647.28621 38 605.45336 658.95390 39 979.59718 605.45336 40 984.20407 979.59718 41 979.56395 984.20407 42 908.91016 979.56395 43 763.93361 908.91016 44 921.81365 763.93361 45 874.54327 921.81365 46 600.06280 874.54327 47 573.32232 600.06280 48 222.08904 573.32232 49 252.10404 222.08904 50 274.89992 252.10404 51 146.59483 274.89992 52 -100.77706 146.59483 53 153.31419 -100.77706 54 202.85523 153.31419 55 260.45952 202.85523 56 144.65070 260.45952 57 33.67203 144.65070 58 -288.98976 33.67203 59 -424.30496 -288.98976 60 -768.28289 -424.30496 61 -637.20458 -768.28289 62 -601.54059 -637.20458 63 -684.18615 -601.54059 64 -809.50500 -684.18615 65 -979.85244 -809.50500 66 -904.73895 -979.85244 67 -819.03445 -904.73895 68 -864.11676 -819.03445 69 -907.32740 -864.11676 70 NA -907.32740 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -139.55831 -270.21392 [2,] -192.53361 -139.55831 [3,] -220.18590 -192.53361 [4,] -135.69672 -220.18590 [5,] -358.14644 -135.69672 [6,] -323.99318 -358.14644 [7,] -317.20746 -323.99318 [8,] -517.47561 -317.20746 [9,] -434.27615 -517.47561 [10,] -590.05029 -434.27615 [11,] -623.32771 -590.05029 [12,] 336.12926 -623.32771 [13,] 414.49450 336.12926 [14,] 315.60624 414.49450 [15,] 107.15925 315.60624 [16,] 28.97130 107.15925 [17,] -33.38292 28.97130 [18,] -207.45597 -33.38292 [19,] -227.96325 -207.45597 [20,] -60.06602 -227.96325 [21,] -77.66295 -60.06602 [22,] -297.76769 -77.66295 [23,] -131.05404 -297.76769 [24,] -167.00770 -131.05404 [25,] -548.78955 -167.00770 [26,] -401.88533 -548.78955 [27,] -328.97921 -401.88533 [28,] 32.80341 -328.97921 [29,] 238.50366 32.80341 [30,] 324.42271 238.50366 [31,] 339.81203 324.42271 [32,] 375.19404 339.81203 [33,] 511.05121 375.19404 [34,] 576.74494 511.05121 [35,] 605.36438 576.74494 [36,] 647.28621 605.36438 [37,] 658.95390 647.28621 [38,] 605.45336 658.95390 [39,] 979.59718 605.45336 [40,] 984.20407 979.59718 [41,] 979.56395 984.20407 [42,] 908.91016 979.56395 [43,] 763.93361 908.91016 [44,] 921.81365 763.93361 [45,] 874.54327 921.81365 [46,] 600.06280 874.54327 [47,] 573.32232 600.06280 [48,] 222.08904 573.32232 [49,] 252.10404 222.08904 [50,] 274.89992 252.10404 [51,] 146.59483 274.89992 [52,] -100.77706 146.59483 [53,] 153.31419 -100.77706 [54,] 202.85523 153.31419 [55,] 260.45952 202.85523 [56,] 144.65070 260.45952 [57,] 33.67203 144.65070 [58,] -288.98976 33.67203 [59,] -424.30496 -288.98976 [60,] -768.28289 -424.30496 [61,] -637.20458 -768.28289 [62,] -601.54059 -637.20458 [63,] -684.18615 -601.54059 [64,] -809.50500 -684.18615 [65,] -979.85244 -809.50500 [66,] -904.73895 -979.85244 [67,] -819.03445 -904.73895 [68,] -864.11676 -819.03445 [69,] -907.32740 -864.11676 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -139.55831 -270.21392 2 -192.53361 -139.55831 3 -220.18590 -192.53361 4 -135.69672 -220.18590 5 -358.14644 -135.69672 6 -323.99318 -358.14644 7 -317.20746 -323.99318 8 -517.47561 -317.20746 9 -434.27615 -517.47561 10 -590.05029 -434.27615 11 -623.32771 -590.05029 12 336.12926 -623.32771 13 414.49450 336.12926 14 315.60624 414.49450 15 107.15925 315.60624 16 28.97130 107.15925 17 -33.38292 28.97130 18 -207.45597 -33.38292 19 -227.96325 -207.45597 20 -60.06602 -227.96325 21 -77.66295 -60.06602 22 -297.76769 -77.66295 23 -131.05404 -297.76769 24 -167.00770 -131.05404 25 -548.78955 -167.00770 26 -401.88533 -548.78955 27 -328.97921 -401.88533 28 32.80341 -328.97921 29 238.50366 32.80341 30 324.42271 238.50366 31 339.81203 324.42271 32 375.19404 339.81203 33 511.05121 375.19404 34 576.74494 511.05121 35 605.36438 576.74494 36 647.28621 605.36438 37 658.95390 647.28621 38 605.45336 658.95390 39 979.59718 605.45336 40 984.20407 979.59718 41 979.56395 984.20407 42 908.91016 979.56395 43 763.93361 908.91016 44 921.81365 763.93361 45 874.54327 921.81365 46 600.06280 874.54327 47 573.32232 600.06280 48 222.08904 573.32232 49 252.10404 222.08904 50 274.89992 252.10404 51 146.59483 274.89992 52 -100.77706 146.59483 53 153.31419 -100.77706 54 202.85523 153.31419 55 260.45952 202.85523 56 144.65070 260.45952 57 33.67203 144.65070 58 -288.98976 33.67203 59 -424.30496 -288.98976 60 -768.28289 -424.30496 61 -637.20458 -768.28289 62 -601.54059 -637.20458 63 -684.18615 -601.54059 64 -809.50500 -684.18615 65 -979.85244 -809.50500 66 -904.73895 -979.85244 67 -819.03445 -904.73895 68 -864.11676 -819.03445 69 -907.32740 -864.11676 > 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/7mfho1291107212.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/8mfho1291107212.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/9mfho1291107212.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/10x6y91291107212.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/110pef1291107212.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/1237v31291107212.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/13n2001291107212.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/14lirz1291107212.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/15p0qn1291107212.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/163ane1291107212.tab") + } > > try(system("convert tmp/185jx1291107212.ps tmp/185jx1291107212.png",intern=TRUE)) character(0) > try(system("convert tmp/285jx1291107212.ps tmp/285jx1291107212.png",intern=TRUE)) character(0) > try(system("convert tmp/3jei01291107212.ps tmp/3jei01291107212.png",intern=TRUE)) character(0) > try(system("convert tmp/4jei01291107212.ps tmp/4jei01291107212.png",intern=TRUE)) character(0) > try(system("convert tmp/5jei01291107212.ps tmp/5jei01291107212.png",intern=TRUE)) character(0) > try(system("convert tmp/6u6il1291107212.ps tmp/6u6il1291107212.png",intern=TRUE)) character(0) > try(system("convert tmp/7mfho1291107212.ps tmp/7mfho1291107212.png",intern=TRUE)) character(0) > try(system("convert tmp/8mfho1291107212.ps tmp/8mfho1291107212.png",intern=TRUE)) character(0) > try(system("convert tmp/9mfho1291107212.ps tmp/9mfho1291107212.png",intern=TRUE)) character(0) > try(system("convert tmp/10x6y91291107212.ps tmp/10x6y91291107212.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.986 2.526 4.451