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(501
+ ,98.1
+ ,509
+ ,510
+ ,517
+ ,519
+ ,507
+ ,104.5
+ ,501
+ ,509
+ ,510
+ ,517
+ ,569
+ ,87.4
+ ,507
+ ,501
+ ,509
+ ,510
+ ,580
+ ,89.9
+ ,569
+ ,507
+ ,501
+ ,509
+ ,578
+ ,109.8
+ ,580
+ ,569
+ ,507
+ ,501
+ ,565
+ ,111.7
+ ,578
+ ,580
+ ,569
+ ,507
+ ,547
+ ,98.6
+ ,565
+ ,578
+ ,580
+ ,569
+ ,555
+ ,96.9
+ ,547
+ ,565
+ ,578
+ ,580
+ ,562
+ ,95.1
+ ,555
+ ,547
+ ,565
+ ,578
+ ,561
+ ,97
+ ,562
+ ,555
+ ,547
+ ,565
+ ,555
+ ,112.7
+ ,561
+ ,562
+ ,555
+ ,547
+ ,544
+ ,102.9
+ ,555
+ ,561
+ ,562
+ ,555
+ ,537
+ ,97.4
+ ,544
+ ,555
+ ,561
+ ,562
+ ,543
+ ,111.4
+ ,537
+ ,544
+ ,555
+ ,561
+ ,594
+ ,87.4
+ ,543
+ ,537
+ ,544
+ ,555
+ ,611
+ ,96.8
+ ,594
+ ,543
+ ,537
+ ,544
+ ,613
+ ,114.1
+ ,611
+ ,594
+ ,543
+ ,537
+ ,611
+ ,110.3
+ ,613
+ ,611
+ ,594
+ ,543
+ ,594
+ ,103.9
+ ,611
+ ,613
+ ,611
+ ,594
+ ,595
+ ,101.6
+ ,594
+ ,611
+ ,613
+ ,611
+ ,591
+ ,94.6
+ ,595
+ ,594
+ ,611
+ ,613
+ ,589
+ ,95.9
+ ,591
+ ,595
+ ,594
+ ,611
+ ,584
+ ,104.7
+ ,589
+ ,591
+ ,595
+ ,594
+ ,573
+ ,102.8
+ ,584
+ ,589
+ ,591
+ ,595
+ ,567
+ ,98.1
+ ,573
+ ,584
+ ,589
+ ,591
+ ,569
+ ,113.9
+ ,567
+ ,573
+ ,584
+ ,589
+ ,621
+ ,80.9
+ ,569
+ ,567
+ ,573
+ ,584
+ ,629
+ ,95.7
+ ,621
+ ,569
+ ,567
+ ,573
+ ,628
+ ,113.2
+ ,629
+ ,621
+ ,569
+ ,567
+ ,612
+ ,105.9
+ ,628
+ ,629
+ ,621
+ ,569
+ ,595
+ ,108.8
+ ,612
+ ,628
+ ,629
+ ,621
+ ,597
+ ,102.3
+ ,595
+ ,612
+ ,628
+ ,629
+ ,593
+ ,99
+ ,597
+ ,595
+ ,612
+ ,628
+ ,590
+ ,100.7
+ ,593
+ ,597
+ ,595
+ ,612
+ ,580
+ ,115.5
+ ,590
+ ,593
+ ,597
+ ,595
+ ,574
+ ,100.7
+ ,580
+ ,590
+ ,593
+ ,597
+ ,573
+ ,109.9
+ ,574
+ ,580
+ ,590
+ ,593
+ ,573
+ ,114.6
+ ,573
+ ,574
+ ,580
+ ,590
+ ,620
+ ,85.4
+ ,573
+ ,573
+ ,574
+ ,580
+ ,626
+ ,100.5
+ ,620
+ ,573
+ ,573
+ ,574
+ ,620
+ ,114.8
+ ,626
+ ,620
+ ,573
+ ,573
+ ,588
+ ,116.5
+ ,620
+ ,626
+ ,620
+ ,573
+ ,566
+ ,112.9
+ ,588
+ ,620
+ ,626
+ ,620
+ ,557
+ ,102
+ ,566
+ ,588
+ ,620
+ ,626
+ ,561
+ ,106
+ ,557
+ ,566
+ ,588
+ ,620
+ ,549
+ ,105.3
+ ,561
+ ,557
+ ,566
+ ,588
+ ,532
+ ,118.8
+ ,549
+ ,561
+ ,557
+ ,566
+ ,526
+ ,106.1
+ ,532
+ ,549
+ ,561
+ ,557
+ ,511
+ ,109.3
+ ,526
+ ,532
+ ,549
+ ,561
+ ,499
+ ,117.2
+ ,511
+ ,526
+ ,532
+ ,549
+ ,555
+ ,92.5
+ ,499
+ ,511
+ ,526
+ ,532
+ ,565
+ ,104.2
+ ,555
+ ,499
+ ,511
+ ,526
+ ,542
+ ,112.5
+ ,565
+ ,555
+ ,499
+ ,511
+ ,527
+ ,122.4
+ ,542
+ ,565
+ ,555
+ ,499
+ ,510
+ ,113.3
+ ,527
+ ,542
+ ,565
+ ,555
+ ,514
+ ,100
+ ,510
+ ,527
+ ,542
+ ,565
+ ,517
+ ,110.7
+ ,514
+ ,510
+ ,527
+ ,542
+ ,508
+ ,112.8
+ ,517
+ ,514
+ ,510
+ ,527
+ ,493
+ ,109.8
+ ,508
+ ,517
+ ,514
+ ,510
+ ,490
+ ,117.3
+ ,493
+ ,508
+ ,517
+ ,514
+ ,469
+ ,109.1
+ ,490
+ ,493
+ ,508
+ ,517
+ ,478
+ ,115.9
+ ,469
+ ,490
+ ,493
+ ,508
+ ,528
+ ,96
+ ,478
+ ,469
+ ,490
+ ,493
+ ,534
+ ,99.8
+ ,528
+ ,478
+ ,469
+ ,490
+ ,518
+ ,116.8
+ ,534
+ ,528
+ ,478
+ ,469
+ ,506
+ ,115.7
+ ,518
+ ,534
+ ,528
+ ,478
+ ,502
+ ,99.4
+ ,506
+ ,518
+ ,534
+ ,528
+ ,516
+ ,94.3
+ ,502
+ ,506
+ ,518
+ ,534)
+ ,dim=c(6
+ ,68)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Yt-1'
+ ,'Yt-2'
+ ,'Yt-3'
+ ,'Yt-4')
+ ,1:68))
> y <- array(NA,dim=c(6,68),dimnames=list(c('Y','X','Yt-1','Yt-2','Yt-3','Yt-4'),1:68))
> 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-2 Yt-3 Yt-4 t
1 501 98.1 509 510 517 519 1
2 507 104.5 501 509 510 517 2
3 569 87.4 507 501 509 510 3
4 580 89.9 569 507 501 509 4
5 578 109.8 580 569 507 501 5
6 565 111.7 578 580 569 507 6
7 547 98.6 565 578 580 569 7
8 555 96.9 547 565 578 580 8
9 562 95.1 555 547 565 578 9
10 561 97.0 562 555 547 565 10
11 555 112.7 561 562 555 547 11
12 544 102.9 555 561 562 555 12
13 537 97.4 544 555 561 562 13
14 543 111.4 537 544 555 561 14
15 594 87.4 543 537 544 555 15
16 611 96.8 594 543 537 544 16
17 613 114.1 611 594 543 537 17
18 611 110.3 613 611 594 543 18
19 594 103.9 611 613 611 594 19
20 595 101.6 594 611 613 611 20
21 591 94.6 595 594 611 613 21
22 589 95.9 591 595 594 611 22
23 584 104.7 589 591 595 594 23
24 573 102.8 584 589 591 595 24
25 567 98.1 573 584 589 591 25
26 569 113.9 567 573 584 589 26
27 621 80.9 569 567 573 584 27
28 629 95.7 621 569 567 573 28
29 628 113.2 629 621 569 567 29
30 612 105.9 628 629 621 569 30
31 595 108.8 612 628 629 621 31
32 597 102.3 595 612 628 629 32
33 593 99.0 597 595 612 628 33
34 590 100.7 593 597 595 612 34
35 580 115.5 590 593 597 595 35
36 574 100.7 580 590 593 597 36
37 573 109.9 574 580 590 593 37
38 573 114.6 573 574 580 590 38
39 620 85.4 573 573 574 580 39
40 626 100.5 620 573 573 574 40
41 620 114.8 626 620 573 573 41
42 588 116.5 620 626 620 573 42
43 566 112.9 588 620 626 620 43
44 557 102.0 566 588 620 626 44
45 561 106.0 557 566 588 620 45
46 549 105.3 561 557 566 588 46
47 532 118.8 549 561 557 566 47
48 526 106.1 532 549 561 557 48
49 511 109.3 526 532 549 561 49
50 499 117.2 511 526 532 549 50
51 555 92.5 499 511 526 532 51
52 565 104.2 555 499 511 526 52
53 542 112.5 565 555 499 511 53
54 527 122.4 542 565 555 499 54
55 510 113.3 527 542 565 555 55
56 514 100.0 510 527 542 565 56
57 517 110.7 514 510 527 542 57
58 508 112.8 517 514 510 527 58
59 493 109.8 508 517 514 510 59
60 490 117.3 493 508 517 514 60
61 469 109.1 490 493 508 517 61
62 478 115.9 469 490 493 508 62
63 528 96.0 478 469 490 493 63
64 534 99.8 528 478 469 490 64
65 518 116.8 534 528 478 469 65
66 506 115.7 518 534 528 478 66
67 502 99.4 506 518 534 528 67
68 516 94.3 502 506 518 534 68
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X `Yt-1` `Yt-2` `Yt-3` `Yt-4`
233.322207 -1.477164 0.958869 0.018595 -0.111140 -0.008305
t
0.023881
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-24.163 -9.475 -1.839 8.872 31.998
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 233.322207 37.136966 6.283 3.89e-08 ***
X -1.477164 0.245121 -6.026 1.06e-07 ***
`Yt-1` 0.958869 0.105429 9.095 5.90e-13 ***
`Yt-2` 0.018595 0.174982 0.106 0.916
`Yt-3` -0.111140 0.156051 -0.712 0.479
`Yt-4` -0.008305 0.100585 -0.083 0.934
t 0.023881 0.103753 0.230 0.819
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.65 on 61 degrees of freedom
Multiple R-squared: 0.8977, Adjusted R-squared: 0.8876
F-statistic: 89.18 on 6 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.6964421 0.6071158 0.30355788
[2,] 0.5379257 0.9241485 0.46207427
[3,] 0.8509633 0.2980734 0.14903672
[4,] 0.9409025 0.1181950 0.05909750
[5,] 0.9312842 0.1374315 0.06871575
[6,] 0.9107688 0.1784625 0.08923123
[7,] 0.8636714 0.2726571 0.13632857
[8,] 0.8199718 0.3600563 0.18002816
[9,] 0.7526560 0.4946881 0.24734404
[10,] 0.7054442 0.5891116 0.29455580
[11,] 0.6403235 0.7193531 0.35967655
[12,] 0.6235878 0.7528245 0.37641225
[13,] 0.6250233 0.7499535 0.37497674
[14,] 0.5745185 0.8509630 0.42548152
[15,] 0.6733667 0.6532666 0.32663329
[16,] 0.8561603 0.2876793 0.14383965
[17,] 0.8178932 0.3642135 0.18210677
[18,] 0.7641166 0.4717667 0.23588335
[19,] 0.7612319 0.4775363 0.23876814
[20,] 0.7051118 0.5897763 0.29488816
[21,] 0.7664539 0.4670921 0.23354606
[22,] 0.7141215 0.5717569 0.28587846
[23,] 0.6496140 0.7007719 0.35038597
[24,] 0.6247953 0.7504094 0.37520468
[25,] 0.6350399 0.7299202 0.36496011
[26,] 0.5635675 0.8728650 0.43643251
[27,] 0.7179240 0.5641521 0.28207603
[28,] 0.6518977 0.6962046 0.34810230
[29,] 0.6007074 0.7985851 0.39929256
[30,] 0.5234014 0.9531972 0.47659860
[31,] 0.4611557 0.9223114 0.53884429
[32,] 0.5651754 0.8696492 0.43482458
[33,] 0.5518191 0.8963617 0.44818087
[34,] 0.5322111 0.9355779 0.46778893
[35,] 0.4786204 0.9572409 0.52137956
[36,] 0.6015398 0.7969205 0.39846024
[37,] 0.6114985 0.7770029 0.38850147
[38,] 0.6402109 0.7195782 0.35978912
[39,] 0.6228977 0.7542045 0.37710226
[40,] 0.6474515 0.7050969 0.35254846
[41,] 0.5612288 0.8775423 0.43877117
[42,] 0.5423626 0.9152748 0.45763740
[43,] 0.4559892 0.9119785 0.54401077
[44,] 0.4583541 0.9167083 0.54164586
[45,] 0.3660744 0.7321487 0.63392564
[46,] 0.2739852 0.5479705 0.72601477
[47,] 0.1890129 0.3780259 0.81098706
[48,] 0.2352799 0.4705599 0.76472006
[49,] 0.3096695 0.6193390 0.69033052
> postscript(file="/var/www/html/rcomp/tmp/16cw61258723095.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/2e2jh1258723095.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/30rfc1258723095.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/4clx41258723096.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/5rkwq1258723096.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 = 68
Frequency = 1
1 2 3 4 5 6
-23.21432373 -0.88939622 30.05348441 -15.73638096 0.53527317 -1.02828534
7 8 9 10 11 12
-24.16306777 -1.32767330 -5.80813253 -12.99473075 5.74119839 -13.14265738
13 14 15 16 17 18
-17.68481380 15.21308108 23.84184074 4.82005898 15.71071602 11.55773775
19 20 21 22 23 24
-10.72650341 3.55357356 -11.65888638 -9.85156193 0.08567064 -9.34954074
25 26 27 28 29 30
-11.93105643 18.76969911 18.92917139 -1.88927618 14.47178842 -5.72938455
31 32 33 34 35 36
-1.78800132 7.14014543 -5.14655302 -3.88322860 10.98699698 -7.68238259
37 38 39 40 41 42
10.45616633 17.30907721 20.42071311 3.47417630 11.93827325 -6.70919859
43 44 45 46 47 48
-2.19830340 -6.25012817 9.06724025 -10.36962481 2.79728973 -5.09284429
49 50 51 52 53 54
-10.62093664 1.53034424 31.99783894 4.06629479 -18.78537197 8.80690100
55 56 57 58 59 60
-5.27197133 -6.83560447 6.56868180 -4.31809794 -14.89605606 8.07582480
61 62 63 64 65 66
-22.88061865 14.59040974 26.47363205 -12.40670196 -9.17588728 -1.96255859
67 68
-17.17817197 -8.40538658
> postscript(file="/var/www/html/rcomp/tmp/64ej71258723096.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 = 68
Frequency = 1
lag(myerror, k = 1) myerror
0 -23.21432373 NA
1 -0.88939622 -23.21432373
2 30.05348441 -0.88939622
3 -15.73638096 30.05348441
4 0.53527317 -15.73638096
5 -1.02828534 0.53527317
6 -24.16306777 -1.02828534
7 -1.32767330 -24.16306777
8 -5.80813253 -1.32767330
9 -12.99473075 -5.80813253
10 5.74119839 -12.99473075
11 -13.14265738 5.74119839
12 -17.68481380 -13.14265738
13 15.21308108 -17.68481380
14 23.84184074 15.21308108
15 4.82005898 23.84184074
16 15.71071602 4.82005898
17 11.55773775 15.71071602
18 -10.72650341 11.55773775
19 3.55357356 -10.72650341
20 -11.65888638 3.55357356
21 -9.85156193 -11.65888638
22 0.08567064 -9.85156193
23 -9.34954074 0.08567064
24 -11.93105643 -9.34954074
25 18.76969911 -11.93105643
26 18.92917139 18.76969911
27 -1.88927618 18.92917139
28 14.47178842 -1.88927618
29 -5.72938455 14.47178842
30 -1.78800132 -5.72938455
31 7.14014543 -1.78800132
32 -5.14655302 7.14014543
33 -3.88322860 -5.14655302
34 10.98699698 -3.88322860
35 -7.68238259 10.98699698
36 10.45616633 -7.68238259
37 17.30907721 10.45616633
38 20.42071311 17.30907721
39 3.47417630 20.42071311
40 11.93827325 3.47417630
41 -6.70919859 11.93827325
42 -2.19830340 -6.70919859
43 -6.25012817 -2.19830340
44 9.06724025 -6.25012817
45 -10.36962481 9.06724025
46 2.79728973 -10.36962481
47 -5.09284429 2.79728973
48 -10.62093664 -5.09284429
49 1.53034424 -10.62093664
50 31.99783894 1.53034424
51 4.06629479 31.99783894
52 -18.78537197 4.06629479
53 8.80690100 -18.78537197
54 -5.27197133 8.80690100
55 -6.83560447 -5.27197133
56 6.56868180 -6.83560447
57 -4.31809794 6.56868180
58 -14.89605606 -4.31809794
59 8.07582480 -14.89605606
60 -22.88061865 8.07582480
61 14.59040974 -22.88061865
62 26.47363205 14.59040974
63 -12.40670196 26.47363205
64 -9.17588728 -12.40670196
65 -1.96255859 -9.17588728
66 -17.17817197 -1.96255859
67 -8.40538658 -17.17817197
68 NA -8.40538658
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.88939622 -23.21432373
[2,] 30.05348441 -0.88939622
[3,] -15.73638096 30.05348441
[4,] 0.53527317 -15.73638096
[5,] -1.02828534 0.53527317
[6,] -24.16306777 -1.02828534
[7,] -1.32767330 -24.16306777
[8,] -5.80813253 -1.32767330
[9,] -12.99473075 -5.80813253
[10,] 5.74119839 -12.99473075
[11,] -13.14265738 5.74119839
[12,] -17.68481380 -13.14265738
[13,] 15.21308108 -17.68481380
[14,] 23.84184074 15.21308108
[15,] 4.82005898 23.84184074
[16,] 15.71071602 4.82005898
[17,] 11.55773775 15.71071602
[18,] -10.72650341 11.55773775
[19,] 3.55357356 -10.72650341
[20,] -11.65888638 3.55357356
[21,] -9.85156193 -11.65888638
[22,] 0.08567064 -9.85156193
[23,] -9.34954074 0.08567064
[24,] -11.93105643 -9.34954074
[25,] 18.76969911 -11.93105643
[26,] 18.92917139 18.76969911
[27,] -1.88927618 18.92917139
[28,] 14.47178842 -1.88927618
[29,] -5.72938455 14.47178842
[30,] -1.78800132 -5.72938455
[31,] 7.14014543 -1.78800132
[32,] -5.14655302 7.14014543
[33,] -3.88322860 -5.14655302
[34,] 10.98699698 -3.88322860
[35,] -7.68238259 10.98699698
[36,] 10.45616633 -7.68238259
[37,] 17.30907721 10.45616633
[38,] 20.42071311 17.30907721
[39,] 3.47417630 20.42071311
[40,] 11.93827325 3.47417630
[41,] -6.70919859 11.93827325
[42,] -2.19830340 -6.70919859
[43,] -6.25012817 -2.19830340
[44,] 9.06724025 -6.25012817
[45,] -10.36962481 9.06724025
[46,] 2.79728973 -10.36962481
[47,] -5.09284429 2.79728973
[48,] -10.62093664 -5.09284429
[49,] 1.53034424 -10.62093664
[50,] 31.99783894 1.53034424
[51,] 4.06629479 31.99783894
[52,] -18.78537197 4.06629479
[53,] 8.80690100 -18.78537197
[54,] -5.27197133 8.80690100
[55,] -6.83560447 -5.27197133
[56,] 6.56868180 -6.83560447
[57,] -4.31809794 6.56868180
[58,] -14.89605606 -4.31809794
[59,] 8.07582480 -14.89605606
[60,] -22.88061865 8.07582480
[61,] 14.59040974 -22.88061865
[62,] 26.47363205 14.59040974
[63,] -12.40670196 26.47363205
[64,] -9.17588728 -12.40670196
[65,] -1.96255859 -9.17588728
[66,] -17.17817197 -1.96255859
[67,] -8.40538658 -17.17817197
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.88939622 -23.21432373
2 30.05348441 -0.88939622
3 -15.73638096 30.05348441
4 0.53527317 -15.73638096
5 -1.02828534 0.53527317
6 -24.16306777 -1.02828534
7 -1.32767330 -24.16306777
8 -5.80813253 -1.32767330
9 -12.99473075 -5.80813253
10 5.74119839 -12.99473075
11 -13.14265738 5.74119839
12 -17.68481380 -13.14265738
13 15.21308108 -17.68481380
14 23.84184074 15.21308108
15 4.82005898 23.84184074
16 15.71071602 4.82005898
17 11.55773775 15.71071602
18 -10.72650341 11.55773775
19 3.55357356 -10.72650341
20 -11.65888638 3.55357356
21 -9.85156193 -11.65888638
22 0.08567064 -9.85156193
23 -9.34954074 0.08567064
24 -11.93105643 -9.34954074
25 18.76969911 -11.93105643
26 18.92917139 18.76969911
27 -1.88927618 18.92917139
28 14.47178842 -1.88927618
29 -5.72938455 14.47178842
30 -1.78800132 -5.72938455
31 7.14014543 -1.78800132
32 -5.14655302 7.14014543
33 -3.88322860 -5.14655302
34 10.98699698 -3.88322860
35 -7.68238259 10.98699698
36 10.45616633 -7.68238259
37 17.30907721 10.45616633
38 20.42071311 17.30907721
39 3.47417630 20.42071311
40 11.93827325 3.47417630
41 -6.70919859 11.93827325
42 -2.19830340 -6.70919859
43 -6.25012817 -2.19830340
44 9.06724025 -6.25012817
45 -10.36962481 9.06724025
46 2.79728973 -10.36962481
47 -5.09284429 2.79728973
48 -10.62093664 -5.09284429
49 1.53034424 -10.62093664
50 31.99783894 1.53034424
51 4.06629479 31.99783894
52 -18.78537197 4.06629479
53 8.80690100 -18.78537197
54 -5.27197133 8.80690100
55 -6.83560447 -5.27197133
56 6.56868180 -6.83560447
57 -4.31809794 6.56868180
58 -14.89605606 -4.31809794
59 8.07582480 -14.89605606
60 -22.88061865 8.07582480
61 14.59040974 -22.88061865
62 26.47363205 14.59040974
63 -12.40670196 26.47363205
64 -9.17588728 -12.40670196
65 -1.96255859 -9.17588728
66 -17.17817197 -1.96255859
67 -8.40538658 -17.17817197
> 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/76h7c1258723096.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/84iw11258723096.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/9l4371258723096.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/10awug1258723096.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/112r001258723096.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/12de691258723096.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/13o1yo1258723096.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/14tvf61258723096.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/15ykei1258723096.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/16w8n01258723096.tab")
+ }
>
> system("convert tmp/16cw61258723095.ps tmp/16cw61258723095.png")
> system("convert tmp/2e2jh1258723095.ps tmp/2e2jh1258723095.png")
> system("convert tmp/30rfc1258723095.ps tmp/30rfc1258723095.png")
> system("convert tmp/4clx41258723096.ps tmp/4clx41258723096.png")
> system("convert tmp/5rkwq1258723096.ps tmp/5rkwq1258723096.png")
> system("convert tmp/64ej71258723096.ps tmp/64ej71258723096.png")
> system("convert tmp/76h7c1258723096.ps tmp/76h7c1258723096.png")
> system("convert tmp/84iw11258723096.ps tmp/84iw11258723096.png")
> system("convert tmp/9l4371258723096.ps tmp/9l4371258723096.png")
> system("convert tmp/10awug1258723096.ps tmp/10awug1258723096.png")
>
>
> proc.time()
user system elapsed
2.578 1.571 2.975