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(101.9
+ ,436
+ ,443
+ ,448
+ ,460
+ ,467
+ ,106.2
+ ,431
+ ,436
+ ,443
+ ,448
+ ,460
+ ,81
+ ,484
+ ,431
+ ,436
+ ,443
+ ,448
+ ,94.7
+ ,510
+ ,484
+ ,431
+ ,436
+ ,443
+ ,101
+ ,513
+ ,510
+ ,484
+ ,431
+ ,436
+ ,109.4
+ ,503
+ ,513
+ ,510
+ ,484
+ ,431
+ ,102.3
+ ,471
+ ,503
+ ,513
+ ,510
+ ,484
+ ,90.7
+ ,471
+ ,471
+ ,503
+ ,513
+ ,510
+ ,96.2
+ ,476
+ ,471
+ ,471
+ ,503
+ ,513
+ ,96.1
+ ,475
+ ,476
+ ,471
+ ,471
+ ,503
+ ,106
+ ,470
+ ,475
+ ,476
+ ,471
+ ,471
+ ,103.1
+ ,461
+ ,470
+ ,475
+ ,476
+ ,471
+ ,102
+ ,455
+ ,461
+ ,470
+ ,475
+ ,476
+ ,104.7
+ ,456
+ ,455
+ ,461
+ ,470
+ ,475
+ ,86
+ ,517
+ ,456
+ ,455
+ ,461
+ ,470
+ ,92.1
+ ,525
+ ,517
+ ,456
+ ,455
+ ,461
+ ,106.9
+ ,523
+ ,525
+ ,517
+ ,456
+ ,455
+ ,112.6
+ ,519
+ ,523
+ ,525
+ ,517
+ ,456
+ ,101.7
+ ,509
+ ,519
+ ,523
+ ,525
+ ,517
+ ,92
+ ,512
+ ,509
+ ,519
+ ,523
+ ,525
+ ,97.4
+ ,519
+ ,512
+ ,509
+ ,519
+ ,523
+ ,97
+ ,517
+ ,519
+ ,512
+ ,509
+ ,519
+ ,105.4
+ ,510
+ ,517
+ ,519
+ ,512
+ ,509
+ ,102.7
+ ,509
+ ,510
+ ,517
+ ,519
+ ,512
+ ,98.1
+ ,501
+ ,509
+ ,510
+ ,517
+ ,519
+ ,104.5
+ ,507
+ ,501
+ ,509
+ ,510
+ ,517
+ ,87.4
+ ,569
+ ,507
+ ,501
+ ,509
+ ,510
+ ,89.9
+ ,580
+ ,569
+ ,507
+ ,501
+ ,509
+ ,109.8
+ ,578
+ ,580
+ ,569
+ ,507
+ ,501
+ ,111.7
+ ,565
+ ,578
+ ,580
+ ,569
+ ,507
+ ,98.6
+ ,547
+ ,565
+ ,578
+ ,580
+ ,569
+ ,96.9
+ ,555
+ ,547
+ ,565
+ ,578
+ ,580
+ ,95.1
+ ,562
+ ,555
+ ,547
+ ,565
+ ,578
+ ,97
+ ,561
+ ,562
+ ,555
+ ,547
+ ,565
+ ,112.7
+ ,555
+ ,561
+ ,562
+ ,555
+ ,547
+ ,102.9
+ ,544
+ ,555
+ ,561
+ ,562
+ ,555
+ ,97.4
+ ,537
+ ,544
+ ,555
+ ,561
+ ,562
+ ,111.4
+ ,543
+ ,537
+ ,544
+ ,555
+ ,561
+ ,87.4
+ ,594
+ ,543
+ ,537
+ ,544
+ ,555
+ ,96.8
+ ,611
+ ,594
+ ,543
+ ,537
+ ,544
+ ,114.1
+ ,613
+ ,611
+ ,594
+ ,543
+ ,537
+ ,110.3
+ ,611
+ ,613
+ ,611
+ ,594
+ ,543
+ ,103.9
+ ,594
+ ,611
+ ,613
+ ,611
+ ,594
+ ,101.6
+ ,595
+ ,594
+ ,611
+ ,613
+ ,611
+ ,94.6
+ ,591
+ ,595
+ ,594
+ ,611
+ ,613
+ ,95.9
+ ,589
+ ,591
+ ,595
+ ,594
+ ,611
+ ,104.7
+ ,584
+ ,589
+ ,591
+ ,595
+ ,594
+ ,102.8
+ ,573
+ ,584
+ ,589
+ ,591
+ ,595
+ ,98.1
+ ,567
+ ,573
+ ,584
+ ,589
+ ,591
+ ,113.9
+ ,569
+ ,567
+ ,573
+ ,584
+ ,589
+ ,80.9
+ ,621
+ ,569
+ ,567
+ ,573
+ ,584
+ ,95.7
+ ,629
+ ,621
+ ,569
+ ,567
+ ,573
+ ,113.2
+ ,628
+ ,629
+ ,621
+ ,569
+ ,567
+ ,105.9
+ ,612
+ ,628
+ ,629
+ ,621
+ ,569
+ ,108.8
+ ,595
+ ,612
+ ,628
+ ,629
+ ,621
+ ,102.3
+ ,597
+ ,595
+ ,612
+ ,628
+ ,629
+ ,99
+ ,593
+ ,597
+ ,595
+ ,612
+ ,628
+ ,100.7
+ ,590
+ ,593
+ ,597
+ ,595
+ ,612
+ ,115.5
+ ,580
+ ,590
+ ,593
+ ,597
+ ,595
+ ,100.7
+ ,574
+ ,580
+ ,590
+ ,593
+ ,597
+ ,109.9
+ ,573
+ ,574
+ ,580
+ ,590
+ ,593
+ ,114.6
+ ,573
+ ,573
+ ,574
+ ,580
+ ,590
+ ,85.4
+ ,620
+ ,573
+ ,573
+ ,574
+ ,580
+ ,100.5
+ ,626
+ ,620
+ ,573
+ ,573
+ ,574
+ ,114.8
+ ,620
+ ,626
+ ,620
+ ,573
+ ,573
+ ,116.5
+ ,588
+ ,620
+ ,626
+ ,620
+ ,573
+ ,112.9
+ ,566
+ ,588
+ ,620
+ ,626
+ ,620
+ ,102
+ ,557
+ ,566
+ ,588
+ ,620
+ ,626)
+ ,dim=c(6
+ ,68)
+ ,dimnames=list(c('X'
+ ,'Y'
+ ,'Y(t-1)'
+ ,'Y(t-2)'
+ ,'Y(t-3)'
+ ,'Y(t-4)')
+ ,1:68))
> y <- array(NA,dim=c(6,68),dimnames=list(c('X','Y','Y(t-1)','Y(t-2)','Y(t-3)','Y(t-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 = '2'
> #'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 Y(t-1) Y(t-2) Y(t-3) Y(t-4) t
1 436 101.9 443 448 460 467 1
2 431 106.2 436 443 448 460 2
3 484 81.0 431 436 443 448 3
4 510 94.7 484 431 436 443 4
5 513 101.0 510 484 431 436 5
6 503 109.4 513 510 484 431 6
7 471 102.3 503 513 510 484 7
8 471 90.7 471 503 513 510 8
9 476 96.2 471 471 503 513 9
10 475 96.1 476 471 471 503 10
11 470 106.0 475 476 471 471 11
12 461 103.1 470 475 476 471 12
13 455 102.0 461 470 475 476 13
14 456 104.7 455 461 470 475 14
15 517 86.0 456 455 461 470 15
16 525 92.1 517 456 455 461 16
17 523 106.9 525 517 456 455 17
18 519 112.6 523 525 517 456 18
19 509 101.7 519 523 525 517 19
20 512 92.0 509 519 523 525 20
21 519 97.4 512 509 519 523 21
22 517 97.0 519 512 509 519 22
23 510 105.4 517 519 512 509 23
24 509 102.7 510 517 519 512 24
25 501 98.1 509 510 517 519 25
26 507 104.5 501 509 510 517 26
27 569 87.4 507 501 509 510 27
28 580 89.9 569 507 501 509 28
29 578 109.8 580 569 507 501 29
30 565 111.7 578 580 569 507 30
31 547 98.6 565 578 580 569 31
32 555 96.9 547 565 578 580 32
33 562 95.1 555 547 565 578 33
34 561 97.0 562 555 547 565 34
35 555 112.7 561 562 555 547 35
36 544 102.9 555 561 562 555 36
37 537 97.4 544 555 561 562 37
38 543 111.4 537 544 555 561 38
39 594 87.4 543 537 544 555 39
40 611 96.8 594 543 537 544 40
41 613 114.1 611 594 543 537 41
42 611 110.3 613 611 594 543 42
43 594 103.9 611 613 611 594 43
44 595 101.6 594 611 613 611 44
45 591 94.6 595 594 611 613 45
46 589 95.9 591 595 594 611 46
47 584 104.7 589 591 595 594 47
48 573 102.8 584 589 591 595 48
49 567 98.1 573 584 589 591 49
50 569 113.9 567 573 584 589 50
51 621 80.9 569 567 573 584 51
52 629 95.7 621 569 567 573 52
53 628 113.2 629 621 569 567 53
54 612 105.9 628 629 621 569 54
55 595 108.8 612 628 629 621 55
56 597 102.3 595 612 628 629 56
57 593 99.0 597 595 612 628 57
58 590 100.7 593 597 595 612 58
59 580 115.5 590 593 597 595 59
60 574 100.7 580 590 593 597 60
61 573 109.9 574 580 590 593 61
62 573 114.6 573 574 580 590 62
63 620 85.4 573 573 574 580 63
64 626 100.5 620 573 573 574 64
65 620 114.8 626 620 573 573 65
66 588 116.5 620 626 620 573 66
67 566 112.9 588 620 626 620 67
68 557 102.0 566 588 620 626 68
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X `Y(t-1)` `Y(t-2)` `Y(t-3)` `Y(t-4)`
335.5712 -1.7048 0.8989 0.0937 -0.2080 -0.1422
t
0.9694
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-22.243 -8.849 -1.828 7.604 31.987
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 335.57116 45.26014 7.414 4.49e-10 ***
X -1.70478 0.22764 -7.489 3.34e-10 ***
`Y(t-1)` 0.89890 0.09667 9.298 2.67e-13 ***
`Y(t-2)` 0.09369 0.15148 0.619 0.539
`Y(t-3)` -0.20801 0.13811 -1.506 0.137
`Y(t-4)` -0.14224 0.10325 -1.378 0.173
t 0.96939 0.21700 4.467 3.49e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12.07 on 61 degrees of freedom
Multiple R-squared: 0.9514, Adjusted R-squared: 0.9466
F-statistic: 199 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.23651048 0.4730210 0.76348952
[2,] 0.16517809 0.3303562 0.83482191
[3,] 0.20342324 0.4068465 0.79657676
[4,] 0.15188746 0.3037749 0.84811254
[5,] 0.09171503 0.1834301 0.90828497
[6,] 0.16150315 0.3230063 0.83849685
[7,] 0.39497391 0.7899478 0.60502609
[8,] 0.34599326 0.6919865 0.65400674
[9,] 0.39060424 0.7812085 0.60939576
[10,] 0.34077989 0.6815598 0.65922011
[11,] 0.27748296 0.5549659 0.72251704
[12,] 0.25651363 0.5130273 0.74348637
[13,] 0.21630095 0.4326019 0.78369905
[14,] 0.17558801 0.3511760 0.82441199
[15,] 0.13374855 0.2674971 0.86625145
[16,] 0.22059628 0.4411926 0.77940372
[17,] 0.21948816 0.4389763 0.78051184
[18,] 0.46439142 0.9287828 0.53560858
[19,] 0.50022224 0.9995555 0.49977776
[20,] 0.55696860 0.8860628 0.44303140
[21,] 0.49483251 0.9896650 0.50516749
[22,] 0.51042310 0.9791538 0.48957690
[23,] 0.54146691 0.9170662 0.45853309
[24,] 0.50019541 0.9996092 0.49980459
[25,] 0.51988757 0.9602249 0.48011243
[26,] 0.51692196 0.9661561 0.48307804
[27,] 0.60709849 0.7858030 0.39290151
[28,] 0.90394436 0.1921113 0.09605564
[29,] 0.90786055 0.1842789 0.09213945
[30,] 0.89529355 0.2094129 0.10470645
[31,] 0.86432213 0.2713557 0.13567787
[32,] 0.85387293 0.2922541 0.14612707
[33,] 0.83312581 0.3337484 0.16687419
[34,] 0.78062909 0.4387418 0.21937091
[35,] 0.77674483 0.4465103 0.22325517
[36,] 0.71729505 0.5654099 0.28270495
[37,] 0.69086066 0.6182787 0.30913934
[38,] 0.61008244 0.7798351 0.38991756
[39,] 0.68296152 0.6340770 0.31703848
[40,] 0.94111545 0.1177691 0.05888455
[41,] 0.90924253 0.1815149 0.09075747
[42,] 0.86706915 0.2658617 0.13293085
[43,] 0.84758593 0.3048281 0.15241407
[44,] 0.76961931 0.4607614 0.23038069
[45,] 0.72170326 0.5565935 0.27829674
[46,] 0.60894792 0.7821042 0.39105208
[47,] 0.75100887 0.4979823 0.24899113
[48,] 0.72053631 0.5589274 0.27946369
[49,] 0.57499995 0.8500001 0.42500005
> postscript(file="/var/www/html/rcomp/tmp/1ul8n1260877241.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/23aiz1260877241.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/3x09i1260877241.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/4uizq1260877241.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/5pa441260877241.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
-4.9022437 -0.2720797 11.2015660 10.2470595 -7.3552323 1.1757961
7 8 9 10 11 12
-22.2425425 -8.9632091 5.7885821 -8.9244015 -2.1378659 -11.4228847
13 14 15 16 17 18
-11.2057340 -1.5178376 23.7134388 -16.3117230 -7.6024658 11.0241815
19 20 21 22 23 24
-4.4033926 -8.8233932 3.5367447 -9.3369818 -2.6427275 -0.8525520
25 26 27 28 29 30
-15.5294618 5.9561004 31.9874578 -11.8203398 3.5484210 7.3350536
31 32 33 34 35 36
-10.9866575 11.6927543 6.1615002 -5.2039467 13.9383577 -6.6567632
37 38 39 40 41 42
-12.7646470 22.0655164 23.3024075 5.9310587 16.6468764 15.2704464
43 44 45 46 47 48
-1.2085999 13.2038834 -3.1365638 -4.2084135 4.7866794 -6.4296717
49 50 51 52 53 54
-12.0401122 21.0255103 11.5635620 -5.9181074 9.4453065 -8.7188544
55 56 57 58 59 60
1.7924273 9.4523667 -4.8181114 -8.2931538 7.0375086 -16.4400215
61 62 63 64 65 66
2.4119010 8.4092463 2.0835890 -10.4535063 -2.9838954 -18.4476788
67 68
-10.2937787 -16.4657479
> postscript(file="/var/www/html/rcomp/tmp/6s28i1260877241.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 -4.9022437 NA
1 -0.2720797 -4.9022437
2 11.2015660 -0.2720797
3 10.2470595 11.2015660
4 -7.3552323 10.2470595
5 1.1757961 -7.3552323
6 -22.2425425 1.1757961
7 -8.9632091 -22.2425425
8 5.7885821 -8.9632091
9 -8.9244015 5.7885821
10 -2.1378659 -8.9244015
11 -11.4228847 -2.1378659
12 -11.2057340 -11.4228847
13 -1.5178376 -11.2057340
14 23.7134388 -1.5178376
15 -16.3117230 23.7134388
16 -7.6024658 -16.3117230
17 11.0241815 -7.6024658
18 -4.4033926 11.0241815
19 -8.8233932 -4.4033926
20 3.5367447 -8.8233932
21 -9.3369818 3.5367447
22 -2.6427275 -9.3369818
23 -0.8525520 -2.6427275
24 -15.5294618 -0.8525520
25 5.9561004 -15.5294618
26 31.9874578 5.9561004
27 -11.8203398 31.9874578
28 3.5484210 -11.8203398
29 7.3350536 3.5484210
30 -10.9866575 7.3350536
31 11.6927543 -10.9866575
32 6.1615002 11.6927543
33 -5.2039467 6.1615002
34 13.9383577 -5.2039467
35 -6.6567632 13.9383577
36 -12.7646470 -6.6567632
37 22.0655164 -12.7646470
38 23.3024075 22.0655164
39 5.9310587 23.3024075
40 16.6468764 5.9310587
41 15.2704464 16.6468764
42 -1.2085999 15.2704464
43 13.2038834 -1.2085999
44 -3.1365638 13.2038834
45 -4.2084135 -3.1365638
46 4.7866794 -4.2084135
47 -6.4296717 4.7866794
48 -12.0401122 -6.4296717
49 21.0255103 -12.0401122
50 11.5635620 21.0255103
51 -5.9181074 11.5635620
52 9.4453065 -5.9181074
53 -8.7188544 9.4453065
54 1.7924273 -8.7188544
55 9.4523667 1.7924273
56 -4.8181114 9.4523667
57 -8.2931538 -4.8181114
58 7.0375086 -8.2931538
59 -16.4400215 7.0375086
60 2.4119010 -16.4400215
61 8.4092463 2.4119010
62 2.0835890 8.4092463
63 -10.4535063 2.0835890
64 -2.9838954 -10.4535063
65 -18.4476788 -2.9838954
66 -10.2937787 -18.4476788
67 -16.4657479 -10.2937787
68 NA -16.4657479
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.2720797 -4.9022437
[2,] 11.2015660 -0.2720797
[3,] 10.2470595 11.2015660
[4,] -7.3552323 10.2470595
[5,] 1.1757961 -7.3552323
[6,] -22.2425425 1.1757961
[7,] -8.9632091 -22.2425425
[8,] 5.7885821 -8.9632091
[9,] -8.9244015 5.7885821
[10,] -2.1378659 -8.9244015
[11,] -11.4228847 -2.1378659
[12,] -11.2057340 -11.4228847
[13,] -1.5178376 -11.2057340
[14,] 23.7134388 -1.5178376
[15,] -16.3117230 23.7134388
[16,] -7.6024658 -16.3117230
[17,] 11.0241815 -7.6024658
[18,] -4.4033926 11.0241815
[19,] -8.8233932 -4.4033926
[20,] 3.5367447 -8.8233932
[21,] -9.3369818 3.5367447
[22,] -2.6427275 -9.3369818
[23,] -0.8525520 -2.6427275
[24,] -15.5294618 -0.8525520
[25,] 5.9561004 -15.5294618
[26,] 31.9874578 5.9561004
[27,] -11.8203398 31.9874578
[28,] 3.5484210 -11.8203398
[29,] 7.3350536 3.5484210
[30,] -10.9866575 7.3350536
[31,] 11.6927543 -10.9866575
[32,] 6.1615002 11.6927543
[33,] -5.2039467 6.1615002
[34,] 13.9383577 -5.2039467
[35,] -6.6567632 13.9383577
[36,] -12.7646470 -6.6567632
[37,] 22.0655164 -12.7646470
[38,] 23.3024075 22.0655164
[39,] 5.9310587 23.3024075
[40,] 16.6468764 5.9310587
[41,] 15.2704464 16.6468764
[42,] -1.2085999 15.2704464
[43,] 13.2038834 -1.2085999
[44,] -3.1365638 13.2038834
[45,] -4.2084135 -3.1365638
[46,] 4.7866794 -4.2084135
[47,] -6.4296717 4.7866794
[48,] -12.0401122 -6.4296717
[49,] 21.0255103 -12.0401122
[50,] 11.5635620 21.0255103
[51,] -5.9181074 11.5635620
[52,] 9.4453065 -5.9181074
[53,] -8.7188544 9.4453065
[54,] 1.7924273 -8.7188544
[55,] 9.4523667 1.7924273
[56,] -4.8181114 9.4523667
[57,] -8.2931538 -4.8181114
[58,] 7.0375086 -8.2931538
[59,] -16.4400215 7.0375086
[60,] 2.4119010 -16.4400215
[61,] 8.4092463 2.4119010
[62,] 2.0835890 8.4092463
[63,] -10.4535063 2.0835890
[64,] -2.9838954 -10.4535063
[65,] -18.4476788 -2.9838954
[66,] -10.2937787 -18.4476788
[67,] -16.4657479 -10.2937787
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.2720797 -4.9022437
2 11.2015660 -0.2720797
3 10.2470595 11.2015660
4 -7.3552323 10.2470595
5 1.1757961 -7.3552323
6 -22.2425425 1.1757961
7 -8.9632091 -22.2425425
8 5.7885821 -8.9632091
9 -8.9244015 5.7885821
10 -2.1378659 -8.9244015
11 -11.4228847 -2.1378659
12 -11.2057340 -11.4228847
13 -1.5178376 -11.2057340
14 23.7134388 -1.5178376
15 -16.3117230 23.7134388
16 -7.6024658 -16.3117230
17 11.0241815 -7.6024658
18 -4.4033926 11.0241815
19 -8.8233932 -4.4033926
20 3.5367447 -8.8233932
21 -9.3369818 3.5367447
22 -2.6427275 -9.3369818
23 -0.8525520 -2.6427275
24 -15.5294618 -0.8525520
25 5.9561004 -15.5294618
26 31.9874578 5.9561004
27 -11.8203398 31.9874578
28 3.5484210 -11.8203398
29 7.3350536 3.5484210
30 -10.9866575 7.3350536
31 11.6927543 -10.9866575
32 6.1615002 11.6927543
33 -5.2039467 6.1615002
34 13.9383577 -5.2039467
35 -6.6567632 13.9383577
36 -12.7646470 -6.6567632
37 22.0655164 -12.7646470
38 23.3024075 22.0655164
39 5.9310587 23.3024075
40 16.6468764 5.9310587
41 15.2704464 16.6468764
42 -1.2085999 15.2704464
43 13.2038834 -1.2085999
44 -3.1365638 13.2038834
45 -4.2084135 -3.1365638
46 4.7866794 -4.2084135
47 -6.4296717 4.7866794
48 -12.0401122 -6.4296717
49 21.0255103 -12.0401122
50 11.5635620 21.0255103
51 -5.9181074 11.5635620
52 9.4453065 -5.9181074
53 -8.7188544 9.4453065
54 1.7924273 -8.7188544
55 9.4523667 1.7924273
56 -4.8181114 9.4523667
57 -8.2931538 -4.8181114
58 7.0375086 -8.2931538
59 -16.4400215 7.0375086
60 2.4119010 -16.4400215
61 8.4092463 2.4119010
62 2.0835890 8.4092463
63 -10.4535063 2.0835890
64 -2.9838954 -10.4535063
65 -18.4476788 -2.9838954
66 -10.2937787 -18.4476788
67 -16.4657479 -10.2937787
> 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/7yyct1260877241.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/8qcee1260877241.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/9x6dq1260877241.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/10yvgu1260877241.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/110svv1260877241.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/12xb181260877241.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/137rgp1260877241.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/14mgib1260877241.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/15exrw1260877241.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/164j3e1260877241.tab")
+ }
>
> try(system("convert tmp/1ul8n1260877241.ps tmp/1ul8n1260877241.png",intern=TRUE))
character(0)
> try(system("convert tmp/23aiz1260877241.ps tmp/23aiz1260877241.png",intern=TRUE))
character(0)
> try(system("convert tmp/3x09i1260877241.ps tmp/3x09i1260877241.png",intern=TRUE))
character(0)
> try(system("convert tmp/4uizq1260877241.ps tmp/4uizq1260877241.png",intern=TRUE))
character(0)
> try(system("convert tmp/5pa441260877241.ps tmp/5pa441260877241.png",intern=TRUE))
character(0)
> try(system("convert tmp/6s28i1260877241.ps tmp/6s28i1260877241.png",intern=TRUE))
character(0)
> try(system("convert tmp/7yyct1260877241.ps tmp/7yyct1260877241.png",intern=TRUE))
character(0)
> try(system("convert tmp/8qcee1260877241.ps tmp/8qcee1260877241.png",intern=TRUE))
character(0)
> try(system("convert tmp/9x6dq1260877241.ps tmp/9x6dq1260877241.png",intern=TRUE))
character(0)
> try(system("convert tmp/10yvgu1260877241.ps tmp/10yvgu1260877241.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.579 1.571 5.650