R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(507
+ ,104.5
+ ,501
+ ,517
+ ,519
+ ,569
+ ,87.4
+ ,507
+ ,510
+ ,517
+ ,580
+ ,89.9
+ ,569
+ ,509
+ ,510
+ ,578
+ ,109.8
+ ,580
+ ,501
+ ,509
+ ,565
+ ,111.7
+ ,578
+ ,507
+ ,501
+ ,547
+ ,98.6
+ ,565
+ ,569
+ ,507
+ ,555
+ ,96.9
+ ,547
+ ,580
+ ,569
+ ,562
+ ,95.1
+ ,555
+ ,578
+ ,580
+ ,561
+ ,97
+ ,562
+ ,565
+ ,578
+ ,555
+ ,112.7
+ ,561
+ ,547
+ ,565
+ ,544
+ ,102.9
+ ,555
+ ,555
+ ,547
+ ,537
+ ,97.4
+ ,544
+ ,562
+ ,555
+ ,543
+ ,111.4
+ ,537
+ ,561
+ ,562
+ ,594
+ ,87.4
+ ,543
+ ,555
+ ,561
+ ,611
+ ,96.8
+ ,594
+ ,544
+ ,555
+ ,613
+ ,114.1
+ ,611
+ ,537
+ ,544
+ ,611
+ ,110.3
+ ,613
+ ,543
+ ,537
+ ,594
+ ,103.9
+ ,611
+ ,594
+ ,543
+ ,595
+ ,101.6
+ ,594
+ ,611
+ ,594
+ ,591
+ ,94.6
+ ,595
+ ,613
+ ,611
+ ,589
+ ,95.9
+ ,591
+ ,611
+ ,613
+ ,584
+ ,104.7
+ ,589
+ ,594
+ ,611
+ ,573
+ ,102.8
+ ,584
+ ,595
+ ,594
+ ,567
+ ,98.1
+ ,573
+ ,591
+ ,595
+ ,569
+ ,113.9
+ ,567
+ ,589
+ ,591
+ ,621
+ ,80.9
+ ,569
+ ,584
+ ,589
+ ,629
+ ,95.7
+ ,621
+ ,573
+ ,584
+ ,628
+ ,113.2
+ ,629
+ ,567
+ ,573
+ ,612
+ ,105.9
+ ,628
+ ,569
+ ,567
+ ,595
+ ,108.8
+ ,612
+ ,621
+ ,569
+ ,597
+ ,102.3
+ ,595
+ ,629
+ ,621
+ ,593
+ ,99
+ ,597
+ ,628
+ ,629
+ ,590
+ ,100.7
+ ,593
+ ,612
+ ,628
+ ,580
+ ,115.5
+ ,590
+ ,595
+ ,612
+ ,574
+ ,100.7
+ ,580
+ ,597
+ ,595
+ ,573
+ ,109.9
+ ,574
+ ,593
+ ,597
+ ,573
+ ,114.6
+ ,573
+ ,590
+ ,593
+ ,620
+ ,85.4
+ ,573
+ ,580
+ ,590
+ ,626
+ ,100.5
+ ,620
+ ,574
+ ,580
+ ,620
+ ,114.8
+ ,626
+ ,573
+ ,574
+ ,588
+ ,116.5
+ ,620
+ ,573
+ ,573
+ ,566
+ ,112.9
+ ,588
+ ,620
+ ,573
+ ,557
+ ,102
+ ,566
+ ,626
+ ,620
+ ,561
+ ,106
+ ,557
+ ,620
+ ,626
+ ,549
+ ,105.3
+ ,561
+ ,588
+ ,620
+ ,532
+ ,118.8
+ ,549
+ ,566
+ ,588
+ ,526
+ ,106.1
+ ,532
+ ,557
+ ,566
+ ,511
+ ,109.3
+ ,526
+ ,561
+ ,557
+ ,499
+ ,117.2
+ ,511
+ ,549
+ ,561
+ ,555
+ ,92.5
+ ,499
+ ,532
+ ,549
+ ,565
+ ,104.2
+ ,555
+ ,526
+ ,532
+ ,542
+ ,112.5
+ ,565
+ ,511
+ ,526
+ ,527
+ ,122.4
+ ,542
+ ,499
+ ,511
+ ,510
+ ,113.3
+ ,527
+ ,555
+ ,499
+ ,514
+ ,100
+ ,510
+ ,565
+ ,555
+ ,517
+ ,110.7
+ ,514
+ ,542
+ ,565
+ ,508
+ ,112.8
+ ,517
+ ,527
+ ,542
+ ,493
+ ,109.8
+ ,508
+ ,510
+ ,527
+ ,490
+ ,117.3
+ ,493
+ ,514
+ ,510
+ ,469
+ ,109.1
+ ,490
+ ,517
+ ,514
+ ,478
+ ,115.9
+ ,469
+ ,508
+ ,517
+ ,528
+ ,96
+ ,478
+ ,493
+ ,508
+ ,534
+ ,99.8
+ ,528
+ ,490
+ ,493
+ ,518
+ ,116.8
+ ,534
+ ,469
+ ,490
+ ,506
+ ,115.7
+ ,518
+ ,478
+ ,469
+ ,502
+ ,99.4
+ ,506
+ ,528
+ ,478
+ ,516
+ ,94.3
+ ,502
+ ,534
+ ,528)
+ ,dim=c(5
+ ,67)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Yt-1'
+ ,'Yt-4'
+ ,'Yt-5')
+ ,1:67))
> y <- array(NA,dim=c(5,67),dimnames=list(c('Y','X','Yt-1','Yt-4','Yt-5'),1:67))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X Yt-1 Yt-4 Yt-5 t
1 507 104.5 501 517 519 1
2 569 87.4 507 510 517 2
3 580 89.9 569 509 510 3
4 578 109.8 580 501 509 4
5 565 111.7 578 507 501 5
6 547 98.6 565 569 507 6
7 555 96.9 547 580 569 7
8 562 95.1 555 578 580 8
9 561 97.0 562 565 578 9
10 555 112.7 561 547 565 10
11 544 102.9 555 555 547 11
12 537 97.4 544 562 555 12
13 543 111.4 537 561 562 13
14 594 87.4 543 555 561 14
15 611 96.8 594 544 555 15
16 613 114.1 611 537 544 16
17 611 110.3 613 543 537 17
18 594 103.9 611 594 543 18
19 595 101.6 594 611 594 19
20 591 94.6 595 613 611 20
21 589 95.9 591 611 613 21
22 584 104.7 589 594 611 22
23 573 102.8 584 595 594 23
24 567 98.1 573 591 595 24
25 569 113.9 567 589 591 25
26 621 80.9 569 584 589 26
27 629 95.7 621 573 584 27
28 628 113.2 629 567 573 28
29 612 105.9 628 569 567 29
30 595 108.8 612 621 569 30
31 597 102.3 595 629 621 31
32 593 99.0 597 628 629 32
33 590 100.7 593 612 628 33
34 580 115.5 590 595 612 34
35 574 100.7 580 597 595 35
36 573 109.9 574 593 597 36
37 573 114.6 573 590 593 37
38 620 85.4 573 580 590 38
39 626 100.5 620 574 580 39
40 620 114.8 626 573 574 40
41 588 116.5 620 573 573 41
42 566 112.9 588 620 573 42
43 557 102.0 566 626 620 43
44 561 106.0 557 620 626 44
45 549 105.3 561 588 620 45
46 532 118.8 549 566 588 46
47 526 106.1 532 557 566 47
48 511 109.3 526 561 557 48
49 499 117.2 511 549 561 49
50 555 92.5 499 532 549 50
51 565 104.2 555 526 532 51
52 542 112.5 565 511 526 52
53 527 122.4 542 499 511 53
54 510 113.3 527 555 499 54
55 514 100.0 510 565 555 55
56 517 110.7 514 542 565 56
57 508 112.8 517 527 542 57
58 493 109.8 508 510 527 58
59 490 117.3 493 514 510 59
60 469 109.1 490 517 514 60
61 478 115.9 469 508 517 61
62 528 96.0 478 493 508 62
63 534 99.8 528 490 493 63
64 518 116.8 534 469 490 64
65 506 115.7 518 478 469 65
66 502 99.4 506 528 478 66
67 516 94.3 502 534 528 67
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X `Yt-1` `Yt-4` `Yt-5` t
239.71051 -1.50324 0.91378 -0.31105 0.25319 -0.03681
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-21.578 -8.567 -1.920 8.315 26.680
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 239.71051 32.72178 7.326 6.38e-10 ***
X -1.50324 0.18744 -8.020 4.07e-11 ***
`Yt-1` 0.91378 0.05281 17.302 < 2e-16 ***
`Yt-4` -0.31105 0.08682 -3.583 0.000675 ***
`Yt-5` 0.25319 0.08215 3.082 0.003084 **
t -0.03681 0.09690 -0.380 0.705311
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12.44 on 61 degrees of freedom
Multiple R-squared: 0.9125, Adjusted R-squared: 0.9053
F-statistic: 127.2 on 5 and 61 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.7121089 0.5757822 0.2878911
[2,] 0.5614256 0.8771488 0.4385744
[3,] 0.5906025 0.8187951 0.4093975
[4,] 0.6010956 0.7978088 0.3989044
[5,] 0.6663693 0.6672613 0.3336307
[6,] 0.7390506 0.5218987 0.2609494
[7,] 0.6567725 0.6864551 0.3432275
[8,] 0.6903235 0.6193530 0.3096765
[9,] 0.6152095 0.7695809 0.3847905
[10,] 0.5395641 0.9208718 0.4604359
[11,] 0.4880227 0.9760455 0.5119773
[12,] 0.4900966 0.9801932 0.5099034
[13,] 0.4745520 0.9491039 0.5254480
[14,] 0.4389701 0.8779401 0.5610299
[15,] 0.5206114 0.9587773 0.4793886
[16,] 0.7336936 0.5326129 0.2663064
[17,] 0.6917393 0.6165215 0.3082607
[18,] 0.6317468 0.7365064 0.3682532
[19,] 0.6182500 0.7634999 0.3817500
[20,] 0.5731598 0.8536804 0.4268402
[21,] 0.6649574 0.6700852 0.3350426
[22,] 0.6033826 0.7932348 0.3966174
[23,] 0.5335334 0.9329333 0.4664666
[24,] 0.5072593 0.9854813 0.4927407
[25,] 0.5228465 0.9543070 0.4771535
[26,] 0.4550523 0.9101046 0.5449477
[27,] 0.6094042 0.7811915 0.3905958
[28,] 0.5369467 0.9261067 0.4630533
[29,] 0.4875440 0.9750880 0.5124560
[30,] 0.4146573 0.8293147 0.5853427
[31,] 0.3559025 0.7118050 0.6440975
[32,] 0.4387869 0.8775737 0.5612131
[33,] 0.4757332 0.9514663 0.5242668
[34,] 0.4741273 0.9482546 0.5258727
[35,] 0.4408484 0.8816969 0.5591516
[36,] 0.5448856 0.9102287 0.4551144
[37,] 0.5908910 0.8182181 0.4091090
[38,] 0.5942645 0.8114710 0.4057355
[39,] 0.5875840 0.8248321 0.4124160
[40,] 0.5979908 0.8040185 0.4020092
[41,] 0.5207979 0.9584042 0.4792021
[42,] 0.5042836 0.9914327 0.4957164
[43,] 0.4711165 0.9422331 0.5288835
[44,] 0.4601881 0.9203761 0.5398119
[45,] 0.3815259 0.7630518 0.6184741
[46,] 0.3793516 0.7587033 0.6206484
[47,] 0.2934818 0.5869635 0.7065182
[48,] 0.3046919 0.6093839 0.6953081
[49,] 0.3924699 0.7849397 0.6075301
[50,] 0.2735380 0.5470760 0.7264620
> postscript(file="/var/www/html/rcomp/tmp/1zglc1258725557.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2dw4x1258725557.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3zyem1258725557.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/40d8e1258725557.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5zglo1258725557.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 67
Frequency = 1
1 2 3 4 5 6
-3.98253488 25.19527120 -15.20308572 0.46128544 -3.92636147 -11.93670409
7 8 9 10 11 12
-2.28348509 -8.66994553 -16.71075308 -0.46681779 -13.63319244 -18.66071916
13 14 15 16 17 18
12.73449642 20.59783821 3.25971453 16.37592283 10.51151163 0.09967075
19 20 21 22 23 24
5.58855932 -13.49317823 -10.97550807 -5.66414972 -10.29929510 -14.77347255
25 26 27 28 29 30
16.88781799 16.44141099 -2.94621238 16.00574006 -7.87604292 9.80897745
31 32 33 34 35 36
6.93162342 -6.15637340 -7.63256082 6.15664814 -7.99027190 8.60842291
37 38 39 40 41 42
16.70383171 17.49521821 3.94868479 15.20715831 -8.46464000 8.02099561
43 44 45 46 47 48
-7.25783642 7.63052023 -17.47457121 -1.91974765 -8.66901125 -9.81623380
49 50 51 52 53 54
-0.94249064 26.68019618 5.57094272 -17.19984971 3.80122893 7.32249779
55 56 57 58 59 60
-4.16751414 1.61271162 -5.77743987 -18.51631860 9.04993422 -21.57803922
61 62 63 64 65 66
13.31118566 22.82249844 -8.25284459 -9.91623265 -0.79600836 -5.02267423
67
-5.79039902
> postscript(file="/var/www/html/rcomp/tmp/6qjf81258725557.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 -3.98253488 NA
1 25.19527120 -3.98253488
2 -15.20308572 25.19527120
3 0.46128544 -15.20308572
4 -3.92636147 0.46128544
5 -11.93670409 -3.92636147
6 -2.28348509 -11.93670409
7 -8.66994553 -2.28348509
8 -16.71075308 -8.66994553
9 -0.46681779 -16.71075308
10 -13.63319244 -0.46681779
11 -18.66071916 -13.63319244
12 12.73449642 -18.66071916
13 20.59783821 12.73449642
14 3.25971453 20.59783821
15 16.37592283 3.25971453
16 10.51151163 16.37592283
17 0.09967075 10.51151163
18 5.58855932 0.09967075
19 -13.49317823 5.58855932
20 -10.97550807 -13.49317823
21 -5.66414972 -10.97550807
22 -10.29929510 -5.66414972
23 -14.77347255 -10.29929510
24 16.88781799 -14.77347255
25 16.44141099 16.88781799
26 -2.94621238 16.44141099
27 16.00574006 -2.94621238
28 -7.87604292 16.00574006
29 9.80897745 -7.87604292
30 6.93162342 9.80897745
31 -6.15637340 6.93162342
32 -7.63256082 -6.15637340
33 6.15664814 -7.63256082
34 -7.99027190 6.15664814
35 8.60842291 -7.99027190
36 16.70383171 8.60842291
37 17.49521821 16.70383171
38 3.94868479 17.49521821
39 15.20715831 3.94868479
40 -8.46464000 15.20715831
41 8.02099561 -8.46464000
42 -7.25783642 8.02099561
43 7.63052023 -7.25783642
44 -17.47457121 7.63052023
45 -1.91974765 -17.47457121
46 -8.66901125 -1.91974765
47 -9.81623380 -8.66901125
48 -0.94249064 -9.81623380
49 26.68019618 -0.94249064
50 5.57094272 26.68019618
51 -17.19984971 5.57094272
52 3.80122893 -17.19984971
53 7.32249779 3.80122893
54 -4.16751414 7.32249779
55 1.61271162 -4.16751414
56 -5.77743987 1.61271162
57 -18.51631860 -5.77743987
58 9.04993422 -18.51631860
59 -21.57803922 9.04993422
60 13.31118566 -21.57803922
61 22.82249844 13.31118566
62 -8.25284459 22.82249844
63 -9.91623265 -8.25284459
64 -0.79600836 -9.91623265
65 -5.02267423 -0.79600836
66 -5.79039902 -5.02267423
67 NA -5.79039902
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 25.19527120 -3.98253488
[2,] -15.20308572 25.19527120
[3,] 0.46128544 -15.20308572
[4,] -3.92636147 0.46128544
[5,] -11.93670409 -3.92636147
[6,] -2.28348509 -11.93670409
[7,] -8.66994553 -2.28348509
[8,] -16.71075308 -8.66994553
[9,] -0.46681779 -16.71075308
[10,] -13.63319244 -0.46681779
[11,] -18.66071916 -13.63319244
[12,] 12.73449642 -18.66071916
[13,] 20.59783821 12.73449642
[14,] 3.25971453 20.59783821
[15,] 16.37592283 3.25971453
[16,] 10.51151163 16.37592283
[17,] 0.09967075 10.51151163
[18,] 5.58855932 0.09967075
[19,] -13.49317823 5.58855932
[20,] -10.97550807 -13.49317823
[21,] -5.66414972 -10.97550807
[22,] -10.29929510 -5.66414972
[23,] -14.77347255 -10.29929510
[24,] 16.88781799 -14.77347255
[25,] 16.44141099 16.88781799
[26,] -2.94621238 16.44141099
[27,] 16.00574006 -2.94621238
[28,] -7.87604292 16.00574006
[29,] 9.80897745 -7.87604292
[30,] 6.93162342 9.80897745
[31,] -6.15637340 6.93162342
[32,] -7.63256082 -6.15637340
[33,] 6.15664814 -7.63256082
[34,] -7.99027190 6.15664814
[35,] 8.60842291 -7.99027190
[36,] 16.70383171 8.60842291
[37,] 17.49521821 16.70383171
[38,] 3.94868479 17.49521821
[39,] 15.20715831 3.94868479
[40,] -8.46464000 15.20715831
[41,] 8.02099561 -8.46464000
[42,] -7.25783642 8.02099561
[43,] 7.63052023 -7.25783642
[44,] -17.47457121 7.63052023
[45,] -1.91974765 -17.47457121
[46,] -8.66901125 -1.91974765
[47,] -9.81623380 -8.66901125
[48,] -0.94249064 -9.81623380
[49,] 26.68019618 -0.94249064
[50,] 5.57094272 26.68019618
[51,] -17.19984971 5.57094272
[52,] 3.80122893 -17.19984971
[53,] 7.32249779 3.80122893
[54,] -4.16751414 7.32249779
[55,] 1.61271162 -4.16751414
[56,] -5.77743987 1.61271162
[57,] -18.51631860 -5.77743987
[58,] 9.04993422 -18.51631860
[59,] -21.57803922 9.04993422
[60,] 13.31118566 -21.57803922
[61,] 22.82249844 13.31118566
[62,] -8.25284459 22.82249844
[63,] -9.91623265 -8.25284459
[64,] -0.79600836 -9.91623265
[65,] -5.02267423 -0.79600836
[66,] -5.79039902 -5.02267423
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 25.19527120 -3.98253488
2 -15.20308572 25.19527120
3 0.46128544 -15.20308572
4 -3.92636147 0.46128544
5 -11.93670409 -3.92636147
6 -2.28348509 -11.93670409
7 -8.66994553 -2.28348509
8 -16.71075308 -8.66994553
9 -0.46681779 -16.71075308
10 -13.63319244 -0.46681779
11 -18.66071916 -13.63319244
12 12.73449642 -18.66071916
13 20.59783821 12.73449642
14 3.25971453 20.59783821
15 16.37592283 3.25971453
16 10.51151163 16.37592283
17 0.09967075 10.51151163
18 5.58855932 0.09967075
19 -13.49317823 5.58855932
20 -10.97550807 -13.49317823
21 -5.66414972 -10.97550807
22 -10.29929510 -5.66414972
23 -14.77347255 -10.29929510
24 16.88781799 -14.77347255
25 16.44141099 16.88781799
26 -2.94621238 16.44141099
27 16.00574006 -2.94621238
28 -7.87604292 16.00574006
29 9.80897745 -7.87604292
30 6.93162342 9.80897745
31 -6.15637340 6.93162342
32 -7.63256082 -6.15637340
33 6.15664814 -7.63256082
34 -7.99027190 6.15664814
35 8.60842291 -7.99027190
36 16.70383171 8.60842291
37 17.49521821 16.70383171
38 3.94868479 17.49521821
39 15.20715831 3.94868479
40 -8.46464000 15.20715831
41 8.02099561 -8.46464000
42 -7.25783642 8.02099561
43 7.63052023 -7.25783642
44 -17.47457121 7.63052023
45 -1.91974765 -17.47457121
46 -8.66901125 -1.91974765
47 -9.81623380 -8.66901125
48 -0.94249064 -9.81623380
49 26.68019618 -0.94249064
50 5.57094272 26.68019618
51 -17.19984971 5.57094272
52 3.80122893 -17.19984971
53 7.32249779 3.80122893
54 -4.16751414 7.32249779
55 1.61271162 -4.16751414
56 -5.77743987 1.61271162
57 -18.51631860 -5.77743987
58 9.04993422 -18.51631860
59 -21.57803922 9.04993422
60 13.31118566 -21.57803922
61 22.82249844 13.31118566
62 -8.25284459 22.82249844
63 -9.91623265 -8.25284459
64 -0.79600836 -9.91623265
65 -5.02267423 -0.79600836
66 -5.79039902 -5.02267423
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/7gok11258725557.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/8oheq1258725557.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9ovte1258725557.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10jar91258725557.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/11p0nt1258725557.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/128nwe1258725557.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/13w36s1258725557.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14dem81258725557.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15hb8v1258725557.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/16oh2r1258725558.tab")
+ }
>
> system("convert tmp/1zglc1258725557.ps tmp/1zglc1258725557.png")
> system("convert tmp/2dw4x1258725557.ps tmp/2dw4x1258725557.png")
> system("convert tmp/3zyem1258725557.ps tmp/3zyem1258725557.png")
> system("convert tmp/40d8e1258725557.ps tmp/40d8e1258725557.png")
> system("convert tmp/5zglo1258725557.ps tmp/5zglo1258725557.png")
> system("convert tmp/6qjf81258725557.ps tmp/6qjf81258725557.png")
> system("convert tmp/7gok11258725557.ps tmp/7gok11258725557.png")
> system("convert tmp/8oheq1258725557.ps tmp/8oheq1258725557.png")
> system("convert tmp/9ovte1258725557.ps tmp/9ovte1258725557.png")
> system("convert tmp/10jar91258725557.ps tmp/10jar91258725557.png")
>
>
> proc.time()
user system elapsed
2.613 1.589 5.976