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