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(6.7
+ ,510
+ ,6.9
+ ,7.0
+ ,6.7
+ ,509
+ ,6.7
+ ,6.9
+ ,6.5
+ ,501
+ ,6.7
+ ,6.7
+ ,6.4
+ ,507
+ ,6.5
+ ,6.7
+ ,6.5
+ ,569
+ ,6.4
+ ,6.5
+ ,6.5
+ ,580
+ ,6.5
+ ,6.4
+ ,6.5
+ ,578
+ ,6.5
+ ,6.5
+ ,6.7
+ ,565
+ ,6.5
+ ,6.5
+ ,6.8
+ ,547
+ ,6.7
+ ,6.5
+ ,7.2
+ ,555
+ ,6.8
+ ,6.7
+ ,7.6
+ ,562
+ ,7.2
+ ,6.8
+ ,7.6
+ ,561
+ ,7.6
+ ,7.2
+ ,7.2
+ ,555
+ ,7.6
+ ,7.6
+ ,6.4
+ ,544
+ ,7.2
+ ,7.6
+ ,6.1
+ ,537
+ ,6.4
+ ,7.2
+ ,6.3
+ ,543
+ ,6.1
+ ,6.4
+ ,7.1
+ ,594
+ ,6.3
+ ,6.1
+ ,7.5
+ ,611
+ ,7.1
+ ,6.3
+ ,7.4
+ ,613
+ ,7.5
+ ,7.1
+ ,7.1
+ ,611
+ ,7.4
+ ,7.5
+ ,6.8
+ ,594
+ ,7.1
+ ,7.4
+ ,6.9
+ ,595
+ ,6.8
+ ,7.1
+ ,7.2
+ ,591
+ ,6.9
+ ,6.8
+ ,7.4
+ ,589
+ ,7.2
+ ,6.9
+ ,7.3
+ ,584
+ ,7.4
+ ,7.2
+ ,6.9
+ ,573
+ ,7.3
+ ,7.4
+ ,6.9
+ ,567
+ ,6.9
+ ,7.3
+ ,6.8
+ ,569
+ ,6.9
+ ,6.9
+ ,7.1
+ ,621
+ ,6.8
+ ,6.9
+ ,7.2
+ ,629
+ ,7.1
+ ,6.8
+ ,7.1
+ ,628
+ ,7.2
+ ,7.1
+ ,7.0
+ ,612
+ ,7.1
+ ,7.2
+ ,6.9
+ ,595
+ ,7.0
+ ,7.1
+ ,7.1
+ ,597
+ ,6.9
+ ,7.0
+ ,7.3
+ ,593
+ ,7.1
+ ,6.9
+ ,7.5
+ ,590
+ ,7.3
+ ,7.1
+ ,7.5
+ ,580
+ ,7.5
+ ,7.3
+ ,7.5
+ ,574
+ ,7.5
+ ,7.5
+ ,7.3
+ ,573
+ ,7.5
+ ,7.5
+ ,7.0
+ ,573
+ ,7.3
+ ,7.5
+ ,6.7
+ ,620
+ ,7.0
+ ,7.3
+ ,6.5
+ ,626
+ ,6.7
+ ,7.0
+ ,6.5
+ ,620
+ ,6.5
+ ,6.7
+ ,6.5
+ ,588
+ ,6.5
+ ,6.5
+ ,6.6
+ ,566
+ ,6.5
+ ,6.5
+ ,6.8
+ ,557
+ ,6.6
+ ,6.5
+ ,6.9
+ ,561
+ ,6.8
+ ,6.6
+ ,6.9
+ ,549
+ ,6.9
+ ,6.8
+ ,6.8
+ ,532
+ ,6.9
+ ,6.9
+ ,6.8
+ ,526
+ ,6.8
+ ,6.9
+ ,6.5
+ ,511
+ ,6.8
+ ,6.8
+ ,6.1
+ ,499
+ ,6.5
+ ,6.8
+ ,6.1
+ ,555
+ ,6.1
+ ,6.5
+ ,5.9
+ ,565
+ ,6.1
+ ,6.1
+ ,5.7
+ ,542
+ ,5.9
+ ,6.1
+ ,5.9
+ ,527
+ ,5.7
+ ,5.9
+ ,5.9
+ ,510
+ ,5.9
+ ,5.7
+ ,6.1
+ ,514
+ ,5.9
+ ,5.9
+ ,6.3
+ ,517
+ ,6.1
+ ,5.9
+ ,6.2
+ ,508
+ ,6.3
+ ,6.1
+ ,5.9
+ ,493
+ ,6.2
+ ,6.3
+ ,5.7
+ ,490
+ ,5.9
+ ,6.2
+ ,5.4
+ ,469
+ ,5.7
+ ,5.9
+ ,5.6
+ ,478
+ ,5.4
+ ,5.7
+ ,6.2
+ ,528
+ ,5.6
+ ,5.4
+ ,6.3
+ ,534
+ ,6.2
+ ,5.6
+ ,6.0
+ ,518
+ ,6.3
+ ,6.2
+ ,5.6
+ ,506
+ ,6.0
+ ,6.3
+ ,5.5
+ ,502
+ ,5.6
+ ,6.0
+ ,5.9
+ ,516
+ ,5.5
+ ,5.6)
+ ,dim=c(4
+ ,70)
+ ,dimnames=list(c('wkgo'
+ ,'werkl'
+ ,'Y1'
+ ,'Y2')
+ ,1:70))
> y <- array(NA,dim=c(4,70),dimnames=list(c('wkgo','werkl','Y1','Y2'),1:70))
> 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
wkgo werkl Y1 Y2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 6.7 510 6.9 7.0 1 0 0 0 0 0 0 0 0 0 0 1
2 6.7 509 6.7 6.9 0 1 0 0 0 0 0 0 0 0 0 2
3 6.5 501 6.7 6.7 0 0 1 0 0 0 0 0 0 0 0 3
4 6.4 507 6.5 6.7 0 0 0 1 0 0 0 0 0 0 0 4
5 6.5 569 6.4 6.5 0 0 0 0 1 0 0 0 0 0 0 5
6 6.5 580 6.5 6.4 0 0 0 0 0 1 0 0 0 0 0 6
7 6.5 578 6.5 6.5 0 0 0 0 0 0 1 0 0 0 0 7
8 6.7 565 6.5 6.5 0 0 0 0 0 0 0 1 0 0 0 8
9 6.8 547 6.7 6.5 0 0 0 0 0 0 0 0 1 0 0 9
10 7.2 555 6.8 6.7 0 0 0 0 0 0 0 0 0 1 0 10
11 7.6 562 7.2 6.8 0 0 0 0 0 0 0 0 0 0 1 11
12 7.6 561 7.6 7.2 0 0 0 0 0 0 0 0 0 0 0 12
13 7.2 555 7.6 7.6 1 0 0 0 0 0 0 0 0 0 0 13
14 6.4 544 7.2 7.6 0 1 0 0 0 0 0 0 0 0 0 14
15 6.1 537 6.4 7.2 0 0 1 0 0 0 0 0 0 0 0 15
16 6.3 543 6.1 6.4 0 0 0 1 0 0 0 0 0 0 0 16
17 7.1 594 6.3 6.1 0 0 0 0 1 0 0 0 0 0 0 17
18 7.5 611 7.1 6.3 0 0 0 0 0 1 0 0 0 0 0 18
19 7.4 613 7.5 7.1 0 0 0 0 0 0 1 0 0 0 0 19
20 7.1 611 7.4 7.5 0 0 0 0 0 0 0 1 0 0 0 20
21 6.8 594 7.1 7.4 0 0 0 0 0 0 0 0 1 0 0 21
22 6.9 595 6.8 7.1 0 0 0 0 0 0 0 0 0 1 0 22
23 7.2 591 6.9 6.8 0 0 0 0 0 0 0 0 0 0 1 23
24 7.4 589 7.2 6.9 0 0 0 0 0 0 0 0 0 0 0 24
25 7.3 584 7.4 7.2 1 0 0 0 0 0 0 0 0 0 0 25
26 6.9 573 7.3 7.4 0 1 0 0 0 0 0 0 0 0 0 26
27 6.9 567 6.9 7.3 0 0 1 0 0 0 0 0 0 0 0 27
28 6.8 569 6.9 6.9 0 0 0 1 0 0 0 0 0 0 0 28
29 7.1 621 6.8 6.9 0 0 0 0 1 0 0 0 0 0 0 29
30 7.2 629 7.1 6.8 0 0 0 0 0 1 0 0 0 0 0 30
31 7.1 628 7.2 7.1 0 0 0 0 0 0 1 0 0 0 0 31
32 7.0 612 7.1 7.2 0 0 0 0 0 0 0 1 0 0 0 32
33 6.9 595 7.0 7.1 0 0 0 0 0 0 0 0 1 0 0 33
34 7.1 597 6.9 7.0 0 0 0 0 0 0 0 0 0 1 0 34
35 7.3 593 7.1 6.9 0 0 0 0 0 0 0 0 0 0 1 35
36 7.5 590 7.3 7.1 0 0 0 0 0 0 0 0 0 0 0 36
37 7.5 580 7.5 7.3 1 0 0 0 0 0 0 0 0 0 0 37
38 7.5 574 7.5 7.5 0 1 0 0 0 0 0 0 0 0 0 38
39 7.3 573 7.5 7.5 0 0 1 0 0 0 0 0 0 0 0 39
40 7.0 573 7.3 7.5 0 0 0 1 0 0 0 0 0 0 0 40
41 6.7 620 7.0 7.3 0 0 0 0 1 0 0 0 0 0 0 41
42 6.5 626 6.7 7.0 0 0 0 0 0 1 0 0 0 0 0 42
43 6.5 620 6.5 6.7 0 0 0 0 0 0 1 0 0 0 0 43
44 6.5 588 6.5 6.5 0 0 0 0 0 0 0 1 0 0 0 44
45 6.6 566 6.5 6.5 0 0 0 0 0 0 0 0 1 0 0 45
46 6.8 557 6.6 6.5 0 0 0 0 0 0 0 0 0 1 0 46
47 6.9 561 6.8 6.6 0 0 0 0 0 0 0 0 0 0 1 47
48 6.9 549 6.9 6.8 0 0 0 0 0 0 0 0 0 0 0 48
49 6.8 532 6.9 6.9 1 0 0 0 0 0 0 0 0 0 0 49
50 6.8 526 6.8 6.9 0 1 0 0 0 0 0 0 0 0 0 50
51 6.5 511 6.8 6.8 0 0 1 0 0 0 0 0 0 0 0 51
52 6.1 499 6.5 6.8 0 0 0 1 0 0 0 0 0 0 0 52
53 6.1 555 6.1 6.5 0 0 0 0 1 0 0 0 0 0 0 53
54 5.9 565 6.1 6.1 0 0 0 0 0 1 0 0 0 0 0 54
55 5.7 542 5.9 6.1 0 0 0 0 0 0 1 0 0 0 0 55
56 5.9 527 5.7 5.9 0 0 0 0 0 0 0 1 0 0 0 56
57 5.9 510 5.9 5.7 0 0 0 0 0 0 0 0 1 0 0 57
58 6.1 514 5.9 5.9 0 0 0 0 0 0 0 0 0 1 0 58
59 6.3 517 6.1 5.9 0 0 0 0 0 0 0 0 0 0 1 59
60 6.2 508 6.3 6.1 0 0 0 0 0 0 0 0 0 0 0 60
61 5.9 493 6.2 6.3 1 0 0 0 0 0 0 0 0 0 0 61
62 5.7 490 5.9 6.2 0 1 0 0 0 0 0 0 0 0 0 62
63 5.4 469 5.7 5.9 0 0 1 0 0 0 0 0 0 0 0 63
64 5.6 478 5.4 5.7 0 0 0 1 0 0 0 0 0 0 0 64
65 6.2 528 5.6 5.4 0 0 0 0 1 0 0 0 0 0 0 65
66 6.3 534 6.2 5.6 0 0 0 0 0 1 0 0 0 0 0 66
67 6.0 518 6.3 6.2 0 0 0 0 0 0 1 0 0 0 0 67
68 5.6 506 6.0 6.3 0 0 0 0 0 0 0 1 0 0 0 68
69 5.5 502 5.6 6.0 0 0 0 0 0 0 0 0 1 0 0 69
70 5.9 516 5.5 5.6 0 0 0 0 0 0 0 0 0 1 0 70
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) werkl Y1 Y2 M1 M2
0.753077 0.006099 1.208670 -0.797725 0.020268 0.077559
M3 M4 M5 M6 M7 M8
0.059439 0.044486 -0.097049 -0.487497 -0.394446 -0.198209
M9 M10 M11 t
-0.160317 0.080424 0.017829 -0.003819
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.334633 -0.061684 -0.001579 0.066097 0.237672
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.7530767 0.3219296 2.339 0.023051 *
werkl 0.0060988 0.0009408 6.483 2.84e-08 ***
Y1 1.2086700 0.0860435 14.047 < 2e-16 ***
Y2 -0.7977251 0.0883421 -9.030 2.23e-12 ***
M1 0.0202679 0.0846541 0.239 0.811686
M2 0.0775591 0.0902316 0.860 0.393833
M3 0.0594394 0.0912283 0.652 0.517458
M4 0.0444861 0.0891022 0.499 0.619617
M5 -0.0970493 0.1001519 -0.969 0.336855
M6 -0.4874970 0.0988603 -4.931 8.17e-06 ***
M7 -0.3944455 0.0909634 -4.336 6.37e-05 ***
M8 -0.1982094 0.0894608 -2.216 0.030955 *
M9 -0.1603166 0.0863128 -1.857 0.068712 .
M10 0.0804244 0.0879450 0.914 0.364527
M11 0.0178293 0.0854357 0.209 0.835478
t -0.0038185 0.0009754 -3.915 0.000256 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1309 on 54 degrees of freedom
Multiple R-squared: 0.9589, Adjusted R-squared: 0.9475
F-statistic: 83.98 on 15 and 54 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.9303359 0.1393281 0.06966407
[2,] 0.8801375 0.2397249 0.11986246
[3,] 0.8208693 0.3582614 0.17913070
[4,] 0.7703090 0.4593820 0.22969102
[5,] 0.7414594 0.5170812 0.25854061
[6,] 0.6666183 0.6667635 0.33338175
[7,] 0.6269762 0.7460476 0.37302378
[8,] 0.6693202 0.6613596 0.33067979
[9,] 0.8489588 0.3020823 0.15104116
[10,] 0.8634215 0.2731569 0.13657846
[11,] 0.8268114 0.3463772 0.17318859
[12,] 0.7673577 0.4652846 0.23264229
[13,] 0.6918678 0.6162644 0.30813219
[14,] 0.6138018 0.7723964 0.38619820
[15,] 0.5234814 0.9530371 0.47651857
[16,] 0.4493195 0.8986391 0.55068045
[17,] 0.3923111 0.7846221 0.60768894
[18,] 0.3667011 0.7334022 0.63329889
[19,] 0.3017192 0.6034384 0.69828081
[20,] 0.3907101 0.7814203 0.60928985
[21,] 0.3533076 0.7066151 0.64669245
[22,] 0.2794894 0.5589788 0.72051060
[23,] 0.5817105 0.8365790 0.41828951
[24,] 0.4872494 0.9744989 0.51275056
[25,] 0.4066730 0.8133459 0.59332705
[26,] 0.7683506 0.4632988 0.23164938
[27,] 0.7202935 0.5594130 0.27970652
[28,] 0.6572745 0.6854511 0.34272553
[29,] 0.6937236 0.6125527 0.30627635
[30,] 0.5813260 0.8373479 0.41867397
[31,] 0.4729435 0.9458870 0.52705650
[32,] 0.4998262 0.9996524 0.50017380
[33,] 0.4654746 0.9309491 0.53452543
> postscript(file="/var/www/html/rcomp/tmp/1cbzq1258988318.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/24tr41258988318.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/3i7si1258988318.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/468gl1258988318.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/59rot1258988318.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 = 70
Frequency = 1
1 2 3 4 5
0.0643419513 0.1789295459 -0.1098868673 0.0140261812 -0.1574230356
6 7 8 9 10
-0.0308830061 -0.0281459266 0.0587208148 -0.0073092162 0.1456559345
11 12 13 14 15
0.1656825435 0.0290511649 -0.0317154574 -0.3346334203 0.0778423820
16 17 18 19 20
-0.0155576318 0.1377063549 0.0209021305 -0.0258164165 -0.0660793603
21 22 23 24 25
-0.0136456768 -0.0333835049 -0.0027591984 -0.0517423056 -0.1401142350
26 27 28 29 30
-0.2460881909 0.2161383308 -0.1963774883 0.0527062394 0.0558086412
31 32 33 34 35
-0.0088750589 -0.0030724278 0.0076272495 -0.0003983613 -0.0310960415
36 37 38 39 40
0.1266591553 0.0890086818 0.2316737588 0.0597108122 0.0202166176
41 42 43 44 45
-0.2180167008 0.0629402888 0.0127165130 -0.1440846985 0.0560144480
46 47 48 49 50
-0.0468859398 -0.0668289462 0.0666823891 0.1336849632 0.2376720270
51 52 53 54 55
-0.0286803284 0.0258779938 0.0738500350 -0.1119616641 -0.0191884375
56 57 58 59 60
0.1620648794 -0.1696089607 -0.0713816357 -0.0649983574 -0.1706504037
61 62 63 64 65
-0.1152059039 -0.0675537206 -0.2151243293 0.1518143275 0.1111771071
66 67 68 69 70
0.0031936098 0.0693093264 -0.0075492076 0.1269221561 0.0063935070
> postscript(file="/var/www/html/rcomp/tmp/6g6571258988318.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 = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 0.0643419513 NA
1 0.1789295459 0.0643419513
2 -0.1098868673 0.1789295459
3 0.0140261812 -0.1098868673
4 -0.1574230356 0.0140261812
5 -0.0308830061 -0.1574230356
6 -0.0281459266 -0.0308830061
7 0.0587208148 -0.0281459266
8 -0.0073092162 0.0587208148
9 0.1456559345 -0.0073092162
10 0.1656825435 0.1456559345
11 0.0290511649 0.1656825435
12 -0.0317154574 0.0290511649
13 -0.3346334203 -0.0317154574
14 0.0778423820 -0.3346334203
15 -0.0155576318 0.0778423820
16 0.1377063549 -0.0155576318
17 0.0209021305 0.1377063549
18 -0.0258164165 0.0209021305
19 -0.0660793603 -0.0258164165
20 -0.0136456768 -0.0660793603
21 -0.0333835049 -0.0136456768
22 -0.0027591984 -0.0333835049
23 -0.0517423056 -0.0027591984
24 -0.1401142350 -0.0517423056
25 -0.2460881909 -0.1401142350
26 0.2161383308 -0.2460881909
27 -0.1963774883 0.2161383308
28 0.0527062394 -0.1963774883
29 0.0558086412 0.0527062394
30 -0.0088750589 0.0558086412
31 -0.0030724278 -0.0088750589
32 0.0076272495 -0.0030724278
33 -0.0003983613 0.0076272495
34 -0.0310960415 -0.0003983613
35 0.1266591553 -0.0310960415
36 0.0890086818 0.1266591553
37 0.2316737588 0.0890086818
38 0.0597108122 0.2316737588
39 0.0202166176 0.0597108122
40 -0.2180167008 0.0202166176
41 0.0629402888 -0.2180167008
42 0.0127165130 0.0629402888
43 -0.1440846985 0.0127165130
44 0.0560144480 -0.1440846985
45 -0.0468859398 0.0560144480
46 -0.0668289462 -0.0468859398
47 0.0666823891 -0.0668289462
48 0.1336849632 0.0666823891
49 0.2376720270 0.1336849632
50 -0.0286803284 0.2376720270
51 0.0258779938 -0.0286803284
52 0.0738500350 0.0258779938
53 -0.1119616641 0.0738500350
54 -0.0191884375 -0.1119616641
55 0.1620648794 -0.0191884375
56 -0.1696089607 0.1620648794
57 -0.0713816357 -0.1696089607
58 -0.0649983574 -0.0713816357
59 -0.1706504037 -0.0649983574
60 -0.1152059039 -0.1706504037
61 -0.0675537206 -0.1152059039
62 -0.2151243293 -0.0675537206
63 0.1518143275 -0.2151243293
64 0.1111771071 0.1518143275
65 0.0031936098 0.1111771071
66 0.0693093264 0.0031936098
67 -0.0075492076 0.0693093264
68 0.1269221561 -0.0075492076
69 0.0063935070 0.1269221561
70 NA 0.0063935070
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.1789295459 0.0643419513
[2,] -0.1098868673 0.1789295459
[3,] 0.0140261812 -0.1098868673
[4,] -0.1574230356 0.0140261812
[5,] -0.0308830061 -0.1574230356
[6,] -0.0281459266 -0.0308830061
[7,] 0.0587208148 -0.0281459266
[8,] -0.0073092162 0.0587208148
[9,] 0.1456559345 -0.0073092162
[10,] 0.1656825435 0.1456559345
[11,] 0.0290511649 0.1656825435
[12,] -0.0317154574 0.0290511649
[13,] -0.3346334203 -0.0317154574
[14,] 0.0778423820 -0.3346334203
[15,] -0.0155576318 0.0778423820
[16,] 0.1377063549 -0.0155576318
[17,] 0.0209021305 0.1377063549
[18,] -0.0258164165 0.0209021305
[19,] -0.0660793603 -0.0258164165
[20,] -0.0136456768 -0.0660793603
[21,] -0.0333835049 -0.0136456768
[22,] -0.0027591984 -0.0333835049
[23,] -0.0517423056 -0.0027591984
[24,] -0.1401142350 -0.0517423056
[25,] -0.2460881909 -0.1401142350
[26,] 0.2161383308 -0.2460881909
[27,] -0.1963774883 0.2161383308
[28,] 0.0527062394 -0.1963774883
[29,] 0.0558086412 0.0527062394
[30,] -0.0088750589 0.0558086412
[31,] -0.0030724278 -0.0088750589
[32,] 0.0076272495 -0.0030724278
[33,] -0.0003983613 0.0076272495
[34,] -0.0310960415 -0.0003983613
[35,] 0.1266591553 -0.0310960415
[36,] 0.0890086818 0.1266591553
[37,] 0.2316737588 0.0890086818
[38,] 0.0597108122 0.2316737588
[39,] 0.0202166176 0.0597108122
[40,] -0.2180167008 0.0202166176
[41,] 0.0629402888 -0.2180167008
[42,] 0.0127165130 0.0629402888
[43,] -0.1440846985 0.0127165130
[44,] 0.0560144480 -0.1440846985
[45,] -0.0468859398 0.0560144480
[46,] -0.0668289462 -0.0468859398
[47,] 0.0666823891 -0.0668289462
[48,] 0.1336849632 0.0666823891
[49,] 0.2376720270 0.1336849632
[50,] -0.0286803284 0.2376720270
[51,] 0.0258779938 -0.0286803284
[52,] 0.0738500350 0.0258779938
[53,] -0.1119616641 0.0738500350
[54,] -0.0191884375 -0.1119616641
[55,] 0.1620648794 -0.0191884375
[56,] -0.1696089607 0.1620648794
[57,] -0.0713816357 -0.1696089607
[58,] -0.0649983574 -0.0713816357
[59,] -0.1706504037 -0.0649983574
[60,] -0.1152059039 -0.1706504037
[61,] -0.0675537206 -0.1152059039
[62,] -0.2151243293 -0.0675537206
[63,] 0.1518143275 -0.2151243293
[64,] 0.1111771071 0.1518143275
[65,] 0.0031936098 0.1111771071
[66,] 0.0693093264 0.0031936098
[67,] -0.0075492076 0.0693093264
[68,] 0.1269221561 -0.0075492076
[69,] 0.0063935070 0.1269221561
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.1789295459 0.0643419513
2 -0.1098868673 0.1789295459
3 0.0140261812 -0.1098868673
4 -0.1574230356 0.0140261812
5 -0.0308830061 -0.1574230356
6 -0.0281459266 -0.0308830061
7 0.0587208148 -0.0281459266
8 -0.0073092162 0.0587208148
9 0.1456559345 -0.0073092162
10 0.1656825435 0.1456559345
11 0.0290511649 0.1656825435
12 -0.0317154574 0.0290511649
13 -0.3346334203 -0.0317154574
14 0.0778423820 -0.3346334203
15 -0.0155576318 0.0778423820
16 0.1377063549 -0.0155576318
17 0.0209021305 0.1377063549
18 -0.0258164165 0.0209021305
19 -0.0660793603 -0.0258164165
20 -0.0136456768 -0.0660793603
21 -0.0333835049 -0.0136456768
22 -0.0027591984 -0.0333835049
23 -0.0517423056 -0.0027591984
24 -0.1401142350 -0.0517423056
25 -0.2460881909 -0.1401142350
26 0.2161383308 -0.2460881909
27 -0.1963774883 0.2161383308
28 0.0527062394 -0.1963774883
29 0.0558086412 0.0527062394
30 -0.0088750589 0.0558086412
31 -0.0030724278 -0.0088750589
32 0.0076272495 -0.0030724278
33 -0.0003983613 0.0076272495
34 -0.0310960415 -0.0003983613
35 0.1266591553 -0.0310960415
36 0.0890086818 0.1266591553
37 0.2316737588 0.0890086818
38 0.0597108122 0.2316737588
39 0.0202166176 0.0597108122
40 -0.2180167008 0.0202166176
41 0.0629402888 -0.2180167008
42 0.0127165130 0.0629402888
43 -0.1440846985 0.0127165130
44 0.0560144480 -0.1440846985
45 -0.0468859398 0.0560144480
46 -0.0668289462 -0.0468859398
47 0.0666823891 -0.0668289462
48 0.1336849632 0.0666823891
49 0.2376720270 0.1336849632
50 -0.0286803284 0.2376720270
51 0.0258779938 -0.0286803284
52 0.0738500350 0.0258779938
53 -0.1119616641 0.0738500350
54 -0.0191884375 -0.1119616641
55 0.1620648794 -0.0191884375
56 -0.1696089607 0.1620648794
57 -0.0713816357 -0.1696089607
58 -0.0649983574 -0.0713816357
59 -0.1706504037 -0.0649983574
60 -0.1152059039 -0.1706504037
61 -0.0675537206 -0.1152059039
62 -0.2151243293 -0.0675537206
63 0.1518143275 -0.2151243293
64 0.1111771071 0.1518143275
65 0.0031936098 0.1111771071
66 0.0693093264 0.0031936098
67 -0.0075492076 0.0693093264
68 0.1269221561 -0.0075492076
69 0.0063935070 0.1269221561
> 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/73fk41258988318.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/8ryow1258988318.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/9vrxi1258988318.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/10rbi81258988318.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/11codz1258988318.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/12s9mz1258988318.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/13cnsl1258988318.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/14ftrk1258988318.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/15ferw1258988318.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/16quxb1258988318.tab")
+ }
>
> system("convert tmp/1cbzq1258988318.ps tmp/1cbzq1258988318.png")
> system("convert tmp/24tr41258988318.ps tmp/24tr41258988318.png")
> system("convert tmp/3i7si1258988318.ps tmp/3i7si1258988318.png")
> system("convert tmp/468gl1258988318.ps tmp/468gl1258988318.png")
> system("convert tmp/59rot1258988318.ps tmp/59rot1258988318.png")
> system("convert tmp/6g6571258988318.ps tmp/6g6571258988318.png")
> system("convert tmp/73fk41258988318.ps tmp/73fk41258988318.png")
> system("convert tmp/8ryow1258988318.ps tmp/8ryow1258988318.png")
> system("convert tmp/9vrxi1258988318.ps tmp/9vrxi1258988318.png")
> system("convert tmp/10rbi81258988318.ps tmp/10rbi81258988318.png")
>
>
> proc.time()
user system elapsed
2.573 1.605 3.541