R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
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(467,0,460,0,448,0,443,0,436,0,431,0,484,0,510,0,513,0,503,0,471,0,471,0,476,0,475,0,470,0,461,0,455,0,456,0,517,0,525,0,523,0,519,0,509,0,512,0,519,0,517,0,510,0,509,0,501,0,507,0,569,0,580,0,578,0,565,0,547,0,555,0,562,0,561,0,555,1,544,1,537,1,543,1,594,1,611,1,613,1,611,1,594,1,595,1,591,1,589,1,584,1,573,1,567,1,569,1,621,1,629,1,628,1,612,1,595,1,597,1,593,1,590,1,580,1,574,1,573,1,573,1,620,1,626,1,620,1,588,1,566,1,557,1,561,1,549,1,532,1,526,1,511,1,499,1,555,1,565,1,542,1,527,1,510,1,514,1),dim=c(2,84),dimnames=list(c('Y','DUM'),1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('Y','DUM'),1:84))
> 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 = 'Include Monthly 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 DUM M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 467 0 1 0 0 0 0 0 0 0 0 0 0 1
2 460 0 0 1 0 0 0 0 0 0 0 0 0 2
3 448 0 0 0 1 0 0 0 0 0 0 0 0 3
4 443 0 0 0 0 1 0 0 0 0 0 0 0 4
5 436 0 0 0 0 0 1 0 0 0 0 0 0 5
6 431 0 0 0 0 0 0 1 0 0 0 0 0 6
7 484 0 0 0 0 0 0 0 1 0 0 0 0 7
8 510 0 0 0 0 0 0 0 0 1 0 0 0 8
9 513 0 0 0 0 0 0 0 0 0 1 0 0 9
10 503 0 0 0 0 0 0 0 0 0 0 1 0 10
11 471 0 0 0 0 0 0 0 0 0 0 0 1 11
12 471 0 0 0 0 0 0 0 0 0 0 0 0 12
13 476 0 1 0 0 0 0 0 0 0 0 0 0 13
14 475 0 0 1 0 0 0 0 0 0 0 0 0 14
15 470 0 0 0 1 0 0 0 0 0 0 0 0 15
16 461 0 0 0 0 1 0 0 0 0 0 0 0 16
17 455 0 0 0 0 0 1 0 0 0 0 0 0 17
18 456 0 0 0 0 0 0 1 0 0 0 0 0 18
19 517 0 0 0 0 0 0 0 1 0 0 0 0 19
20 525 0 0 0 0 0 0 0 0 1 0 0 0 20
21 523 0 0 0 0 0 0 0 0 0 1 0 0 21
22 519 0 0 0 0 0 0 0 0 0 0 1 0 22
23 509 0 0 0 0 0 0 0 0 0 0 0 1 23
24 512 0 0 0 0 0 0 0 0 0 0 0 0 24
25 519 0 1 0 0 0 0 0 0 0 0 0 0 25
26 517 0 0 1 0 0 0 0 0 0 0 0 0 26
27 510 0 0 0 1 0 0 0 0 0 0 0 0 27
28 509 0 0 0 0 1 0 0 0 0 0 0 0 28
29 501 0 0 0 0 0 1 0 0 0 0 0 0 29
30 507 0 0 0 0 0 0 1 0 0 0 0 0 30
31 569 0 0 0 0 0 0 0 1 0 0 0 0 31
32 580 0 0 0 0 0 0 0 0 1 0 0 0 32
33 578 0 0 0 0 0 0 0 0 0 1 0 0 33
34 565 0 0 0 0 0 0 0 0 0 0 1 0 34
35 547 0 0 0 0 0 0 0 0 0 0 0 1 35
36 555 0 0 0 0 0 0 0 0 0 0 0 0 36
37 562 0 1 0 0 0 0 0 0 0 0 0 0 37
38 561 0 0 1 0 0 0 0 0 0 0 0 0 38
39 555 1 0 0 1 0 0 0 0 0 0 0 0 39
40 544 1 0 0 0 1 0 0 0 0 0 0 0 40
41 537 1 0 0 0 0 1 0 0 0 0 0 0 41
42 543 1 0 0 0 0 0 1 0 0 0 0 0 42
43 594 1 0 0 0 0 0 0 1 0 0 0 0 43
44 611 1 0 0 0 0 0 0 0 1 0 0 0 44
45 613 1 0 0 0 0 0 0 0 0 1 0 0 45
46 611 1 0 0 0 0 0 0 0 0 0 1 0 46
47 594 1 0 0 0 0 0 0 0 0 0 0 1 47
48 595 1 0 0 0 0 0 0 0 0 0 0 0 48
49 591 1 1 0 0 0 0 0 0 0 0 0 0 49
50 589 1 0 1 0 0 0 0 0 0 0 0 0 50
51 584 1 0 0 1 0 0 0 0 0 0 0 0 51
52 573 1 0 0 0 1 0 0 0 0 0 0 0 52
53 567 1 0 0 0 0 1 0 0 0 0 0 0 53
54 569 1 0 0 0 0 0 1 0 0 0 0 0 54
55 621 1 0 0 0 0 0 0 1 0 0 0 0 55
56 629 1 0 0 0 0 0 0 0 1 0 0 0 56
57 628 1 0 0 0 0 0 0 0 0 1 0 0 57
58 612 1 0 0 0 0 0 0 0 0 0 1 0 58
59 595 1 0 0 0 0 0 0 0 0 0 0 1 59
60 597 1 0 0 0 0 0 0 0 0 0 0 0 60
61 593 1 1 0 0 0 0 0 0 0 0 0 0 61
62 590 1 0 1 0 0 0 0 0 0 0 0 0 62
63 580 1 0 0 1 0 0 0 0 0 0 0 0 63
64 574 1 0 0 0 1 0 0 0 0 0 0 0 64
65 573 1 0 0 0 0 1 0 0 0 0 0 0 65
66 573 1 0 0 0 0 0 1 0 0 0 0 0 66
67 620 1 0 0 0 0 0 0 1 0 0 0 0 67
68 626 1 0 0 0 0 0 0 0 1 0 0 0 68
69 620 1 0 0 0 0 0 0 0 0 1 0 0 69
70 588 1 0 0 0 0 0 0 0 0 0 1 0 70
71 566 1 0 0 0 0 0 0 0 0 0 0 1 71
72 557 1 0 0 0 0 0 0 0 0 0 0 0 72
73 561 1 1 0 0 0 0 0 0 0 0 0 0 73
74 549 1 0 1 0 0 0 0 0 0 0 0 0 74
75 532 1 0 0 1 0 0 0 0 0 0 0 0 75
76 526 1 0 0 0 1 0 0 0 0 0 0 0 76
77 511 1 0 0 0 0 1 0 0 0 0 0 0 77
78 499 1 0 0 0 0 0 1 0 0 0 0 0 78
79 555 1 0 0 0 0 0 0 1 0 0 0 0 79
80 565 1 0 0 0 0 0 0 0 1 0 0 0 80
81 542 1 0 0 0 0 0 0 0 0 1 0 0 81
82 527 1 0 0 0 0 0 0 0 0 0 1 0 82
83 510 1 0 0 0 0 0 0 0 0 0 0 1 83
84 514 1 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) DUM M1 M2 M3 M4
498.6111 64.8056 6.3725 2.2192 -16.0491 -23.2024
M5 M6 M7 M8 M9 M10
-30.4985 -30.9375 23.4807 35.6131 31.3170 18.0208
M11 t
-1.1324 0.1533
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-67.006 -26.113 4.655 24.266 54.345
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 498.6111 15.0173 33.203 < 2e-16 ***
DUM 64.8056 14.5776 4.446 3.22e-05 ***
M1 6.3725 17.7651 0.359 0.7209
M2 2.2192 17.7422 0.125 0.9008
M3 -16.0491 17.8766 -0.898 0.3724
M4 -23.2024 17.8336 -1.301 0.1975
M5 -30.4985 17.7956 -1.714 0.0910 .
M6 -30.9375 17.7626 -1.742 0.0859 .
M7 23.4807 17.7346 1.324 0.1898
M8 35.6131 17.7116 2.011 0.0482 *
M9 31.3170 17.6937 1.770 0.0811 .
M10 18.0208 17.6810 1.019 0.3116
M11 -1.1324 17.6733 -0.064 0.9491
t 0.1533 0.3006 0.510 0.6117
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 33.06 on 70 degrees of freedom
Multiple R-squared: 0.6572, Adjusted R-squared: 0.5936
F-statistic: 10.33 on 13 and 70 DF, p-value: 1.001e-11
> 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,] 4.436822e-03 8.873644e-03 0.9955632
[2,] 1.515829e-03 3.031659e-03 0.9984842
[3,] 1.714401e-03 3.428802e-03 0.9982856
[4,] 4.964630e-04 9.929260e-04 0.9995035
[5,] 2.373371e-04 4.746743e-04 0.9997627
[6,] 6.308356e-05 1.261671e-04 0.9999369
[7,] 1.749548e-04 3.499096e-04 0.9998250
[8,] 3.504417e-04 7.008834e-04 0.9996496
[9,] 3.374961e-04 6.749923e-04 0.9996625
[10,] 3.118306e-04 6.236612e-04 0.9996882
[11,] 2.239047e-04 4.478093e-04 0.9997761
[12,] 2.317868e-04 4.635735e-04 0.9997682
[13,] 1.716414e-04 3.432829e-04 0.9998284
[14,] 2.252780e-04 4.505559e-04 0.9997747
[15,] 3.440320e-04 6.880640e-04 0.9996560
[16,] 2.448558e-04 4.897117e-04 0.9997551
[17,] 1.367222e-04 2.734443e-04 0.9998633
[18,] 5.832701e-05 1.166540e-04 0.9999417
[19,] 2.631508e-05 5.263017e-05 0.9999737
[20,] 1.553211e-05 3.106422e-05 0.9999845
[21,] 6.534950e-06 1.306990e-05 0.9999935
[22,] 2.864318e-06 5.728635e-06 0.9999971
[23,] 2.142067e-06 4.284133e-06 0.9999979
[24,] 2.397879e-06 4.795759e-06 0.9999976
[25,] 3.571870e-06 7.143739e-06 0.9999964
[26,] 4.877430e-06 9.754859e-06 0.9999951
[27,] 1.140335e-05 2.280670e-05 0.9999886
[28,] 2.384241e-05 4.768483e-05 0.9999762
[29,] 3.512894e-05 7.025788e-05 0.9999649
[30,] 2.822914e-05 5.645828e-05 0.9999718
[31,] 2.657927e-05 5.315854e-05 0.9999734
[32,] 2.238547e-05 4.477094e-05 0.9999776
[33,] 4.440393e-05 8.880787e-05 0.9999556
[34,] 7.237852e-05 1.447570e-04 0.9999276
[35,] 7.215081e-05 1.443016e-04 0.9999278
[36,] 1.541561e-04 3.083122e-04 0.9998458
[37,] 3.966526e-04 7.933052e-04 0.9996033
[38,] 7.326793e-04 1.465359e-03 0.9992673
[39,] 2.359511e-03 4.719023e-03 0.9976405
[40,] 1.905965e-02 3.811930e-02 0.9809404
[41,] 5.182582e-02 1.036516e-01 0.9481742
[42,] 1.259888e-01 2.519777e-01 0.8740112
[43,] 2.346742e-01 4.693485e-01 0.7653258
[44,] 3.813401e-01 7.626803e-01 0.6186599
[45,] 6.642027e-01 6.715945e-01 0.3357973
[46,] 7.794388e-01 4.411223e-01 0.2205612
[47,] 8.034683e-01 3.930634e-01 0.1965317
[48,] 8.227101e-01 3.545798e-01 0.1772899
[49,] 7.304022e-01 5.391956e-01 0.2695978
[50,] 6.970121e-01 6.059758e-01 0.3029879
[51,] 5.726259e-01 8.547481e-01 0.4273741
> postscript(file="/var/www/html/rcomp/tmp/1tne01229364084.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/2ad0u1229364084.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/3r7t61229364084.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/433ht1229364084.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/5iejk1229364084.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 = 84
Frequency = 1
1 2 3 4 5 6
-38.1369048 -41.1369048 -35.0218254 -33.0218254 -32.8789683 -37.5932540
7 8 9 10 11 12
-39.1646825 -25.4503968 -18.3075397 -15.1646825 -28.1646825 -29.4503968
13 14 15 16 17 18
-30.9761905 -27.9761905 -14.8611111 -16.8611111 -15.7182540 -14.4325397
19 20 21 22 23 24
-8.0039683 -12.2896825 -10.1468254 -1.0039683 7.9960317 9.7103175
25 26 27 28 29 30
10.1845238 12.1845238 23.2996032 29.2996032 28.4424603 34.7281746
31 32 33 34 35 36
42.1567460 40.8710317 43.0138889 43.1567460 44.1567460 50.8710317
37 38 39 40 41 42
51.3452381 54.3452381 1.6547619 -2.3452381 -2.2023810 4.0833333
43 44 45 46 47 48
0.5119048 5.2261905 11.3690476 22.5119048 24.5119048 24.2261905
49 50 51 52 53 54
13.7003968 15.7003968 28.8154762 24.8154762 25.9583333 28.2440476
55 56 57 58 59 60
25.6726190 21.3869048 24.5297619 21.6726190 23.6726190 24.3869048
61 62 63 64 65 66
13.8611111 14.8611111 22.9761905 23.9761905 30.1190476 30.4047619
67 68 69 70 71 72
22.8333333 16.5476190 14.6904762 -4.1666667 -7.1666667 -17.4523810
73 74 75 76 77 78
-19.9781746 -27.9781746 -26.8630952 -25.8630952 -33.7202381 -45.4345238
79 80 81 82 83 84
-44.0059524 -46.2916667 -65.1488095 -67.0059524 -65.0059524 -62.2916667
> postscript(file="/var/www/html/rcomp/tmp/6ilpv1229364084.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 -38.1369048 NA
1 -41.1369048 -38.1369048
2 -35.0218254 -41.1369048
3 -33.0218254 -35.0218254
4 -32.8789683 -33.0218254
5 -37.5932540 -32.8789683
6 -39.1646825 -37.5932540
7 -25.4503968 -39.1646825
8 -18.3075397 -25.4503968
9 -15.1646825 -18.3075397
10 -28.1646825 -15.1646825
11 -29.4503968 -28.1646825
12 -30.9761905 -29.4503968
13 -27.9761905 -30.9761905
14 -14.8611111 -27.9761905
15 -16.8611111 -14.8611111
16 -15.7182540 -16.8611111
17 -14.4325397 -15.7182540
18 -8.0039683 -14.4325397
19 -12.2896825 -8.0039683
20 -10.1468254 -12.2896825
21 -1.0039683 -10.1468254
22 7.9960317 -1.0039683
23 9.7103175 7.9960317
24 10.1845238 9.7103175
25 12.1845238 10.1845238
26 23.2996032 12.1845238
27 29.2996032 23.2996032
28 28.4424603 29.2996032
29 34.7281746 28.4424603
30 42.1567460 34.7281746
31 40.8710317 42.1567460
32 43.0138889 40.8710317
33 43.1567460 43.0138889
34 44.1567460 43.1567460
35 50.8710317 44.1567460
36 51.3452381 50.8710317
37 54.3452381 51.3452381
38 1.6547619 54.3452381
39 -2.3452381 1.6547619
40 -2.2023810 -2.3452381
41 4.0833333 -2.2023810
42 0.5119048 4.0833333
43 5.2261905 0.5119048
44 11.3690476 5.2261905
45 22.5119048 11.3690476
46 24.5119048 22.5119048
47 24.2261905 24.5119048
48 13.7003968 24.2261905
49 15.7003968 13.7003968
50 28.8154762 15.7003968
51 24.8154762 28.8154762
52 25.9583333 24.8154762
53 28.2440476 25.9583333
54 25.6726190 28.2440476
55 21.3869048 25.6726190
56 24.5297619 21.3869048
57 21.6726190 24.5297619
58 23.6726190 21.6726190
59 24.3869048 23.6726190
60 13.8611111 24.3869048
61 14.8611111 13.8611111
62 22.9761905 14.8611111
63 23.9761905 22.9761905
64 30.1190476 23.9761905
65 30.4047619 30.1190476
66 22.8333333 30.4047619
67 16.5476190 22.8333333
68 14.6904762 16.5476190
69 -4.1666667 14.6904762
70 -7.1666667 -4.1666667
71 -17.4523810 -7.1666667
72 -19.9781746 -17.4523810
73 -27.9781746 -19.9781746
74 -26.8630952 -27.9781746
75 -25.8630952 -26.8630952
76 -33.7202381 -25.8630952
77 -45.4345238 -33.7202381
78 -44.0059524 -45.4345238
79 -46.2916667 -44.0059524
80 -65.1488095 -46.2916667
81 -67.0059524 -65.1488095
82 -65.0059524 -67.0059524
83 -62.2916667 -65.0059524
84 NA -62.2916667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -41.1369048 -38.1369048
[2,] -35.0218254 -41.1369048
[3,] -33.0218254 -35.0218254
[4,] -32.8789683 -33.0218254
[5,] -37.5932540 -32.8789683
[6,] -39.1646825 -37.5932540
[7,] -25.4503968 -39.1646825
[8,] -18.3075397 -25.4503968
[9,] -15.1646825 -18.3075397
[10,] -28.1646825 -15.1646825
[11,] -29.4503968 -28.1646825
[12,] -30.9761905 -29.4503968
[13,] -27.9761905 -30.9761905
[14,] -14.8611111 -27.9761905
[15,] -16.8611111 -14.8611111
[16,] -15.7182540 -16.8611111
[17,] -14.4325397 -15.7182540
[18,] -8.0039683 -14.4325397
[19,] -12.2896825 -8.0039683
[20,] -10.1468254 -12.2896825
[21,] -1.0039683 -10.1468254
[22,] 7.9960317 -1.0039683
[23,] 9.7103175 7.9960317
[24,] 10.1845238 9.7103175
[25,] 12.1845238 10.1845238
[26,] 23.2996032 12.1845238
[27,] 29.2996032 23.2996032
[28,] 28.4424603 29.2996032
[29,] 34.7281746 28.4424603
[30,] 42.1567460 34.7281746
[31,] 40.8710317 42.1567460
[32,] 43.0138889 40.8710317
[33,] 43.1567460 43.0138889
[34,] 44.1567460 43.1567460
[35,] 50.8710317 44.1567460
[36,] 51.3452381 50.8710317
[37,] 54.3452381 51.3452381
[38,] 1.6547619 54.3452381
[39,] -2.3452381 1.6547619
[40,] -2.2023810 -2.3452381
[41,] 4.0833333 -2.2023810
[42,] 0.5119048 4.0833333
[43,] 5.2261905 0.5119048
[44,] 11.3690476 5.2261905
[45,] 22.5119048 11.3690476
[46,] 24.5119048 22.5119048
[47,] 24.2261905 24.5119048
[48,] 13.7003968 24.2261905
[49,] 15.7003968 13.7003968
[50,] 28.8154762 15.7003968
[51,] 24.8154762 28.8154762
[52,] 25.9583333 24.8154762
[53,] 28.2440476 25.9583333
[54,] 25.6726190 28.2440476
[55,] 21.3869048 25.6726190
[56,] 24.5297619 21.3869048
[57,] 21.6726190 24.5297619
[58,] 23.6726190 21.6726190
[59,] 24.3869048 23.6726190
[60,] 13.8611111 24.3869048
[61,] 14.8611111 13.8611111
[62,] 22.9761905 14.8611111
[63,] 23.9761905 22.9761905
[64,] 30.1190476 23.9761905
[65,] 30.4047619 30.1190476
[66,] 22.8333333 30.4047619
[67,] 16.5476190 22.8333333
[68,] 14.6904762 16.5476190
[69,] -4.1666667 14.6904762
[70,] -7.1666667 -4.1666667
[71,] -17.4523810 -7.1666667
[72,] -19.9781746 -17.4523810
[73,] -27.9781746 -19.9781746
[74,] -26.8630952 -27.9781746
[75,] -25.8630952 -26.8630952
[76,] -33.7202381 -25.8630952
[77,] -45.4345238 -33.7202381
[78,] -44.0059524 -45.4345238
[79,] -46.2916667 -44.0059524
[80,] -65.1488095 -46.2916667
[81,] -67.0059524 -65.1488095
[82,] -65.0059524 -67.0059524
[83,] -62.2916667 -65.0059524
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -41.1369048 -38.1369048
2 -35.0218254 -41.1369048
3 -33.0218254 -35.0218254
4 -32.8789683 -33.0218254
5 -37.5932540 -32.8789683
6 -39.1646825 -37.5932540
7 -25.4503968 -39.1646825
8 -18.3075397 -25.4503968
9 -15.1646825 -18.3075397
10 -28.1646825 -15.1646825
11 -29.4503968 -28.1646825
12 -30.9761905 -29.4503968
13 -27.9761905 -30.9761905
14 -14.8611111 -27.9761905
15 -16.8611111 -14.8611111
16 -15.7182540 -16.8611111
17 -14.4325397 -15.7182540
18 -8.0039683 -14.4325397
19 -12.2896825 -8.0039683
20 -10.1468254 -12.2896825
21 -1.0039683 -10.1468254
22 7.9960317 -1.0039683
23 9.7103175 7.9960317
24 10.1845238 9.7103175
25 12.1845238 10.1845238
26 23.2996032 12.1845238
27 29.2996032 23.2996032
28 28.4424603 29.2996032
29 34.7281746 28.4424603
30 42.1567460 34.7281746
31 40.8710317 42.1567460
32 43.0138889 40.8710317
33 43.1567460 43.0138889
34 44.1567460 43.1567460
35 50.8710317 44.1567460
36 51.3452381 50.8710317
37 54.3452381 51.3452381
38 1.6547619 54.3452381
39 -2.3452381 1.6547619
40 -2.2023810 -2.3452381
41 4.0833333 -2.2023810
42 0.5119048 4.0833333
43 5.2261905 0.5119048
44 11.3690476 5.2261905
45 22.5119048 11.3690476
46 24.5119048 22.5119048
47 24.2261905 24.5119048
48 13.7003968 24.2261905
49 15.7003968 13.7003968
50 28.8154762 15.7003968
51 24.8154762 28.8154762
52 25.9583333 24.8154762
53 28.2440476 25.9583333
54 25.6726190 28.2440476
55 21.3869048 25.6726190
56 24.5297619 21.3869048
57 21.6726190 24.5297619
58 23.6726190 21.6726190
59 24.3869048 23.6726190
60 13.8611111 24.3869048
61 14.8611111 13.8611111
62 22.9761905 14.8611111
63 23.9761905 22.9761905
64 30.1190476 23.9761905
65 30.4047619 30.1190476
66 22.8333333 30.4047619
67 16.5476190 22.8333333
68 14.6904762 16.5476190
69 -4.1666667 14.6904762
70 -7.1666667 -4.1666667
71 -17.4523810 -7.1666667
72 -19.9781746 -17.4523810
73 -27.9781746 -19.9781746
74 -26.8630952 -27.9781746
75 -25.8630952 -26.8630952
76 -33.7202381 -25.8630952
77 -45.4345238 -33.7202381
78 -44.0059524 -45.4345238
79 -46.2916667 -44.0059524
80 -65.1488095 -46.2916667
81 -67.0059524 -65.1488095
82 -65.0059524 -67.0059524
83 -62.2916667 -65.0059524
> 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/769uc1229364084.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/839j01229364084.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/9yc471229364084.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/104p5p1229364084.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/11zfbc1229364084.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/12min81229364085.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/133w291229364085.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/14i0xe1229364085.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/150qwx1229364085.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/16rsnr1229364085.tab")
+ }
>
> system("convert tmp/1tne01229364084.ps tmp/1tne01229364084.png")
> system("convert tmp/2ad0u1229364084.ps tmp/2ad0u1229364084.png")
> system("convert tmp/3r7t61229364084.ps tmp/3r7t61229364084.png")
> system("convert tmp/433ht1229364084.ps tmp/433ht1229364084.png")
> system("convert tmp/5iejk1229364084.ps tmp/5iejk1229364084.png")
> system("convert tmp/6ilpv1229364084.ps tmp/6ilpv1229364084.png")
> system("convert tmp/769uc1229364084.ps tmp/769uc1229364084.png")
> system("convert tmp/839j01229364084.ps tmp/839j01229364084.png")
> system("convert tmp/9yc471229364084.ps tmp/9yc471229364084.png")
> system("convert tmp/104p5p1229364084.ps tmp/104p5p1229364084.png")
>
>
> proc.time()
user system elapsed
2.732 1.578 3.199