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(500857
+ ,1.1
+ ,509127
+ ,509933
+ ,517009
+ ,519164
+ ,506971
+ ,1.6
+ ,500857
+ ,509127
+ ,509933
+ ,517009
+ ,569323
+ ,1.5
+ ,506971
+ ,500857
+ ,509127
+ ,509933
+ ,579714
+ ,1.6
+ ,569323
+ ,506971
+ ,500857
+ ,509127
+ ,577992
+ ,1.7
+ ,579714
+ ,569323
+ ,506971
+ ,500857
+ ,565464
+ ,1.6
+ ,577992
+ ,579714
+ ,569323
+ ,506971
+ ,547344
+ ,1.7
+ ,565464
+ ,577992
+ ,579714
+ ,569323
+ ,554788
+ ,1.6
+ ,547344
+ ,565464
+ ,577992
+ ,579714
+ ,562325
+ ,1.6
+ ,554788
+ ,547344
+ ,565464
+ ,577992
+ ,560854
+ ,1.3
+ ,562325
+ ,554788
+ ,547344
+ ,565464
+ ,555332
+ ,1.1
+ ,560854
+ ,562325
+ ,554788
+ ,547344
+ ,543599
+ ,1.6
+ ,555332
+ ,560854
+ ,562325
+ ,554788
+ ,536662
+ ,1.9
+ ,543599
+ ,555332
+ ,560854
+ ,562325
+ ,542722
+ ,1.6
+ ,536662
+ ,543599
+ ,555332
+ ,560854
+ ,593530
+ ,1.7
+ ,542722
+ ,536662
+ ,543599
+ ,555332
+ ,610763
+ ,1.6
+ ,593530
+ ,542722
+ ,536662
+ ,543599
+ ,612613
+ ,1.4
+ ,610763
+ ,593530
+ ,542722
+ ,536662
+ ,611324
+ ,2.1
+ ,612613
+ ,610763
+ ,593530
+ ,542722
+ ,594167
+ ,1.9
+ ,611324
+ ,612613
+ ,610763
+ ,593530
+ ,595454
+ ,1.7
+ ,594167
+ ,611324
+ ,612613
+ ,610763
+ ,590865
+ ,1.8
+ ,595454
+ ,594167
+ ,611324
+ ,612613
+ ,589379
+ ,2
+ ,590865
+ ,595454
+ ,594167
+ ,611324
+ ,584428
+ ,2.5
+ ,589379
+ ,590865
+ ,595454
+ ,594167
+ ,573100
+ ,2.1
+ ,584428
+ ,589379
+ ,590865
+ ,595454
+ ,567456
+ ,2.1
+ ,573100
+ ,584428
+ ,589379
+ ,590865
+ ,569028
+ ,2.3
+ ,567456
+ ,573100
+ ,584428
+ ,589379
+ ,620735
+ ,2.4
+ ,569028
+ ,567456
+ ,573100
+ ,584428
+ ,628884
+ ,2.4
+ ,620735
+ ,569028
+ ,567456
+ ,573100
+ ,628232
+ ,2.3
+ ,628884
+ ,620735
+ ,569028
+ ,567456
+ ,612117
+ ,1.7
+ ,628232
+ ,628884
+ ,620735
+ ,569028
+ ,595404
+ ,2
+ ,612117
+ ,628232
+ ,628884
+ ,620735
+ ,597141
+ ,2.3
+ ,595404
+ ,612117
+ ,628232
+ ,628884
+ ,593408
+ ,2
+ ,597141
+ ,595404
+ ,612117
+ ,628232
+ ,590072
+ ,2
+ ,593408
+ ,597141
+ ,595404
+ ,612117
+ ,579799
+ ,1.3
+ ,590072
+ ,593408
+ ,597141
+ ,595404
+ ,574205
+ ,1.7
+ ,579799
+ ,590072
+ ,593408
+ ,597141
+ ,572775
+ ,1.9
+ ,574205
+ ,579799
+ ,590072
+ ,593408
+ ,572942
+ ,1.7
+ ,572775
+ ,574205
+ ,579799
+ ,590072
+ ,619567
+ ,1.6
+ ,572942
+ ,572775
+ ,574205
+ ,579799
+ ,625809
+ ,1.7
+ ,619567
+ ,572942
+ ,572775
+ ,574205
+ ,619916
+ ,1.8
+ ,625809
+ ,619567
+ ,572942
+ ,572775
+ ,587625
+ ,1.9
+ ,619916
+ ,625809
+ ,619567
+ ,572942
+ ,565742
+ ,1.9
+ ,587625
+ ,619916
+ ,625809
+ ,619567
+ ,557274
+ ,1.9
+ ,565742
+ ,587625
+ ,619916
+ ,625809
+ ,560576
+ ,2
+ ,557274
+ ,565742
+ ,587625
+ ,619916
+ ,548854
+ ,2.1
+ ,560576
+ ,557274
+ ,565742
+ ,587625
+ ,531673
+ ,1.9
+ ,548854
+ ,560576
+ ,557274
+ ,565742
+ ,525919
+ ,1.9
+ ,531673
+ ,548854
+ ,560576
+ ,557274
+ ,511038
+ ,1.3
+ ,525919
+ ,531673
+ ,548854
+ ,560576
+ ,498662
+ ,1.3
+ ,511038
+ ,525919
+ ,531673
+ ,548854
+ ,555362
+ ,1.4
+ ,498662
+ ,511038
+ ,525919
+ ,531673
+ ,564591
+ ,1.2
+ ,555362
+ ,498662
+ ,511038
+ ,525919
+ ,541657
+ ,1.3
+ ,564591
+ ,555362
+ ,498662
+ ,511038
+ ,527070
+ ,1.8
+ ,541657
+ ,564591
+ ,555362
+ ,498662
+ ,509846
+ ,2.2
+ ,527070
+ ,541657
+ ,564591
+ ,555362
+ ,514258
+ ,2.6
+ ,509846
+ ,527070
+ ,541657
+ ,564591
+ ,516922
+ ,2.8
+ ,514258
+ ,509846
+ ,527070
+ ,541657
+ ,507561
+ ,3.1
+ ,516922
+ ,514258
+ ,509846
+ ,527070
+ ,492622
+ ,3.9
+ ,507561
+ ,516922
+ ,514258
+ ,509846
+ ,490243
+ ,3.7
+ ,492622
+ ,507561
+ ,516922
+ ,514258
+ ,469357
+ ,4.6
+ ,490243
+ ,492622
+ ,507561
+ ,516922
+ ,477580
+ ,5.1
+ ,469357
+ ,490243
+ ,492622
+ ,507561
+ ,528379
+ ,5.2
+ ,477580
+ ,469357
+ ,490243
+ ,492622
+ ,533590
+ ,4.9
+ ,528379
+ ,477580
+ ,469357
+ ,490243
+ ,517945
+ ,5.1
+ ,533590
+ ,528379
+ ,477580
+ ,469357
+ ,506174
+ ,4.8
+ ,517945
+ ,533590
+ ,528379
+ ,477580
+ ,501866
+ ,3.9
+ ,506174
+ ,517945
+ ,533590
+ ,528379
+ ,516141
+ ,3.5
+ ,501866
+ ,506174
+ ,517945
+ ,533590)
+ ,dim=c(6
+ ,68)
+ ,dimnames=list(c('TWIB'
+ ,'GI'
+ ,'TWIB1'
+ ,'TWIB2'
+ ,'TWIB3'
+ ,'TWIB4')
+ ,1:68))
> y <- array(NA,dim=c(6,68),dimnames=list(c('TWIB','GI','TWIB1','TWIB2','TWIB3','TWIB4'),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 = '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
TWIB GI TWIB1 TWIB2 TWIB3 TWIB4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 500857 1.1 509127 509933 517009 519164 1 0 0 0 0 0 0 0 0 0 0 1
2 506971 1.6 500857 509127 509933 517009 0 1 0 0 0 0 0 0 0 0 0 2
3 569323 1.5 506971 500857 509127 509933 0 0 1 0 0 0 0 0 0 0 0 3
4 579714 1.6 569323 506971 500857 509127 0 0 0 1 0 0 0 0 0 0 0 4
5 577992 1.7 579714 569323 506971 500857 0 0 0 0 1 0 0 0 0 0 0 5
6 565464 1.6 577992 579714 569323 506971 0 0 0 0 0 1 0 0 0 0 0 6
7 547344 1.7 565464 577992 579714 569323 0 0 0 0 0 0 1 0 0 0 0 7
8 554788 1.6 547344 565464 577992 579714 0 0 0 0 0 0 0 1 0 0 0 8
9 562325 1.6 554788 547344 565464 577992 0 0 0 0 0 0 0 0 1 0 0 9
10 560854 1.3 562325 554788 547344 565464 0 0 0 0 0 0 0 0 0 1 0 10
11 555332 1.1 560854 562325 554788 547344 0 0 0 0 0 0 0 0 0 0 1 11
12 543599 1.6 555332 560854 562325 554788 0 0 0 0 0 0 0 0 0 0 0 12
13 536662 1.9 543599 555332 560854 562325 1 0 0 0 0 0 0 0 0 0 0 13
14 542722 1.6 536662 543599 555332 560854 0 1 0 0 0 0 0 0 0 0 0 14
15 593530 1.7 542722 536662 543599 555332 0 0 1 0 0 0 0 0 0 0 0 15
16 610763 1.6 593530 542722 536662 543599 0 0 0 1 0 0 0 0 0 0 0 16
17 612613 1.4 610763 593530 542722 536662 0 0 0 0 1 0 0 0 0 0 0 17
18 611324 2.1 612613 610763 593530 542722 0 0 0 0 0 1 0 0 0 0 0 18
19 594167 1.9 611324 612613 610763 593530 0 0 0 0 0 0 1 0 0 0 0 19
20 595454 1.7 594167 611324 612613 610763 0 0 0 0 0 0 0 1 0 0 0 20
21 590865 1.8 595454 594167 611324 612613 0 0 0 0 0 0 0 0 1 0 0 21
22 589379 2.0 590865 595454 594167 611324 0 0 0 0 0 0 0 0 0 1 0 22
23 584428 2.5 589379 590865 595454 594167 0 0 0 0 0 0 0 0 0 0 1 23
24 573100 2.1 584428 589379 590865 595454 0 0 0 0 0 0 0 0 0 0 0 24
25 567456 2.1 573100 584428 589379 590865 1 0 0 0 0 0 0 0 0 0 0 25
26 569028 2.3 567456 573100 584428 589379 0 1 0 0 0 0 0 0 0 0 0 26
27 620735 2.4 569028 567456 573100 584428 0 0 1 0 0 0 0 0 0 0 0 27
28 628884 2.4 620735 569028 567456 573100 0 0 0 1 0 0 0 0 0 0 0 28
29 628232 2.3 628884 620735 569028 567456 0 0 0 0 1 0 0 0 0 0 0 29
30 612117 1.7 628232 628884 620735 569028 0 0 0 0 0 1 0 0 0 0 0 30
31 595404 2.0 612117 628232 628884 620735 0 0 0 0 0 0 1 0 0 0 0 31
32 597141 2.3 595404 612117 628232 628884 0 0 0 0 0 0 0 1 0 0 0 32
33 593408 2.0 597141 595404 612117 628232 0 0 0 0 0 0 0 0 1 0 0 33
34 590072 2.0 593408 597141 595404 612117 0 0 0 0 0 0 0 0 0 1 0 34
35 579799 1.3 590072 593408 597141 595404 0 0 0 0 0 0 0 0 0 0 1 35
36 574205 1.7 579799 590072 593408 597141 0 0 0 0 0 0 0 0 0 0 0 36
37 572775 1.9 574205 579799 590072 593408 1 0 0 0 0 0 0 0 0 0 0 37
38 572942 1.7 572775 574205 579799 590072 0 1 0 0 0 0 0 0 0 0 0 38
39 619567 1.6 572942 572775 574205 579799 0 0 1 0 0 0 0 0 0 0 0 39
40 625809 1.7 619567 572942 572775 574205 0 0 0 1 0 0 0 0 0 0 0 40
41 619916 1.8 625809 619567 572942 572775 0 0 0 0 1 0 0 0 0 0 0 41
42 587625 1.9 619916 625809 619567 572942 0 0 0 0 0 1 0 0 0 0 0 42
43 565742 1.9 587625 619916 625809 619567 0 0 0 0 0 0 1 0 0 0 0 43
44 557274 1.9 565742 587625 619916 625809 0 0 0 0 0 0 0 1 0 0 0 44
45 560576 2.0 557274 565742 587625 619916 0 0 0 0 0 0 0 0 1 0 0 45
46 548854 2.1 560576 557274 565742 587625 0 0 0 0 0 0 0 0 0 1 0 46
47 531673 1.9 548854 560576 557274 565742 0 0 0 0 0 0 0 0 0 0 1 47
48 525919 1.9 531673 548854 560576 557274 0 0 0 0 0 0 0 0 0 0 0 48
49 511038 1.3 525919 531673 548854 560576 1 0 0 0 0 0 0 0 0 0 0 49
50 498662 1.3 511038 525919 531673 548854 0 1 0 0 0 0 0 0 0 0 0 50
51 555362 1.4 498662 511038 525919 531673 0 0 1 0 0 0 0 0 0 0 0 51
52 564591 1.2 555362 498662 511038 525919 0 0 0 1 0 0 0 0 0 0 0 52
53 541657 1.3 564591 555362 498662 511038 0 0 0 0 1 0 0 0 0 0 0 53
54 527070 1.8 541657 564591 555362 498662 0 0 0 0 0 1 0 0 0 0 0 54
55 509846 2.2 527070 541657 564591 555362 0 0 0 0 0 0 1 0 0 0 0 55
56 514258 2.6 509846 527070 541657 564591 0 0 0 0 0 0 0 1 0 0 0 56
57 516922 2.8 514258 509846 527070 541657 0 0 0 0 0 0 0 0 1 0 0 57
58 507561 3.1 516922 514258 509846 527070 0 0 0 0 0 0 0 0 0 1 0 58
59 492622 3.9 507561 516922 514258 509846 0 0 0 0 0 0 0 0 0 0 1 59
60 490243 3.7 492622 507561 516922 514258 0 0 0 0 0 0 0 0 0 0 0 60
61 469357 4.6 490243 492622 507561 516922 1 0 0 0 0 0 0 0 0 0 0 61
62 477580 5.1 469357 490243 492622 507561 0 1 0 0 0 0 0 0 0 0 0 62
63 528379 5.2 477580 469357 490243 492622 0 0 1 0 0 0 0 0 0 0 0 63
64 533590 4.9 528379 477580 469357 490243 0 0 0 1 0 0 0 0 0 0 0 64
65 517945 5.1 533590 528379 477580 469357 0 0 0 0 1 0 0 0 0 0 0 65
66 506174 4.8 517945 533590 528379 477580 0 0 0 0 0 1 0 0 0 0 0 66
67 501866 3.9 506174 517945 533590 528379 0 0 0 0 0 0 1 0 0 0 0 67
68 516141 3.5 501866 506174 517945 533590 0 0 0 0 0 0 0 1 0 0 0 68
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) GI TWIB1 TWIB2 TWIB3 TWIB4
1.728e+04 1.316e+03 9.939e-01 -7.127e-02 6.746e-02 -2.877e-02
M1 M2 M3 M4 M5 M6
-3.746e+03 7.599e+03 5.872e+04 1.614e+04 2.817e+03 -7.256e+03
M7 M8 M9 M10 M11 t
-7.886e+03 1.126e+04 8.263e+03 2.729e+03 -2.820e+03 -1.717e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-14508.78 -3598.96 -19.11 3355.82 12371.26
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.728e+04 1.756e+04 0.984 0.3300
GI 1.316e+03 1.187e+03 1.109 0.2726
TWIB1 9.939e-01 1.474e-01 6.741 1.53e-08 ***
TWIB2 -7.127e-02 2.032e-01 -0.351 0.7273
TWIB3 6.746e-02 2.041e-01 0.330 0.7424
TWIB4 -2.877e-02 1.637e-01 -0.176 0.8612
M1 -3.746e+03 4.271e+03 -0.877 0.3845
M2 7.599e+03 4.475e+03 1.698 0.0957 .
M3 5.872e+04 4.549e+03 12.908 < 2e-16 ***
M4 1.614e+04 1.013e+04 1.593 0.1175
M5 2.817e+03 1.077e+04 0.262 0.7947
M6 -7.256e+03 1.009e+04 -0.719 0.4753
M7 -7.886e+03 4.242e+03 -1.859 0.0689 .
M8 1.126e+04 4.639e+03 2.428 0.0188 *
M9 8.263e+03 5.706e+03 1.448 0.1538
M10 2.729e+03 5.404e+03 0.505 0.6158
M11 -2.820e+03 4.350e+03 -0.648 0.5197
t -1.717e+02 7.776e+01 -2.208 0.0319 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6547 on 50 degrees of freedom
Multiple R-squared: 0.9807, Adjusted R-squared: 0.9742
F-statistic: 149.7 on 17 and 50 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.244755102 0.489510204 0.7552449
[2,] 0.131062324 0.262124649 0.8689377
[3,] 0.131429522 0.262859045 0.8685705
[4,] 0.095927835 0.191855670 0.9040722
[5,] 0.064499332 0.128998664 0.9355007
[6,] 0.070463214 0.140926428 0.9295368
[7,] 0.072112129 0.144224259 0.9278879
[8,] 0.063360463 0.126720925 0.9366395
[9,] 0.046981888 0.093963776 0.9530181
[10,] 0.049642275 0.099284550 0.9503577
[11,] 0.028902637 0.057805275 0.9710974
[12,] 0.015604163 0.031208327 0.9843958
[13,] 0.010730726 0.021461451 0.9892693
[14,] 0.005970051 0.011940103 0.9940299
[15,] 0.003422450 0.006844899 0.9965776
[16,] 0.003103764 0.006207528 0.9968962
[17,] 0.012834720 0.025669440 0.9871653
[18,] 0.015370644 0.030741288 0.9846294
[19,] 0.010284530 0.020569061 0.9897155
[20,] 0.006556352 0.013112704 0.9934436
[21,] 0.174451306 0.348902612 0.8255487
[22,] 0.379965149 0.759930298 0.6200349
[23,] 0.297880621 0.595761243 0.7021194
[24,] 0.366985552 0.733971104 0.6330144
[25,] 0.309924113 0.619848226 0.6900759
[26,] 0.199407649 0.398815298 0.8005924
[27,] 0.117564038 0.235128076 0.8824360
> postscript(file="/var/www/html/rcomp/tmp/15d3p1258743727.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/2c6fp1258743727.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/36w6n1258743727.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/490l11258743727.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/58or71258743727.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
-3550.745018 -691.863167 4026.298722 -3959.595230 1146.608727
6 7 8 9 10
-2583.013136 -6611.572470 -480.389728 2332.429067 863.427418
11 12 13 14 15
2301.349951 -7649.234096 520.740179 2189.802003 -3968.787075
16 17 18 19 20
6215.598531 7707.888283 11879.053383 -2500.871970 -2595.318882
21 22 23 24 25
-6507.599516 3221.617515 3903.049854 -4385.795546 4762.311796
26 27 28 29 30
-9.829637 -727.683237 -1047.802456 7243.187519 -51.217214
31 32 33 34 35
550.228874 -1342.710419 -3360.038502 3507.199764 2328.020042
36 37 38 39 40
3833.031480 11003.226716 1878.771742 -2501.078997 -28.381557
41 42 43 44 45
4507.542950 -14508.779238 -2996.776303 -10415.714085 4789.877662
46 47 48 49 50
-4696.280865 -4065.935004 3305.355643 -1487.768639 -9836.615974
51 52 53 54 55
6914.848719 2765.012690 -11531.650559 2738.470219 -338.461240
56 57 58 59 60
2462.876055 2745.331289 -2895.963832 -4466.484843 4896.642520
61 62 63 64 65
-11247.765035 6469.735032 -3743.598132 -3944.831977 -9073.576920
66 67 68
2525.485986 11897.453109 12371.257059
> postscript(file="/var/www/html/rcomp/tmp/6dqwj1258743727.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 -3550.745018 NA
1 -691.863167 -3550.745018
2 4026.298722 -691.863167
3 -3959.595230 4026.298722
4 1146.608727 -3959.595230
5 -2583.013136 1146.608727
6 -6611.572470 -2583.013136
7 -480.389728 -6611.572470
8 2332.429067 -480.389728
9 863.427418 2332.429067
10 2301.349951 863.427418
11 -7649.234096 2301.349951
12 520.740179 -7649.234096
13 2189.802003 520.740179
14 -3968.787075 2189.802003
15 6215.598531 -3968.787075
16 7707.888283 6215.598531
17 11879.053383 7707.888283
18 -2500.871970 11879.053383
19 -2595.318882 -2500.871970
20 -6507.599516 -2595.318882
21 3221.617515 -6507.599516
22 3903.049854 3221.617515
23 -4385.795546 3903.049854
24 4762.311796 -4385.795546
25 -9.829637 4762.311796
26 -727.683237 -9.829637
27 -1047.802456 -727.683237
28 7243.187519 -1047.802456
29 -51.217214 7243.187519
30 550.228874 -51.217214
31 -1342.710419 550.228874
32 -3360.038502 -1342.710419
33 3507.199764 -3360.038502
34 2328.020042 3507.199764
35 3833.031480 2328.020042
36 11003.226716 3833.031480
37 1878.771742 11003.226716
38 -2501.078997 1878.771742
39 -28.381557 -2501.078997
40 4507.542950 -28.381557
41 -14508.779238 4507.542950
42 -2996.776303 -14508.779238
43 -10415.714085 -2996.776303
44 4789.877662 -10415.714085
45 -4696.280865 4789.877662
46 -4065.935004 -4696.280865
47 3305.355643 -4065.935004
48 -1487.768639 3305.355643
49 -9836.615974 -1487.768639
50 6914.848719 -9836.615974
51 2765.012690 6914.848719
52 -11531.650559 2765.012690
53 2738.470219 -11531.650559
54 -338.461240 2738.470219
55 2462.876055 -338.461240
56 2745.331289 2462.876055
57 -2895.963832 2745.331289
58 -4466.484843 -2895.963832
59 4896.642520 -4466.484843
60 -11247.765035 4896.642520
61 6469.735032 -11247.765035
62 -3743.598132 6469.735032
63 -3944.831977 -3743.598132
64 -9073.576920 -3944.831977
65 2525.485986 -9073.576920
66 11897.453109 2525.485986
67 12371.257059 11897.453109
68 NA 12371.257059
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -691.863167 -3550.745018
[2,] 4026.298722 -691.863167
[3,] -3959.595230 4026.298722
[4,] 1146.608727 -3959.595230
[5,] -2583.013136 1146.608727
[6,] -6611.572470 -2583.013136
[7,] -480.389728 -6611.572470
[8,] 2332.429067 -480.389728
[9,] 863.427418 2332.429067
[10,] 2301.349951 863.427418
[11,] -7649.234096 2301.349951
[12,] 520.740179 -7649.234096
[13,] 2189.802003 520.740179
[14,] -3968.787075 2189.802003
[15,] 6215.598531 -3968.787075
[16,] 7707.888283 6215.598531
[17,] 11879.053383 7707.888283
[18,] -2500.871970 11879.053383
[19,] -2595.318882 -2500.871970
[20,] -6507.599516 -2595.318882
[21,] 3221.617515 -6507.599516
[22,] 3903.049854 3221.617515
[23,] -4385.795546 3903.049854
[24,] 4762.311796 -4385.795546
[25,] -9.829637 4762.311796
[26,] -727.683237 -9.829637
[27,] -1047.802456 -727.683237
[28,] 7243.187519 -1047.802456
[29,] -51.217214 7243.187519
[30,] 550.228874 -51.217214
[31,] -1342.710419 550.228874
[32,] -3360.038502 -1342.710419
[33,] 3507.199764 -3360.038502
[34,] 2328.020042 3507.199764
[35,] 3833.031480 2328.020042
[36,] 11003.226716 3833.031480
[37,] 1878.771742 11003.226716
[38,] -2501.078997 1878.771742
[39,] -28.381557 -2501.078997
[40,] 4507.542950 -28.381557
[41,] -14508.779238 4507.542950
[42,] -2996.776303 -14508.779238
[43,] -10415.714085 -2996.776303
[44,] 4789.877662 -10415.714085
[45,] -4696.280865 4789.877662
[46,] -4065.935004 -4696.280865
[47,] 3305.355643 -4065.935004
[48,] -1487.768639 3305.355643
[49,] -9836.615974 -1487.768639
[50,] 6914.848719 -9836.615974
[51,] 2765.012690 6914.848719
[52,] -11531.650559 2765.012690
[53,] 2738.470219 -11531.650559
[54,] -338.461240 2738.470219
[55,] 2462.876055 -338.461240
[56,] 2745.331289 2462.876055
[57,] -2895.963832 2745.331289
[58,] -4466.484843 -2895.963832
[59,] 4896.642520 -4466.484843
[60,] -11247.765035 4896.642520
[61,] 6469.735032 -11247.765035
[62,] -3743.598132 6469.735032
[63,] -3944.831977 -3743.598132
[64,] -9073.576920 -3944.831977
[65,] 2525.485986 -9073.576920
[66,] 11897.453109 2525.485986
[67,] 12371.257059 11897.453109
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -691.863167 -3550.745018
2 4026.298722 -691.863167
3 -3959.595230 4026.298722
4 1146.608727 -3959.595230
5 -2583.013136 1146.608727
6 -6611.572470 -2583.013136
7 -480.389728 -6611.572470
8 2332.429067 -480.389728
9 863.427418 2332.429067
10 2301.349951 863.427418
11 -7649.234096 2301.349951
12 520.740179 -7649.234096
13 2189.802003 520.740179
14 -3968.787075 2189.802003
15 6215.598531 -3968.787075
16 7707.888283 6215.598531
17 11879.053383 7707.888283
18 -2500.871970 11879.053383
19 -2595.318882 -2500.871970
20 -6507.599516 -2595.318882
21 3221.617515 -6507.599516
22 3903.049854 3221.617515
23 -4385.795546 3903.049854
24 4762.311796 -4385.795546
25 -9.829637 4762.311796
26 -727.683237 -9.829637
27 -1047.802456 -727.683237
28 7243.187519 -1047.802456
29 -51.217214 7243.187519
30 550.228874 -51.217214
31 -1342.710419 550.228874
32 -3360.038502 -1342.710419
33 3507.199764 -3360.038502
34 2328.020042 3507.199764
35 3833.031480 2328.020042
36 11003.226716 3833.031480
37 1878.771742 11003.226716
38 -2501.078997 1878.771742
39 -28.381557 -2501.078997
40 4507.542950 -28.381557
41 -14508.779238 4507.542950
42 -2996.776303 -14508.779238
43 -10415.714085 -2996.776303
44 4789.877662 -10415.714085
45 -4696.280865 4789.877662
46 -4065.935004 -4696.280865
47 3305.355643 -4065.935004
48 -1487.768639 3305.355643
49 -9836.615974 -1487.768639
50 6914.848719 -9836.615974
51 2765.012690 6914.848719
52 -11531.650559 2765.012690
53 2738.470219 -11531.650559
54 -338.461240 2738.470219
55 2462.876055 -338.461240
56 2745.331289 2462.876055
57 -2895.963832 2745.331289
58 -4466.484843 -2895.963832
59 4896.642520 -4466.484843
60 -11247.765035 4896.642520
61 6469.735032 -11247.765035
62 -3743.598132 6469.735032
63 -3944.831977 -3743.598132
64 -9073.576920 -3944.831977
65 2525.485986 -9073.576920
66 11897.453109 2525.485986
67 12371.257059 11897.453109
> 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/7waz51258743727.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/8mt951258743727.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/9usy21258743727.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/101l921258743727.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/11xls81258743727.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/12diwz1258743727.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/13twh81258743727.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/143hbz1258743727.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/15f5mn1258743727.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/16tjbc1258743727.tab")
+ }
>
> system("convert tmp/15d3p1258743727.ps tmp/15d3p1258743727.png")
> system("convert tmp/2c6fp1258743727.ps tmp/2c6fp1258743727.png")
> system("convert tmp/36w6n1258743727.ps tmp/36w6n1258743727.png")
> system("convert tmp/490l11258743727.ps tmp/490l11258743727.png")
> system("convert tmp/58or71258743727.ps tmp/58or71258743727.png")
> system("convert tmp/6dqwj1258743727.ps tmp/6dqwj1258743727.png")
> system("convert tmp/7waz51258743727.ps tmp/7waz51258743727.png")
> system("convert tmp/8mt951258743727.ps tmp/8mt951258743727.png")
> system("convert tmp/9usy21258743727.ps tmp/9usy21258743727.png")
> system("convert tmp/101l921258743727.ps tmp/101l921258743727.png")
>
>
> proc.time()
user system elapsed
2.518 1.543 2.935