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(591 + ,0 + ,595 + ,594 + ,611 + ,613 + ,562 + ,519 + ,589 + ,0 + ,591 + ,595 + ,594 + ,611 + ,561 + ,517 + ,584 + ,0 + ,589 + ,591 + ,595 + ,594 + ,555 + ,510 + ,573 + ,0 + ,584 + ,589 + ,591 + ,595 + ,544 + ,509 + ,567 + ,0 + ,573 + ,584 + ,589 + ,591 + ,537 + ,501 + ,569 + ,0 + ,567 + ,573 + ,584 + ,589 + ,543 + ,507 + ,621 + ,0 + ,569 + ,567 + ,573 + ,584 + ,594 + ,569 + ,629 + ,0 + ,621 + ,569 + ,567 + ,573 + ,611 + ,580 + ,628 + ,0 + ,629 + ,621 + ,569 + ,567 + ,613 + ,578 + ,612 + ,0 + ,628 + ,629 + ,621 + ,569 + ,611 + ,565 + ,595 + ,0 + ,612 + ,628 + ,629 + ,621 + ,594 + ,547 + ,597 + ,0 + ,595 + ,612 + ,628 + ,629 + ,595 + ,555 + ,593 + ,0 + ,597 + ,595 + ,612 + ,628 + ,591 + ,562 + ,590 + ,0 + ,593 + ,597 + ,595 + ,612 + ,589 + ,561 + ,580 + ,0 + ,590 + ,593 + ,597 + ,595 + ,584 + ,555 + ,574 + ,0 + ,580 + ,590 + ,593 + ,597 + ,573 + ,544 + ,573 + ,0 + ,574 + ,580 + ,590 + ,593 + ,567 + ,537 + ,573 + ,0 + ,573 + ,574 + ,580 + ,590 + ,569 + ,543 + ,620 + ,0 + ,573 + ,573 + ,574 + ,580 + ,621 + ,594 + ,626 + ,0 + ,620 + ,573 + ,573 + ,574 + ,629 + ,611 + ,620 + ,0 + ,626 + ,620 + ,573 + ,573 + ,628 + ,613 + ,588 + ,0 + ,620 + ,626 + ,620 + ,573 + ,612 + ,611 + ,566 + ,0 + ,588 + ,620 + ,626 + ,620 + ,595 + ,594 + ,557 + ,0 + ,566 + ,588 + ,620 + ,626 + ,597 + ,595 + ,561 + ,0 + ,557 + ,566 + ,588 + ,620 + ,593 + ,591 + ,549 + ,0 + ,561 + ,557 + ,566 + ,588 + ,590 + ,589 + ,532 + ,0 + ,549 + ,561 + ,557 + ,566 + ,580 + ,584 + ,526 + ,0 + ,532 + ,549 + ,561 + ,557 + ,574 + ,573 + ,511 + ,0 + ,526 + ,532 + ,549 + ,561 + ,573 + ,567 + ,499 + ,0 + ,511 + ,526 + ,532 + ,549 + ,573 + ,569 + ,555 + ,1 + ,499 + ,511 + ,526 + ,532 + ,620 + ,621 + ,565 + ,1 + ,555 + ,499 + ,511 + ,526 + ,626 + ,629 + ,542 + ,1 + ,565 + ,555 + ,499 + ,511 + ,620 + ,628 + ,527 + ,1 + ,542 + ,565 + ,555 + ,499 + ,588 + ,612 + ,510 + ,1 + ,527 + ,542 + ,565 + ,555 + ,566 + ,595 + ,514 + ,1 + ,510 + ,527 + ,542 + ,565 + ,557 + ,597 + ,517 + ,1 + ,514 + ,510 + ,527 + ,542 + ,561 + ,593 + ,508 + ,1 + ,517 + ,514 + ,510 + ,527 + ,549 + ,590 + ,493 + ,1 + ,508 + ,517 + ,514 + ,510 + ,532 + ,580 + ,490 + ,1 + ,493 + ,508 + ,517 + ,514 + ,526 + ,574 + ,469 + ,1 + ,490 + ,493 + ,508 + ,517 + ,511 + ,573 + ,478 + ,1 + ,469 + ,490 + ,493 + ,508 + ,499 + ,573 + ,528 + ,1 + ,478 + ,469 + ,490 + ,493 + ,555 + ,620 + ,534 + ,1 + ,528 + ,478 + ,469 + ,490 + ,565 + ,626 + ,518 + ,1 + ,534 + ,528 + ,478 + ,469 + ,542 + ,620 + ,506 + ,1 + ,518 + ,534 + ,528 + ,478 + ,527 + ,588 + ,502 + ,1 + ,506 + ,518 + ,534 + ,528 + ,510 + ,566 + ,516 + ,1 + ,502 + ,506 + ,518 + ,534 + ,514 + ,557 + ,528 + ,1 + ,516 + ,502 + ,506 + ,518 + ,517 + ,561 + ,533 + ,1 + ,528 + ,516 + ,502 + ,506 + ,508 + ,549 + ,536 + ,1 + ,533 + ,528 + ,516 + ,502 + ,493 + ,532 + ,537 + ,1 + ,536 + ,533 + ,528 + ,516 + ,490 + ,526 + ,524 + ,1 + ,537 + ,536 + ,533 + ,528 + ,469 + ,511 + ,536 + ,1 + ,524 + ,537 + ,536 + ,533 + ,478 + ,499 + ,587 + ,1 + ,536 + ,524 + ,537 + ,536 + ,528 + ,555 + ,597 + ,1 + ,587 + ,536 + ,524 + ,537 + ,534 + ,565 + ,581 + ,1 + ,597 + ,587 + ,536 + ,524 + ,518 + ,542 + ,564 + ,1 + ,581 + ,597 + ,587 + ,536 + ,506 + ,527 + ,558 + ,1 + ,564 + ,581 + ,597 + ,587 + ,502 + ,510 + ,575 + ,1 + ,558 + ,564 + ,581 + ,597 + ,516 + ,514) + ,dim=c(8 + ,60) + ,dimnames=list(c('Totaal' + ,'crisis' + ,'t-1' + ,'t-2' + ,'t-3' + ,'t-4' + ,'t-12' + ,'t-24 ') + ,1:60)) > y <- array(NA,dim=c(8,60),dimnames=list(c('Totaal','crisis','t-1','t-2','t-3','t-4','t-12','t-24 '),1:60)) > 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' > #'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 Totaal crisis t-1 t-2 t-3 t-4 t-12 t-24\r 1 591 0 595 594 611 613 562 519 2 589 0 591 595 594 611 561 517 3 584 0 589 591 595 594 555 510 4 573 0 584 589 591 595 544 509 5 567 0 573 584 589 591 537 501 6 569 0 567 573 584 589 543 507 7 621 0 569 567 573 584 594 569 8 629 0 621 569 567 573 611 580 9 628 0 629 621 569 567 613 578 10 612 0 628 629 621 569 611 565 11 595 0 612 628 629 621 594 547 12 597 0 595 612 628 629 595 555 13 593 0 597 595 612 628 591 562 14 590 0 593 597 595 612 589 561 15 580 0 590 593 597 595 584 555 16 574 0 580 590 593 597 573 544 17 573 0 574 580 590 593 567 537 18 573 0 573 574 580 590 569 543 19 620 0 573 573 574 580 621 594 20 626 0 620 573 573 574 629 611 21 620 0 626 620 573 573 628 613 22 588 0 620 626 620 573 612 611 23 566 0 588 620 626 620 595 594 24 557 0 566 588 620 626 597 595 25 561 0 557 566 588 620 593 591 26 549 0 561 557 566 588 590 589 27 532 0 549 561 557 566 580 584 28 526 0 532 549 561 557 574 573 29 511 0 526 532 549 561 573 567 30 499 0 511 526 532 549 573 569 31 555 1 499 511 526 532 620 621 32 565 1 555 499 511 526 626 629 33 542 1 565 555 499 511 620 628 34 527 1 542 565 555 499 588 612 35 510 1 527 542 565 555 566 595 36 514 1 510 527 542 565 557 597 37 517 1 514 510 527 542 561 593 38 508 1 517 514 510 527 549 590 39 493 1 508 517 514 510 532 580 40 490 1 493 508 517 514 526 574 41 469 1 490 493 508 517 511 573 42 478 1 469 490 493 508 499 573 43 528 1 478 469 490 493 555 620 44 534 1 528 478 469 490 565 626 45 518 1 534 528 478 469 542 620 46 506 1 518 534 528 478 527 588 47 502 1 506 518 534 528 510 566 48 516 1 502 506 518 534 514 557 49 528 1 516 502 506 518 517 561 50 533 1 528 516 502 506 508 549 51 536 1 533 528 516 502 493 532 52 537 1 536 533 528 516 490 526 53 524 1 537 536 533 528 469 511 54 536 1 524 537 536 533 478 499 55 587 1 536 524 537 536 528 555 56 597 1 587 536 524 537 534 565 57 581 1 597 587 536 524 518 542 58 564 1 581 597 587 536 506 527 59 558 1 564 581 597 587 502 510 60 575 1 558 564 581 597 516 514 > k <- length(x[1,]) > df <- as.data.frame(x) > (mylm <- lm(df)) Call: lm(formula = df) Coefficients: (Intercept) crisis `t-1` `t-2` `t-3` `t-4` 121.35039 15.44953 1.07951 -0.50928 0.07596 0.06722 `t-12` `t-24\r` 0.44007 -0.38210 > (mysum <- summary(mylm)) Call: lm(formula = df) Residuals: Min 1Q Median 3Q Max -27.214 -6.920 -1.709 4.263 47.400 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 121.35039 87.27716 1.390 0.17033 crisis 15.44953 8.97103 1.722 0.09098 . `t-1` 1.07951 0.14209 7.598 5.53e-10 *** `t-2` -0.50928 0.20721 -2.458 0.01735 * `t-3` 0.07596 0.20522 0.370 0.71280 `t-4` 0.06722 0.14843 0.453 0.65254 `t-12` 0.44007 0.15909 2.766 0.00783 ** `t-24\r` -0.38210 0.15264 -2.503 0.01549 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 16.2 on 52 degrees of freedom Multiple R-squared: 0.8599, Adjusted R-squared: 0.8411 F-statistic: 45.61 on 7 and 52 DF, p-value: < 2.2e-16 > if (n > n25) { + kp3 <- k + 3 + nmkm3 <- n - k - 3 + gqarr <- array(NA, dim=c(nmkm3-kp3+1,3)) + numgqtests <- 0 + numsignificant1 <- 0 + numsignificant5 <- 0 + numsignificant10 <- 0 + for (mypoint in kp3:nmkm3) { + j <- 0 + numgqtests <- numgqtests + 1 + for (myalt in c('greater', 'two.sided', 'less')) { + j <- j + 1 + gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value + } + if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1 + if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1 + if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1 + } + gqarr + } [,1] [,2] [,3] [1,] 0.110948631 0.221897262 0.88905137 [2,] 0.048232948 0.096465895 0.95176705 [3,] 0.019550305 0.039100610 0.98044970 [4,] 0.012738650 0.025477299 0.98726135 [5,] 0.015736962 0.031473924 0.98426304 [6,] 0.011679301 0.023358601 0.98832070 [7,] 0.008131131 0.016262262 0.99186887 [8,] 0.006030130 0.012060260 0.99396987 [9,] 0.009405712 0.018811424 0.99059429 [10,] 0.004393253 0.008786506 0.99560675 [11,] 0.004333118 0.008666235 0.99566688 [12,] 0.004849237 0.009698474 0.99515076 [13,] 0.002414770 0.004829540 0.99758523 [14,] 0.002075811 0.004151623 0.99792419 [15,] 0.002606788 0.005213576 0.99739321 [16,] 0.016250004 0.032500008 0.98375000 [17,] 0.034000635 0.068001269 0.96599937 [18,] 0.037293817 0.074587634 0.96270618 [19,] 0.088727378 0.177454756 0.91127262 [20,] 0.110324888 0.220649776 0.88967511 [21,] 0.175724343 0.351448686 0.82427566 [22,] 0.252778693 0.505557387 0.74722131 [23,] 0.255466219 0.510932439 0.74453378 [24,] 0.266284582 0.532569165 0.73371542 [25,] 0.278588352 0.557176704 0.72141165 [26,] 0.302086138 0.604172276 0.69791386 [27,] 0.295089748 0.590179496 0.70491025 [28,] 0.320700927 0.641401853 0.67929907 [29,] 0.407317994 0.814635988 0.59268201 [30,] 0.497682566 0.995365133 0.50231743 [31,] 0.693741499 0.612517001 0.30625850 [32,] 0.747616618 0.504766765 0.25238338 [33,] 0.870739091 0.258521817 0.12926091 [34,] 0.937082321 0.125835359 0.06291768 [35,] 0.930383974 0.139232051 0.06961603 [36,] 0.900177813 0.199644374 0.09982219 [37,] 0.824322291 0.351355417 0.17567771 [38,] 0.715682493 0.568635015 0.28431751 [39,] 0.624371888 0.751256223 0.37562811 > postscript(file="/var/www/html/freestat/rcomp/tmp/1fdh91292753668.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/2fdh91292753668.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/385yd1292753668.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/485yd1292753668.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/585yd1292753668.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 = 60 Frequency = 1 1 2 3 4 5 6 -6.7752642 -2.8463778 -6.6919657 -8.6176628 -4.8449132 -1.8035323 7 8 9 10 11 12 47.3996196 -1.7993722 13.6539222 -5.3636235 -9.1001556 5.2580094 13 14 15 16 17 18 -3.8412165 1.3601464 -6.5398656 -2.4654590 -1.6186532 -1.2211554 19 20 21 22 23 24 43.0005752 1.7179334 14.4482959 -5.3119829 1.5473491 -0.4458667 25 26 27 28 29 30 5.1314731 -11.3920406 -8.7482068 -3.7693406 -22.1598535 -18.1608163 31 32 33 34 35 36 28.4886135 -26.1162381 -27.2138494 -7.7705530 -21.6289751 -1.1167773 37 38 39 40 41 42 -11.6958335 -15.4631711 -14.7206279 -6.2603397 -24.9599975 12.2070338 43 44 45 46 47 48 36.3472199 -7.3561848 4.1875559 2.4866258 -1.4491074 6.3704838 49 50 51 52 53 54 3.4153601 3.0770363 6.1017358 3.5847195 -6.6432586 10.7898951 55 56 57 58 59 60 41.3313179 4.4883806 1.8816351 2.1156000 -2.6040907 8.1257910 > postscript(file="/var/www/html/freestat/rcomp/tmp/6iwyf1292753668.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 = 60 Frequency = 1 lag(myerror, k = 1) myerror 0 -6.7752642 NA 1 -2.8463778 -6.7752642 2 -6.6919657 -2.8463778 3 -8.6176628 -6.6919657 4 -4.8449132 -8.6176628 5 -1.8035323 -4.8449132 6 47.3996196 -1.8035323 7 -1.7993722 47.3996196 8 13.6539222 -1.7993722 9 -5.3636235 13.6539222 10 -9.1001556 -5.3636235 11 5.2580094 -9.1001556 12 -3.8412165 5.2580094 13 1.3601464 -3.8412165 14 -6.5398656 1.3601464 15 -2.4654590 -6.5398656 16 -1.6186532 -2.4654590 17 -1.2211554 -1.6186532 18 43.0005752 -1.2211554 19 1.7179334 43.0005752 20 14.4482959 1.7179334 21 -5.3119829 14.4482959 22 1.5473491 -5.3119829 23 -0.4458667 1.5473491 24 5.1314731 -0.4458667 25 -11.3920406 5.1314731 26 -8.7482068 -11.3920406 27 -3.7693406 -8.7482068 28 -22.1598535 -3.7693406 29 -18.1608163 -22.1598535 30 28.4886135 -18.1608163 31 -26.1162381 28.4886135 32 -27.2138494 -26.1162381 33 -7.7705530 -27.2138494 34 -21.6289751 -7.7705530 35 -1.1167773 -21.6289751 36 -11.6958335 -1.1167773 37 -15.4631711 -11.6958335 38 -14.7206279 -15.4631711 39 -6.2603397 -14.7206279 40 -24.9599975 -6.2603397 41 12.2070338 -24.9599975 42 36.3472199 12.2070338 43 -7.3561848 36.3472199 44 4.1875559 -7.3561848 45 2.4866258 4.1875559 46 -1.4491074 2.4866258 47 6.3704838 -1.4491074 48 3.4153601 6.3704838 49 3.0770363 3.4153601 50 6.1017358 3.0770363 51 3.5847195 6.1017358 52 -6.6432586 3.5847195 53 10.7898951 -6.6432586 54 41.3313179 10.7898951 55 4.4883806 41.3313179 56 1.8816351 4.4883806 57 2.1156000 1.8816351 58 -2.6040907 2.1156000 59 8.1257910 -2.6040907 60 NA 8.1257910 > dum1 <- dum[2:length(myerror),] > dum1 lag(myerror, k = 1) myerror [1,] -2.8463778 -6.7752642 [2,] -6.6919657 -2.8463778 [3,] -8.6176628 -6.6919657 [4,] -4.8449132 -8.6176628 [5,] -1.8035323 -4.8449132 [6,] 47.3996196 -1.8035323 [7,] -1.7993722 47.3996196 [8,] 13.6539222 -1.7993722 [9,] -5.3636235 13.6539222 [10,] -9.1001556 -5.3636235 [11,] 5.2580094 -9.1001556 [12,] -3.8412165 5.2580094 [13,] 1.3601464 -3.8412165 [14,] -6.5398656 1.3601464 [15,] -2.4654590 -6.5398656 [16,] -1.6186532 -2.4654590 [17,] -1.2211554 -1.6186532 [18,] 43.0005752 -1.2211554 [19,] 1.7179334 43.0005752 [20,] 14.4482959 1.7179334 [21,] -5.3119829 14.4482959 [22,] 1.5473491 -5.3119829 [23,] -0.4458667 1.5473491 [24,] 5.1314731 -0.4458667 [25,] -11.3920406 5.1314731 [26,] -8.7482068 -11.3920406 [27,] -3.7693406 -8.7482068 [28,] -22.1598535 -3.7693406 [29,] -18.1608163 -22.1598535 [30,] 28.4886135 -18.1608163 [31,] -26.1162381 28.4886135 [32,] -27.2138494 -26.1162381 [33,] -7.7705530 -27.2138494 [34,] -21.6289751 -7.7705530 [35,] -1.1167773 -21.6289751 [36,] -11.6958335 -1.1167773 [37,] -15.4631711 -11.6958335 [38,] -14.7206279 -15.4631711 [39,] -6.2603397 -14.7206279 [40,] -24.9599975 -6.2603397 [41,] 12.2070338 -24.9599975 [42,] 36.3472199 12.2070338 [43,] -7.3561848 36.3472199 [44,] 4.1875559 -7.3561848 [45,] 2.4866258 4.1875559 [46,] -1.4491074 2.4866258 [47,] 6.3704838 -1.4491074 [48,] 3.4153601 6.3704838 [49,] 3.0770363 3.4153601 [50,] 6.1017358 3.0770363 [51,] 3.5847195 6.1017358 [52,] -6.6432586 3.5847195 [53,] 10.7898951 -6.6432586 [54,] 41.3313179 10.7898951 [55,] 4.4883806 41.3313179 [56,] 1.8816351 4.4883806 [57,] 2.1156000 1.8816351 [58,] -2.6040907 2.1156000 [59,] 8.1257910 -2.6040907 > z <- as.data.frame(dum1) > z lag(myerror, k = 1) myerror 1 -2.8463778 -6.7752642 2 -6.6919657 -2.8463778 3 -8.6176628 -6.6919657 4 -4.8449132 -8.6176628 5 -1.8035323 -4.8449132 6 47.3996196 -1.8035323 7 -1.7993722 47.3996196 8 13.6539222 -1.7993722 9 -5.3636235 13.6539222 10 -9.1001556 -5.3636235 11 5.2580094 -9.1001556 12 -3.8412165 5.2580094 13 1.3601464 -3.8412165 14 -6.5398656 1.3601464 15 -2.4654590 -6.5398656 16 -1.6186532 -2.4654590 17 -1.2211554 -1.6186532 18 43.0005752 -1.2211554 19 1.7179334 43.0005752 20 14.4482959 1.7179334 21 -5.3119829 14.4482959 22 1.5473491 -5.3119829 23 -0.4458667 1.5473491 24 5.1314731 -0.4458667 25 -11.3920406 5.1314731 26 -8.7482068 -11.3920406 27 -3.7693406 -8.7482068 28 -22.1598535 -3.7693406 29 -18.1608163 -22.1598535 30 28.4886135 -18.1608163 31 -26.1162381 28.4886135 32 -27.2138494 -26.1162381 33 -7.7705530 -27.2138494 34 -21.6289751 -7.7705530 35 -1.1167773 -21.6289751 36 -11.6958335 -1.1167773 37 -15.4631711 -11.6958335 38 -14.7206279 -15.4631711 39 -6.2603397 -14.7206279 40 -24.9599975 -6.2603397 41 12.2070338 -24.9599975 42 36.3472199 12.2070338 43 -7.3561848 36.3472199 44 4.1875559 -7.3561848 45 2.4866258 4.1875559 46 -1.4491074 2.4866258 47 6.3704838 -1.4491074 48 3.4153601 6.3704838 49 3.0770363 3.4153601 50 6.1017358 3.0770363 51 3.5847195 6.1017358 52 -6.6432586 3.5847195 53 10.7898951 -6.6432586 54 41.3313179 10.7898951 55 4.4883806 41.3313179 56 1.8816351 4.4883806 57 2.1156000 1.8816351 58 -2.6040907 2.1156000 59 8.1257910 -2.6040907 > 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/7tnx01292753668.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/8tnx01292753668.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/9tnx01292753668.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/104xe41292753668.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/117fvr1292753668.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/12bytx1292753668.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/13ppr61292753668.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/14a8qc1292753668.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/15er6i1292753668.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/16z9n61292753668.tab") + } > > try(system("convert tmp/1fdh91292753668.ps tmp/1fdh91292753668.png",intern=TRUE)) character(0) > try(system("convert tmp/2fdh91292753668.ps tmp/2fdh91292753668.png",intern=TRUE)) character(0) > try(system("convert tmp/385yd1292753668.ps tmp/385yd1292753668.png",intern=TRUE)) character(0) > try(system("convert tmp/485yd1292753668.ps tmp/485yd1292753668.png",intern=TRUE)) character(0) > try(system("convert tmp/585yd1292753668.ps tmp/585yd1292753668.png",intern=TRUE)) character(0) > try(system("convert tmp/6iwyf1292753668.ps tmp/6iwyf1292753668.png",intern=TRUE)) character(0) > try(system("convert tmp/7tnx01292753668.ps tmp/7tnx01292753668.png",intern=TRUE)) character(0) > try(system("convert tmp/8tnx01292753668.ps tmp/8tnx01292753668.png",intern=TRUE)) character(0) > try(system("convert tmp/9tnx01292753668.ps tmp/9tnx01292753668.png",intern=TRUE)) character(0) > try(system("convert tmp/104xe41292753668.ps tmp/104xe41292753668.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.844 2.416 4.166