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(431
+ ,106.2
+ ,436
+ ,443
+ ,448
+ ,460
+ ,467
+ ,484
+ ,81
+ ,431
+ ,436
+ ,443
+ ,448
+ ,460
+ ,510
+ ,94.7
+ ,484
+ ,431
+ ,436
+ ,443
+ ,448
+ ,513
+ ,101
+ ,510
+ ,484
+ ,431
+ ,436
+ ,443
+ ,503
+ ,109.4
+ ,513
+ ,510
+ ,484
+ ,431
+ ,436
+ ,471
+ ,102.3
+ ,503
+ ,513
+ ,510
+ ,484
+ ,431
+ ,471
+ ,90.7
+ ,471
+ ,503
+ ,513
+ ,510
+ ,484
+ ,476
+ ,96.2
+ ,471
+ ,471
+ ,503
+ ,513
+ ,510
+ ,475
+ ,96.1
+ ,476
+ ,471
+ ,471
+ ,503
+ ,513
+ ,470
+ ,106
+ ,475
+ ,476
+ ,471
+ ,471
+ ,503
+ ,461
+ ,103.1
+ ,470
+ ,475
+ ,476
+ ,471
+ ,471
+ ,455
+ ,102
+ ,461
+ ,470
+ ,475
+ ,476
+ ,471
+ ,456
+ ,104.7
+ ,455
+ ,461
+ ,470
+ ,475
+ ,476
+ ,517
+ ,86
+ ,456
+ ,455
+ ,461
+ ,470
+ ,475
+ ,525
+ ,92.1
+ ,517
+ ,456
+ ,455
+ ,461
+ ,470
+ ,523
+ ,106.9
+ ,525
+ ,517
+ ,456
+ ,455
+ ,461
+ ,519
+ ,112.6
+ ,523
+ ,525
+ ,517
+ ,456
+ ,455
+ ,509
+ ,101.7
+ ,519
+ ,523
+ ,525
+ ,517
+ ,456
+ ,512
+ ,92
+ ,509
+ ,519
+ ,523
+ ,525
+ ,517
+ ,519
+ ,97.4
+ ,512
+ ,509
+ ,519
+ ,523
+ ,525
+ ,517
+ ,97
+ ,519
+ ,512
+ ,509
+ ,519
+ ,523
+ ,510
+ ,105.4
+ ,517
+ ,519
+ ,512
+ ,509
+ ,519
+ ,509
+ ,102.7
+ ,510
+ ,517
+ ,519
+ ,512
+ ,509
+ ,501
+ ,98.1
+ ,509
+ ,510
+ ,517
+ ,519
+ ,512
+ ,507
+ ,104.5
+ ,501
+ ,509
+ ,510
+ ,517
+ ,519
+ ,569
+ ,87.4
+ ,507
+ ,501
+ ,509
+ ,510
+ ,517
+ ,580
+ ,89.9
+ ,569
+ ,507
+ ,501
+ ,509
+ ,510
+ ,578
+ ,109.8
+ ,580
+ ,569
+ ,507
+ ,501
+ ,509
+ ,565
+ ,111.7
+ ,578
+ ,580
+ ,569
+ ,507
+ ,501
+ ,547
+ ,98.6
+ ,565
+ ,578
+ ,580
+ ,569
+ ,507
+ ,555
+ ,96.9
+ ,547
+ ,565
+ ,578
+ ,580
+ ,569
+ ,562
+ ,95.1
+ ,555
+ ,547
+ ,565
+ ,578
+ ,580
+ ,561
+ ,97
+ ,562
+ ,555
+ ,547
+ ,565
+ ,578
+ ,555
+ ,112.7
+ ,561
+ ,562
+ ,555
+ ,547
+ ,565
+ ,544
+ ,102.9
+ ,555
+ ,561
+ ,562
+ ,555
+ ,547
+ ,537
+ ,97.4
+ ,544
+ ,555
+ ,561
+ ,562
+ ,555
+ ,543
+ ,111.4
+ ,537
+ ,544
+ ,555
+ ,561
+ ,562
+ ,594
+ ,87.4
+ ,543
+ ,537
+ ,544
+ ,555
+ ,561
+ ,611
+ ,96.8
+ ,594
+ ,543
+ ,537
+ ,544
+ ,555
+ ,613
+ ,114.1
+ ,611
+ ,594
+ ,543
+ ,537
+ ,544
+ ,611
+ ,110.3
+ ,613
+ ,611
+ ,594
+ ,543
+ ,537
+ ,594
+ ,103.9
+ ,611
+ ,613
+ ,611
+ ,594
+ ,543
+ ,595
+ ,101.6
+ ,594
+ ,611
+ ,613
+ ,611
+ ,594
+ ,591
+ ,94.6
+ ,595
+ ,594
+ ,611
+ ,613
+ ,611
+ ,589
+ ,95.9
+ ,591
+ ,595
+ ,594
+ ,611
+ ,613
+ ,584
+ ,104.7
+ ,589
+ ,591
+ ,595
+ ,594
+ ,611
+ ,573
+ ,102.8
+ ,584
+ ,589
+ ,591
+ ,595
+ ,594
+ ,567
+ ,98.1
+ ,573
+ ,584
+ ,589
+ ,591
+ ,595
+ ,569
+ ,113.9
+ ,567
+ ,573
+ ,584
+ ,589
+ ,591
+ ,621
+ ,80.9
+ ,569
+ ,567
+ ,573
+ ,584
+ ,589
+ ,629
+ ,95.7
+ ,621
+ ,569
+ ,567
+ ,573
+ ,584
+ ,628
+ ,113.2
+ ,629
+ ,621
+ ,569
+ ,567
+ ,573
+ ,612
+ ,105.9
+ ,628
+ ,629
+ ,621
+ ,569
+ ,567
+ ,595
+ ,108.8
+ ,612
+ ,628
+ ,629
+ ,621
+ ,569
+ ,597
+ ,102.3
+ ,595
+ ,612
+ ,628
+ ,629
+ ,621
+ ,593
+ ,99
+ ,597
+ ,595
+ ,612
+ ,628
+ ,629
+ ,590
+ ,100.7
+ ,593
+ ,597
+ ,595
+ ,612
+ ,628
+ ,580
+ ,115.5
+ ,590
+ ,593
+ ,597
+ ,595
+ ,612
+ ,574
+ ,100.7
+ ,580
+ ,590
+ ,593
+ ,597
+ ,595
+ ,573
+ ,109.9
+ ,574
+ ,580
+ ,590
+ ,593
+ ,597
+ ,573
+ ,114.6
+ ,573
+ ,574
+ ,580
+ ,590
+ ,593
+ ,620
+ ,85.4
+ ,573
+ ,573
+ ,574
+ ,580
+ ,590
+ ,626
+ ,100.5
+ ,620
+ ,573
+ ,573
+ ,574
+ ,580
+ ,620
+ ,114.8
+ ,626
+ ,620
+ ,573
+ ,573
+ ,574
+ ,588
+ ,116.5
+ ,620
+ ,626
+ ,620
+ ,573
+ ,573
+ ,566
+ ,112.9
+ ,588
+ ,620
+ ,626
+ ,620
+ ,573
+ ,557
+ ,102
+ ,566
+ ,588
+ ,620
+ ,626
+ ,620)
+ ,dim=c(7
+ ,67)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y(t-1)'
+ ,'Y(t-2)'
+ ,'Y(t-3)'
+ ,'Y(t-4)'
+ ,'Y(t-5)
')
+ ,1:67))
> y <- array(NA,dim=c(7,67),dimnames=list(c('Y','X','Y(t-1)','Y(t-2)','Y(t-3)','Y(t-4)','Y(t-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 Y(t-1) Y(t-2) Y(t-3) Y(t-4) Y(t-5)\r t
1 431 106.2 436 443 448 460 467 1
2 484 81.0 431 436 443 448 460 2
3 510 94.7 484 431 436 443 448 3
4 513 101.0 510 484 431 436 443 4
5 503 109.4 513 510 484 431 436 5
6 471 102.3 503 513 510 484 431 6
7 471 90.7 471 503 513 510 484 7
8 476 96.2 471 471 503 513 510 8
9 475 96.1 476 471 471 503 513 9
10 470 106.0 475 476 471 471 503 10
11 461 103.1 470 475 476 471 471 11
12 455 102.0 461 470 475 476 471 12
13 456 104.7 455 461 470 475 476 13
14 517 86.0 456 455 461 470 475 14
15 525 92.1 517 456 455 461 470 15
16 523 106.9 525 517 456 455 461 16
17 519 112.6 523 525 517 456 455 17
18 509 101.7 519 523 525 517 456 18
19 512 92.0 509 519 523 525 517 19
20 519 97.4 512 509 519 523 525 20
21 517 97.0 519 512 509 519 523 21
22 510 105.4 517 519 512 509 519 22
23 509 102.7 510 517 519 512 509 23
24 501 98.1 509 510 517 519 512 24
25 507 104.5 501 509 510 517 519 25
26 569 87.4 507 501 509 510 517 26
27 580 89.9 569 507 501 509 510 27
28 578 109.8 580 569 507 501 509 28
29 565 111.7 578 580 569 507 501 29
30 547 98.6 565 578 580 569 507 30
31 555 96.9 547 565 578 580 569 31
32 562 95.1 555 547 565 578 580 32
33 561 97.0 562 555 547 565 578 33
34 555 112.7 561 562 555 547 565 34
35 544 102.9 555 561 562 555 547 35
36 537 97.4 544 555 561 562 555 36
37 543 111.4 537 544 555 561 562 37
38 594 87.4 543 537 544 555 561 38
39 611 96.8 594 543 537 544 555 39
40 613 114.1 611 594 543 537 544 40
41 611 110.3 613 611 594 543 537 41
42 594 103.9 611 613 611 594 543 42
43 595 101.6 594 611 613 611 594 43
44 591 94.6 595 594 611 613 611 44
45 589 95.9 591 595 594 611 613 45
46 584 104.7 589 591 595 594 611 46
47 573 102.8 584 589 591 595 594 47
48 567 98.1 573 584 589 591 595 48
49 569 113.9 567 573 584 589 591 49
50 621 80.9 569 567 573 584 589 50
51 629 95.7 621 569 567 573 584 51
52 628 113.2 629 621 569 567 573 52
53 612 105.9 628 629 621 569 567 53
54 595 108.8 612 628 629 621 569 54
55 597 102.3 595 612 628 629 621 55
56 593 99.0 597 595 612 628 629 56
57 590 100.7 593 597 595 612 628 57
58 580 115.5 590 593 597 595 612 58
59 574 100.7 580 590 593 597 595 59
60 573 109.9 574 580 590 593 597 60
61 573 114.6 573 574 580 590 593 61
62 620 85.4 573 573 574 580 590 62
63 626 100.5 620 573 573 574 580 63
64 620 114.8 626 620 573 573 574 64
65 588 116.5 620 626 620 573 573 65
66 566 112.9 588 620 626 620 573 66
67 557 102.0 566 588 620 626 620 67
> 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)`
302.5215 -1.7212 0.8668 0.1345 -0.1257 -0.3806
`Y(t-5)\r` t
0.2269 0.7647
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-19.255 -8.842 -1.233 8.167 31.007
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 302.52150 46.18559 6.550 1.55e-08 ***
X -1.72124 0.22301 -7.718 1.64e-10 ***
`Y(t-1)` 0.86675 0.09457 9.165 6.06e-13 ***
`Y(t-2)` 0.13447 0.14831 0.907 0.36827
`Y(t-3)` -0.12571 0.13799 -0.911 0.36601
`Y(t-4)` -0.38060 0.14179 -2.684 0.00942 **
`Y(t-5)\r` 0.22689 0.09416 2.410 0.01911 *
t 0.76469 0.22622 3.380 0.00129 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.69 on 59 degrees of freedom
Multiple R-squared: 0.9525, Adjusted R-squared: 0.9469
F-statistic: 169.2 on 7 and 59 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.02400894 0.04801787 0.9759911
[2,] 0.01292731 0.02585461 0.9870727
[3,] 0.03034712 0.06069423 0.9696529
[4,] 0.18829084 0.37658168 0.8117092
[5,] 0.44366696 0.88733392 0.5563330
[6,] 0.41447789 0.82895578 0.5855221
[7,] 0.41682836 0.83365672 0.5831716
[8,] 0.47038332 0.94076664 0.5296167
[9,] 0.40314339 0.80628677 0.5968566
[10,] 0.34417659 0.68835319 0.6558234
[11,] 0.29868946 0.59737892 0.7013105
[12,] 0.24454228 0.48908455 0.7554577
[13,] 0.18857628 0.37715256 0.8114237
[14,] 0.27035111 0.54070221 0.7296489
[15,] 0.26915717 0.53831434 0.7308428
[16,] 0.52410342 0.95179315 0.4758966
[17,] 0.52984414 0.94031171 0.4701559
[18,] 0.56562648 0.86874705 0.4343735
[19,] 0.48836376 0.97672753 0.5116362
[20,] 0.46622555 0.93245109 0.5337745
[21,] 0.44740193 0.89480385 0.5525981
[22,] 0.37816301 0.75632602 0.6218370
[23,] 0.41016519 0.82033038 0.5898348
[24,] 0.37733183 0.75466366 0.6226682
[25,] 0.45605288 0.91210577 0.5439471
[26,] 0.81929204 0.36141593 0.1807080
[27,] 0.84926054 0.30147891 0.1507395
[28,] 0.82969623 0.34060753 0.1703038
[29,] 0.78768229 0.42463541 0.2123177
[30,] 0.78182900 0.43634200 0.2181710
[31,] 0.74947977 0.50104045 0.2505202
[32,] 0.69657978 0.60684044 0.3034202
[33,] 0.67812349 0.64375303 0.3218765
[34,] 0.62633519 0.74732962 0.3736648
[35,] 0.60245498 0.79509003 0.3975450
[36,] 0.52118012 0.95763977 0.4788199
[37,] 0.59621193 0.80757615 0.4037881
[38,] 0.90708191 0.18583619 0.0929181
[39,] 0.86558976 0.26882049 0.1344102
[40,] 0.80478335 0.39043329 0.1952166
[41,] 0.78498573 0.43002853 0.2150143
[42,] 0.68773141 0.62453717 0.3122686
[43,] 0.62987646 0.74024708 0.3701235
[44,] 0.49838162 0.99676325 0.5016184
[45,] 0.61660979 0.76678042 0.3833902
[46,] 0.55844551 0.88310897 0.4415545
> postscript(file="/var/www/html/rcomp/tmp/1u14h1260713218.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/2w3nw1260713218.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/3cdmn1260713218.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/4vrsm1260713218.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/5i8t31260713218.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
-1.53176890 8.99592457 12.48629472 -6.25526609 -2.31027114 -14.45702445
7 8 9 10 11 12
-7.85984975 4.13076108 -10.64910192 -9.08935430 -11.48829429 -9.89592190
13 14 15 16 17 18
-0.74612209 24.43463140 -13.88218846 -6.42497008 8.68922773 6.89429523
19 20 21 22 23 24
-9.40804518 1.78720054 -10.46234475 -5.49769419 -1.28283469 -14.42508577
25 26 27 28 29 30
5.66526629 31.00650887 -8.79869773 2.75429703 4.40687390 -1.75073157
31 32 33 34 35 36
9.77598360 3.50824504 -8.88583396 8.40287349 -6.88614441 -13.05318394
37 38 39 40 41 42
21.10280316 22.32967634 6.02814110 16.03405536 12.99208512 5.86225138
43 44 45 46 47 48
12.29256274 -6.44898292 -6.99550302 -1.23260410 -7.92998716 -14.57855060
49 50 51 52 53 54
20.04983536 10.72561051 -5.71131021 9.18295075 -11.69646050 9.87599978
55 56 57 58 59 60
7.93025588 -6.16912109 -11.80932038 3.44996308 -15.60254839 2.66005988
61 62 63 64 65 66
10.16745992 2.39754238 -7.25430694 0.05495899 -19.25479439 -1.03058652
67
-15.31978940
> postscript(file="/var/www/html/rcomp/tmp/6ctnt1260713218.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 -1.53176890 NA
1 8.99592457 -1.53176890
2 12.48629472 8.99592457
3 -6.25526609 12.48629472
4 -2.31027114 -6.25526609
5 -14.45702445 -2.31027114
6 -7.85984975 -14.45702445
7 4.13076108 -7.85984975
8 -10.64910192 4.13076108
9 -9.08935430 -10.64910192
10 -11.48829429 -9.08935430
11 -9.89592190 -11.48829429
12 -0.74612209 -9.89592190
13 24.43463140 -0.74612209
14 -13.88218846 24.43463140
15 -6.42497008 -13.88218846
16 8.68922773 -6.42497008
17 6.89429523 8.68922773
18 -9.40804518 6.89429523
19 1.78720054 -9.40804518
20 -10.46234475 1.78720054
21 -5.49769419 -10.46234475
22 -1.28283469 -5.49769419
23 -14.42508577 -1.28283469
24 5.66526629 -14.42508577
25 31.00650887 5.66526629
26 -8.79869773 31.00650887
27 2.75429703 -8.79869773
28 4.40687390 2.75429703
29 -1.75073157 4.40687390
30 9.77598360 -1.75073157
31 3.50824504 9.77598360
32 -8.88583396 3.50824504
33 8.40287349 -8.88583396
34 -6.88614441 8.40287349
35 -13.05318394 -6.88614441
36 21.10280316 -13.05318394
37 22.32967634 21.10280316
38 6.02814110 22.32967634
39 16.03405536 6.02814110
40 12.99208512 16.03405536
41 5.86225138 12.99208512
42 12.29256274 5.86225138
43 -6.44898292 12.29256274
44 -6.99550302 -6.44898292
45 -1.23260410 -6.99550302
46 -7.92998716 -1.23260410
47 -14.57855060 -7.92998716
48 20.04983536 -14.57855060
49 10.72561051 20.04983536
50 -5.71131021 10.72561051
51 9.18295075 -5.71131021
52 -11.69646050 9.18295075
53 9.87599978 -11.69646050
54 7.93025588 9.87599978
55 -6.16912109 7.93025588
56 -11.80932038 -6.16912109
57 3.44996308 -11.80932038
58 -15.60254839 3.44996308
59 2.66005988 -15.60254839
60 10.16745992 2.66005988
61 2.39754238 10.16745992
62 -7.25430694 2.39754238
63 0.05495899 -7.25430694
64 -19.25479439 0.05495899
65 -1.03058652 -19.25479439
66 -15.31978940 -1.03058652
67 NA -15.31978940
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 8.99592457 -1.53176890
[2,] 12.48629472 8.99592457
[3,] -6.25526609 12.48629472
[4,] -2.31027114 -6.25526609
[5,] -14.45702445 -2.31027114
[6,] -7.85984975 -14.45702445
[7,] 4.13076108 -7.85984975
[8,] -10.64910192 4.13076108
[9,] -9.08935430 -10.64910192
[10,] -11.48829429 -9.08935430
[11,] -9.89592190 -11.48829429
[12,] -0.74612209 -9.89592190
[13,] 24.43463140 -0.74612209
[14,] -13.88218846 24.43463140
[15,] -6.42497008 -13.88218846
[16,] 8.68922773 -6.42497008
[17,] 6.89429523 8.68922773
[18,] -9.40804518 6.89429523
[19,] 1.78720054 -9.40804518
[20,] -10.46234475 1.78720054
[21,] -5.49769419 -10.46234475
[22,] -1.28283469 -5.49769419
[23,] -14.42508577 -1.28283469
[24,] 5.66526629 -14.42508577
[25,] 31.00650887 5.66526629
[26,] -8.79869773 31.00650887
[27,] 2.75429703 -8.79869773
[28,] 4.40687390 2.75429703
[29,] -1.75073157 4.40687390
[30,] 9.77598360 -1.75073157
[31,] 3.50824504 9.77598360
[32,] -8.88583396 3.50824504
[33,] 8.40287349 -8.88583396
[34,] -6.88614441 8.40287349
[35,] -13.05318394 -6.88614441
[36,] 21.10280316 -13.05318394
[37,] 22.32967634 21.10280316
[38,] 6.02814110 22.32967634
[39,] 16.03405536 6.02814110
[40,] 12.99208512 16.03405536
[41,] 5.86225138 12.99208512
[42,] 12.29256274 5.86225138
[43,] -6.44898292 12.29256274
[44,] -6.99550302 -6.44898292
[45,] -1.23260410 -6.99550302
[46,] -7.92998716 -1.23260410
[47,] -14.57855060 -7.92998716
[48,] 20.04983536 -14.57855060
[49,] 10.72561051 20.04983536
[50,] -5.71131021 10.72561051
[51,] 9.18295075 -5.71131021
[52,] -11.69646050 9.18295075
[53,] 9.87599978 -11.69646050
[54,] 7.93025588 9.87599978
[55,] -6.16912109 7.93025588
[56,] -11.80932038 -6.16912109
[57,] 3.44996308 -11.80932038
[58,] -15.60254839 3.44996308
[59,] 2.66005988 -15.60254839
[60,] 10.16745992 2.66005988
[61,] 2.39754238 10.16745992
[62,] -7.25430694 2.39754238
[63,] 0.05495899 -7.25430694
[64,] -19.25479439 0.05495899
[65,] -1.03058652 -19.25479439
[66,] -15.31978940 -1.03058652
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 8.99592457 -1.53176890
2 12.48629472 8.99592457
3 -6.25526609 12.48629472
4 -2.31027114 -6.25526609
5 -14.45702445 -2.31027114
6 -7.85984975 -14.45702445
7 4.13076108 -7.85984975
8 -10.64910192 4.13076108
9 -9.08935430 -10.64910192
10 -11.48829429 -9.08935430
11 -9.89592190 -11.48829429
12 -0.74612209 -9.89592190
13 24.43463140 -0.74612209
14 -13.88218846 24.43463140
15 -6.42497008 -13.88218846
16 8.68922773 -6.42497008
17 6.89429523 8.68922773
18 -9.40804518 6.89429523
19 1.78720054 -9.40804518
20 -10.46234475 1.78720054
21 -5.49769419 -10.46234475
22 -1.28283469 -5.49769419
23 -14.42508577 -1.28283469
24 5.66526629 -14.42508577
25 31.00650887 5.66526629
26 -8.79869773 31.00650887
27 2.75429703 -8.79869773
28 4.40687390 2.75429703
29 -1.75073157 4.40687390
30 9.77598360 -1.75073157
31 3.50824504 9.77598360
32 -8.88583396 3.50824504
33 8.40287349 -8.88583396
34 -6.88614441 8.40287349
35 -13.05318394 -6.88614441
36 21.10280316 -13.05318394
37 22.32967634 21.10280316
38 6.02814110 22.32967634
39 16.03405536 6.02814110
40 12.99208512 16.03405536
41 5.86225138 12.99208512
42 12.29256274 5.86225138
43 -6.44898292 12.29256274
44 -6.99550302 -6.44898292
45 -1.23260410 -6.99550302
46 -7.92998716 -1.23260410
47 -14.57855060 -7.92998716
48 20.04983536 -14.57855060
49 10.72561051 20.04983536
50 -5.71131021 10.72561051
51 9.18295075 -5.71131021
52 -11.69646050 9.18295075
53 9.87599978 -11.69646050
54 7.93025588 9.87599978
55 -6.16912109 7.93025588
56 -11.80932038 -6.16912109
57 3.44996308 -11.80932038
58 -15.60254839 3.44996308
59 2.66005988 -15.60254839
60 10.16745992 2.66005988
61 2.39754238 10.16745992
62 -7.25430694 2.39754238
63 0.05495899 -7.25430694
64 -19.25479439 0.05495899
65 -1.03058652 -19.25479439
66 -15.31978940 -1.03058652
> 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/72rtv1260713218.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/8lb511260713218.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/9qmhe1260713218.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/10f72h1260713218.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/1149il1260713218.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/12r1hc1260713218.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/1398ab1260713218.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/14bzpe1260713219.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/15paxb1260713219.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/16jmbc1260713219.tab")
+ }
>
> try(system("convert tmp/1u14h1260713218.ps tmp/1u14h1260713218.png",intern=TRUE))
character(0)
> try(system("convert tmp/2w3nw1260713218.ps tmp/2w3nw1260713218.png",intern=TRUE))
character(0)
> try(system("convert tmp/3cdmn1260713218.ps tmp/3cdmn1260713218.png",intern=TRUE))
character(0)
> try(system("convert tmp/4vrsm1260713218.ps tmp/4vrsm1260713218.png",intern=TRUE))
character(0)
> try(system("convert tmp/5i8t31260713218.ps tmp/5i8t31260713218.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ctnt1260713218.ps tmp/6ctnt1260713218.png",intern=TRUE))
character(0)
> try(system("convert tmp/72rtv1260713218.ps tmp/72rtv1260713218.png",intern=TRUE))
character(0)
> try(system("convert tmp/8lb511260713218.ps tmp/8lb511260713218.png",intern=TRUE))
character(0)
> try(system("convert tmp/9qmhe1260713218.ps tmp/9qmhe1260713218.png",intern=TRUE))
character(0)
> try(system("convert tmp/10f72h1260713218.ps tmp/10f72h1260713218.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.575 1.570 2.936