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
+ ,506971
+ ,1.6
+ ,500857
+ ,509127
+ ,569323
+ ,1.5
+ ,506971
+ ,500857
+ ,579714
+ ,1.6
+ ,569323
+ ,506971
+ ,577992
+ ,1.7
+ ,579714
+ ,569323
+ ,565464
+ ,1.6
+ ,577992
+ ,579714
+ ,547344
+ ,1.7
+ ,565464
+ ,577992
+ ,554788
+ ,1.6
+ ,547344
+ ,565464
+ ,562325
+ ,1.6
+ ,554788
+ ,547344
+ ,560854
+ ,1.3
+ ,562325
+ ,554788
+ ,555332
+ ,1.1
+ ,560854
+ ,562325
+ ,543599
+ ,1.6
+ ,555332
+ ,560854
+ ,536662
+ ,1.9
+ ,543599
+ ,555332
+ ,542722
+ ,1.6
+ ,536662
+ ,543599
+ ,593530
+ ,1.7
+ ,542722
+ ,536662
+ ,610763
+ ,1.6
+ ,593530
+ ,542722
+ ,612613
+ ,1.4
+ ,610763
+ ,593530
+ ,611324
+ ,2.1
+ ,612613
+ ,610763
+ ,594167
+ ,1.9
+ ,611324
+ ,612613
+ ,595454
+ ,1.7
+ ,594167
+ ,611324
+ ,590865
+ ,1.8
+ ,595454
+ ,594167
+ ,589379
+ ,2
+ ,590865
+ ,595454
+ ,584428
+ ,2.5
+ ,589379
+ ,590865
+ ,573100
+ ,2.1
+ ,584428
+ ,589379
+ ,567456
+ ,2.1
+ ,573100
+ ,584428
+ ,569028
+ ,2.3
+ ,567456
+ ,573100
+ ,620735
+ ,2.4
+ ,569028
+ ,567456
+ ,628884
+ ,2.4
+ ,620735
+ ,569028
+ ,628232
+ ,2.3
+ ,628884
+ ,620735
+ ,612117
+ ,1.7
+ ,628232
+ ,628884
+ ,595404
+ ,2
+ ,612117
+ ,628232
+ ,597141
+ ,2.3
+ ,595404
+ ,612117
+ ,593408
+ ,2
+ ,597141
+ ,595404
+ ,590072
+ ,2
+ ,593408
+ ,597141
+ ,579799
+ ,1.3
+ ,590072
+ ,593408
+ ,574205
+ ,1.7
+ ,579799
+ ,590072
+ ,572775
+ ,1.9
+ ,574205
+ ,579799
+ ,572942
+ ,1.7
+ ,572775
+ ,574205
+ ,619567
+ ,1.6
+ ,572942
+ ,572775
+ ,625809
+ ,1.7
+ ,619567
+ ,572942
+ ,619916
+ ,1.8
+ ,625809
+ ,619567
+ ,587625
+ ,1.9
+ ,619916
+ ,625809
+ ,565742
+ ,1.9
+ ,587625
+ ,619916
+ ,557274
+ ,1.9
+ ,565742
+ ,587625
+ ,560576
+ ,2
+ ,557274
+ ,565742
+ ,548854
+ ,2.1
+ ,560576
+ ,557274
+ ,531673
+ ,1.9
+ ,548854
+ ,560576
+ ,525919
+ ,1.9
+ ,531673
+ ,548854
+ ,511038
+ ,1.3
+ ,525919
+ ,531673
+ ,498662
+ ,1.3
+ ,511038
+ ,525919
+ ,555362
+ ,1.4
+ ,498662
+ ,511038
+ ,564591
+ ,1.2
+ ,555362
+ ,498662
+ ,541657
+ ,1.3
+ ,564591
+ ,555362
+ ,527070
+ ,1.8
+ ,541657
+ ,564591
+ ,509846
+ ,2.2
+ ,527070
+ ,541657
+ ,514258
+ ,2.6
+ ,509846
+ ,527070
+ ,516922
+ ,2.8
+ ,514258
+ ,509846
+ ,507561
+ ,3.1
+ ,516922
+ ,514258
+ ,492622
+ ,3.9
+ ,507561
+ ,516922
+ ,490243
+ ,3.7
+ ,492622
+ ,507561
+ ,469357
+ ,4.6
+ ,490243
+ ,492622
+ ,477580
+ ,5.1
+ ,469357
+ ,490243
+ ,528379
+ ,5.2
+ ,477580
+ ,469357
+ ,533590
+ ,4.9
+ ,528379
+ ,477580
+ ,517945
+ ,5.1
+ ,533590
+ ,528379
+ ,506174
+ ,4.8
+ ,517945
+ ,533590
+ ,501866
+ ,3.9
+ ,506174
+ ,517945
+ ,516141
+ ,3.5
+ ,501866
+ ,506174)
+ ,dim=c(4
+ ,68)
+ ,dimnames=list(c('TWIB'
+ ,'GI'
+ ,'TWIB1'
+ ,'TWIB2
')
+ ,1:68))
> y <- array(NA,dim=c(4,68),dimnames=list(c('TWIB','GI','TWIB1','TWIB2
'),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\r M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 500857 1.1 509127 509933 1 0 0 0 0 0 0 0 0 0 0 1
2 506971 1.6 500857 509127 0 1 0 0 0 0 0 0 0 0 0 2
3 569323 1.5 506971 500857 0 0 1 0 0 0 0 0 0 0 0 3
4 579714 1.6 569323 506971 0 0 0 1 0 0 0 0 0 0 0 4
5 577992 1.7 579714 569323 0 0 0 0 1 0 0 0 0 0 0 5
6 565464 1.6 577992 579714 0 0 0 0 0 1 0 0 0 0 0 6
7 547344 1.7 565464 577992 0 0 0 0 0 0 1 0 0 0 0 7
8 554788 1.6 547344 565464 0 0 0 0 0 0 0 1 0 0 0 8
9 562325 1.6 554788 547344 0 0 0 0 0 0 0 0 1 0 0 9
10 560854 1.3 562325 554788 0 0 0 0 0 0 0 0 0 1 0 10
11 555332 1.1 560854 562325 0 0 0 0 0 0 0 0 0 0 1 11
12 543599 1.6 555332 560854 0 0 0 0 0 0 0 0 0 0 0 12
13 536662 1.9 543599 555332 1 0 0 0 0 0 0 0 0 0 0 13
14 542722 1.6 536662 543599 0 1 0 0 0 0 0 0 0 0 0 14
15 593530 1.7 542722 536662 0 0 1 0 0 0 0 0 0 0 0 15
16 610763 1.6 593530 542722 0 0 0 1 0 0 0 0 0 0 0 16
17 612613 1.4 610763 593530 0 0 0 0 1 0 0 0 0 0 0 17
18 611324 2.1 612613 610763 0 0 0 0 0 1 0 0 0 0 0 18
19 594167 1.9 611324 612613 0 0 0 0 0 0 1 0 0 0 0 19
20 595454 1.7 594167 611324 0 0 0 0 0 0 0 1 0 0 0 20
21 590865 1.8 595454 594167 0 0 0 0 0 0 0 0 1 0 0 21
22 589379 2.0 590865 595454 0 0 0 0 0 0 0 0 0 1 0 22
23 584428 2.5 589379 590865 0 0 0 0 0 0 0 0 0 0 1 23
24 573100 2.1 584428 589379 0 0 0 0 0 0 0 0 0 0 0 24
25 567456 2.1 573100 584428 1 0 0 0 0 0 0 0 0 0 0 25
26 569028 2.3 567456 573100 0 1 0 0 0 0 0 0 0 0 0 26
27 620735 2.4 569028 567456 0 0 1 0 0 0 0 0 0 0 0 27
28 628884 2.4 620735 569028 0 0 0 1 0 0 0 0 0 0 0 28
29 628232 2.3 628884 620735 0 0 0 0 1 0 0 0 0 0 0 29
30 612117 1.7 628232 628884 0 0 0 0 0 1 0 0 0 0 0 30
31 595404 2.0 612117 628232 0 0 0 0 0 0 1 0 0 0 0 31
32 597141 2.3 595404 612117 0 0 0 0 0 0 0 1 0 0 0 32
33 593408 2.0 597141 595404 0 0 0 0 0 0 0 0 1 0 0 33
34 590072 2.0 593408 597141 0 0 0 0 0 0 0 0 0 1 0 34
35 579799 1.3 590072 593408 0 0 0 0 0 0 0 0 0 0 1 35
36 574205 1.7 579799 590072 0 0 0 0 0 0 0 0 0 0 0 36
37 572775 1.9 574205 579799 1 0 0 0 0 0 0 0 0 0 0 37
38 572942 1.7 572775 574205 0 1 0 0 0 0 0 0 0 0 0 38
39 619567 1.6 572942 572775 0 0 1 0 0 0 0 0 0 0 0 39
40 625809 1.7 619567 572942 0 0 0 1 0 0 0 0 0 0 0 40
41 619916 1.8 625809 619567 0 0 0 0 1 0 0 0 0 0 0 41
42 587625 1.9 619916 625809 0 0 0 0 0 1 0 0 0 0 0 42
43 565742 1.9 587625 619916 0 0 0 0 0 0 1 0 0 0 0 43
44 557274 1.9 565742 587625 0 0 0 0 0 0 0 1 0 0 0 44
45 560576 2.0 557274 565742 0 0 0 0 0 0 0 0 1 0 0 45
46 548854 2.1 560576 557274 0 0 0 0 0 0 0 0 0 1 0 46
47 531673 1.9 548854 560576 0 0 0 0 0 0 0 0 0 0 1 47
48 525919 1.9 531673 548854 0 0 0 0 0 0 0 0 0 0 0 48
49 511038 1.3 525919 531673 1 0 0 0 0 0 0 0 0 0 0 49
50 498662 1.3 511038 525919 0 1 0 0 0 0 0 0 0 0 0 50
51 555362 1.4 498662 511038 0 0 1 0 0 0 0 0 0 0 0 51
52 564591 1.2 555362 498662 0 0 0 1 0 0 0 0 0 0 0 52
53 541657 1.3 564591 555362 0 0 0 0 1 0 0 0 0 0 0 53
54 527070 1.8 541657 564591 0 0 0 0 0 1 0 0 0 0 0 54
55 509846 2.2 527070 541657 0 0 0 0 0 0 1 0 0 0 0 55
56 514258 2.6 509846 527070 0 0 0 0 0 0 0 1 0 0 0 56
57 516922 2.8 514258 509846 0 0 0 0 0 0 0 0 1 0 0 57
58 507561 3.1 516922 514258 0 0 0 0 0 0 0 0 0 1 0 58
59 492622 3.9 507561 516922 0 0 0 0 0 0 0 0 0 0 1 59
60 490243 3.7 492622 507561 0 0 0 0 0 0 0 0 0 0 0 60
61 469357 4.6 490243 492622 1 0 0 0 0 0 0 0 0 0 0 61
62 477580 5.1 469357 490243 0 1 0 0 0 0 0 0 0 0 0 62
63 528379 5.2 477580 469357 0 0 1 0 0 0 0 0 0 0 0 63
64 533590 4.9 528379 477580 0 0 0 1 0 0 0 0 0 0 0 64
65 517945 5.1 533590 528379 0 0 0 0 1 0 0 0 0 0 0 65
66 506174 4.8 517945 533590 0 0 0 0 0 1 0 0 0 0 0 66
67 501866 3.9 506174 517945 0 0 0 0 0 0 1 0 0 0 0 67
68 516141 3.5 501866 506174 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\r` M1 M2
1.745e+04 1.308e+03 9.943e-01 -3.293e-02 -3.737e+03 7.320e+03
M3 M4 M5 M6 M7 M8
5.868e+04 1.554e+04 5.578e+02 -6.335e+03 -7.568e+03 1.137e+04
M9 M10 M11 t
8.282e+03 1.911e+03 -3.067e+03 -1.702e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-14589.6 -3468.9 252.8 3547.8 12063.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.745e+04 1.715e+04 1.017 0.31372
GI 1.308e+03 1.079e+03 1.213 0.23072
TWIB1 9.943e-01 1.436e-01 6.926 6.48e-09 ***
`TWIB2\r` -3.293e-02 1.408e-01 -0.234 0.81598
M1 -3.737e+03 3.957e+03 -0.944 0.34932
M2 7.320e+03 3.956e+03 1.850 0.06995 .
M3 5.868e+04 4.269e+03 13.746 < 2e-16 ***
M4 1.554e+04 9.777e+03 1.589 0.11812
M5 5.578e+02 4.929e+03 0.113 0.91032
M6 -6.335e+03 4.044e+03 -1.567 0.12329
M7 -7.568e+03 3.979e+03 -1.902 0.06274 .
M8 1.137e+04 3.961e+03 2.871 0.00590 **
M9 8.282e+03 4.393e+03 1.885 0.06501 .
M10 1.911e+03 4.389e+03 0.435 0.66512
M11 -3.067e+03 4.134e+03 -0.742 0.46160
t -1.702e+02 5.867e+01 -2.901 0.00544 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6427 on 52 degrees of freedom
Multiple R-squared: 0.9807, Adjusted R-squared: 0.9751
F-statistic: 176.1 on 15 and 52 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.394236688 0.788473376 0.6057633
[2,] 0.277816185 0.555632370 0.7221838
[3,] 0.196706716 0.393413432 0.8032933
[4,] 0.124075750 0.248151500 0.8759243
[5,] 0.134113604 0.268227208 0.8658864
[6,] 0.096349674 0.192699348 0.9036503
[7,] 0.053714232 0.107428464 0.9462858
[8,] 0.046554723 0.093109447 0.9534453
[9,] 0.037932369 0.075864738 0.9620676
[10,] 0.035074419 0.070148838 0.9649256
[11,] 0.026144628 0.052289255 0.9738554
[12,] 0.028803861 0.057607721 0.9711961
[13,] 0.016726653 0.033453306 0.9832733
[14,] 0.009317009 0.018634019 0.9906830
[15,] 0.005860576 0.011721152 0.9941394
[16,] 0.003396607 0.006793213 0.9966034
[17,] 0.002056732 0.004113465 0.9979433
[18,] 0.002035185 0.004070370 0.9979648
[19,] 0.010793936 0.021587872 0.9892061
[20,] 0.012320432 0.024640864 0.9876796
[21,] 0.008401169 0.016802338 0.9915988
[22,] 0.005019935 0.010039870 0.9949801
[23,] 0.134450778 0.268901556 0.8655492
[24,] 0.325274229 0.650548458 0.6747258
[25,] 0.258524734 0.517049468 0.7414753
[26,] 0.315604327 0.631208654 0.6843957
[27,] 0.271008938 0.542017876 0.7289911
[28,] 0.215372496 0.430744992 0.7846275
[29,] 0.151371393 0.302742787 0.8486286
[30,] 0.179832646 0.359665291 0.8201674
[31,] 0.888596439 0.222807122 0.1114036
> postscript(file="/var/www/html/rcomp/tmp/152yf1258744037.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/2v3ta1258744037.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/398ny1258744037.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/4wl971258744037.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/5lu4p1258744037.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 68
Frequency = 1
1 2 3 4 5 6
-3547.83012 -778.80644 4167.02946 -4055.70344 960.73239 -2319.51138
7 8 9 10 11 12
-6767.36870 -360.28723 2440.68957 655.10874 2253.04353 -7588.41194
13 14 15 16 17 18
473.33695 2549.68250 -4212.52818 6144.59672 7942.44362 11528.39645
19 20 21 22 23 24
-2621.12198 -2827.60247 -6129.93117 3269.31532 4138.11612 -4689.21012
25 26 27 28 29 30
4674.11089 336.09777 -1022.19148 -921.54008 7305.15952 -45.71066
31 32 33 34 35 36
253.52570 -1086.74381 -3442.60776 3531.98699 2516.21177 3606.99473
37 38 39 40 41 42
11046.20250 1825.38012 -2817.53405 251.97652 4704.78138 -14589.57595
43 44 45 46 47 48
-3156.96075 -9702.03173 4430.23517 -4442.83113 -4450.96720 3595.43019
49 50 51 52 53 54
-1438.22507 -10094.85356 7103.92760 3122.43455 -12103.90064 2824.41073
55 56 57 58 59 60
228.79034 1991.10073 2701.61420 -3013.57993 -4456.40421 5075.19715
61 62 63 64 65 66
-11207.59516 6162.49961 -3218.70335 -4541.76426 -8809.21626 2601.99081
67 68
12063.13539 11985.56449
> postscript(file="/var/www/html/rcomp/tmp/62t1q1258744037.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 -3547.83012 NA
1 -778.80644 -3547.83012
2 4167.02946 -778.80644
3 -4055.70344 4167.02946
4 960.73239 -4055.70344
5 -2319.51138 960.73239
6 -6767.36870 -2319.51138
7 -360.28723 -6767.36870
8 2440.68957 -360.28723
9 655.10874 2440.68957
10 2253.04353 655.10874
11 -7588.41194 2253.04353
12 473.33695 -7588.41194
13 2549.68250 473.33695
14 -4212.52818 2549.68250
15 6144.59672 -4212.52818
16 7942.44362 6144.59672
17 11528.39645 7942.44362
18 -2621.12198 11528.39645
19 -2827.60247 -2621.12198
20 -6129.93117 -2827.60247
21 3269.31532 -6129.93117
22 4138.11612 3269.31532
23 -4689.21012 4138.11612
24 4674.11089 -4689.21012
25 336.09777 4674.11089
26 -1022.19148 336.09777
27 -921.54008 -1022.19148
28 7305.15952 -921.54008
29 -45.71066 7305.15952
30 253.52570 -45.71066
31 -1086.74381 253.52570
32 -3442.60776 -1086.74381
33 3531.98699 -3442.60776
34 2516.21177 3531.98699
35 3606.99473 2516.21177
36 11046.20250 3606.99473
37 1825.38012 11046.20250
38 -2817.53405 1825.38012
39 251.97652 -2817.53405
40 4704.78138 251.97652
41 -14589.57595 4704.78138
42 -3156.96075 -14589.57595
43 -9702.03173 -3156.96075
44 4430.23517 -9702.03173
45 -4442.83113 4430.23517
46 -4450.96720 -4442.83113
47 3595.43019 -4450.96720
48 -1438.22507 3595.43019
49 -10094.85356 -1438.22507
50 7103.92760 -10094.85356
51 3122.43455 7103.92760
52 -12103.90064 3122.43455
53 2824.41073 -12103.90064
54 228.79034 2824.41073
55 1991.10073 228.79034
56 2701.61420 1991.10073
57 -3013.57993 2701.61420
58 -4456.40421 -3013.57993
59 5075.19715 -4456.40421
60 -11207.59516 5075.19715
61 6162.49961 -11207.59516
62 -3218.70335 6162.49961
63 -4541.76426 -3218.70335
64 -8809.21626 -4541.76426
65 2601.99081 -8809.21626
66 12063.13539 2601.99081
67 11985.56449 12063.13539
68 NA 11985.56449
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -778.80644 -3547.83012
[2,] 4167.02946 -778.80644
[3,] -4055.70344 4167.02946
[4,] 960.73239 -4055.70344
[5,] -2319.51138 960.73239
[6,] -6767.36870 -2319.51138
[7,] -360.28723 -6767.36870
[8,] 2440.68957 -360.28723
[9,] 655.10874 2440.68957
[10,] 2253.04353 655.10874
[11,] -7588.41194 2253.04353
[12,] 473.33695 -7588.41194
[13,] 2549.68250 473.33695
[14,] -4212.52818 2549.68250
[15,] 6144.59672 -4212.52818
[16,] 7942.44362 6144.59672
[17,] 11528.39645 7942.44362
[18,] -2621.12198 11528.39645
[19,] -2827.60247 -2621.12198
[20,] -6129.93117 -2827.60247
[21,] 3269.31532 -6129.93117
[22,] 4138.11612 3269.31532
[23,] -4689.21012 4138.11612
[24,] 4674.11089 -4689.21012
[25,] 336.09777 4674.11089
[26,] -1022.19148 336.09777
[27,] -921.54008 -1022.19148
[28,] 7305.15952 -921.54008
[29,] -45.71066 7305.15952
[30,] 253.52570 -45.71066
[31,] -1086.74381 253.52570
[32,] -3442.60776 -1086.74381
[33,] 3531.98699 -3442.60776
[34,] 2516.21177 3531.98699
[35,] 3606.99473 2516.21177
[36,] 11046.20250 3606.99473
[37,] 1825.38012 11046.20250
[38,] -2817.53405 1825.38012
[39,] 251.97652 -2817.53405
[40,] 4704.78138 251.97652
[41,] -14589.57595 4704.78138
[42,] -3156.96075 -14589.57595
[43,] -9702.03173 -3156.96075
[44,] 4430.23517 -9702.03173
[45,] -4442.83113 4430.23517
[46,] -4450.96720 -4442.83113
[47,] 3595.43019 -4450.96720
[48,] -1438.22507 3595.43019
[49,] -10094.85356 -1438.22507
[50,] 7103.92760 -10094.85356
[51,] 3122.43455 7103.92760
[52,] -12103.90064 3122.43455
[53,] 2824.41073 -12103.90064
[54,] 228.79034 2824.41073
[55,] 1991.10073 228.79034
[56,] 2701.61420 1991.10073
[57,] -3013.57993 2701.61420
[58,] -4456.40421 -3013.57993
[59,] 5075.19715 -4456.40421
[60,] -11207.59516 5075.19715
[61,] 6162.49961 -11207.59516
[62,] -3218.70335 6162.49961
[63,] -4541.76426 -3218.70335
[64,] -8809.21626 -4541.76426
[65,] 2601.99081 -8809.21626
[66,] 12063.13539 2601.99081
[67,] 11985.56449 12063.13539
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -778.80644 -3547.83012
2 4167.02946 -778.80644
3 -4055.70344 4167.02946
4 960.73239 -4055.70344
5 -2319.51138 960.73239
6 -6767.36870 -2319.51138
7 -360.28723 -6767.36870
8 2440.68957 -360.28723
9 655.10874 2440.68957
10 2253.04353 655.10874
11 -7588.41194 2253.04353
12 473.33695 -7588.41194
13 2549.68250 473.33695
14 -4212.52818 2549.68250
15 6144.59672 -4212.52818
16 7942.44362 6144.59672
17 11528.39645 7942.44362
18 -2621.12198 11528.39645
19 -2827.60247 -2621.12198
20 -6129.93117 -2827.60247
21 3269.31532 -6129.93117
22 4138.11612 3269.31532
23 -4689.21012 4138.11612
24 4674.11089 -4689.21012
25 336.09777 4674.11089
26 -1022.19148 336.09777
27 -921.54008 -1022.19148
28 7305.15952 -921.54008
29 -45.71066 7305.15952
30 253.52570 -45.71066
31 -1086.74381 253.52570
32 -3442.60776 -1086.74381
33 3531.98699 -3442.60776
34 2516.21177 3531.98699
35 3606.99473 2516.21177
36 11046.20250 3606.99473
37 1825.38012 11046.20250
38 -2817.53405 1825.38012
39 251.97652 -2817.53405
40 4704.78138 251.97652
41 -14589.57595 4704.78138
42 -3156.96075 -14589.57595
43 -9702.03173 -3156.96075
44 4430.23517 -9702.03173
45 -4442.83113 4430.23517
46 -4450.96720 -4442.83113
47 3595.43019 -4450.96720
48 -1438.22507 3595.43019
49 -10094.85356 -1438.22507
50 7103.92760 -10094.85356
51 3122.43455 7103.92760
52 -12103.90064 3122.43455
53 2824.41073 -12103.90064
54 228.79034 2824.41073
55 1991.10073 228.79034
56 2701.61420 1991.10073
57 -3013.57993 2701.61420
58 -4456.40421 -3013.57993
59 5075.19715 -4456.40421
60 -11207.59516 5075.19715
61 6162.49961 -11207.59516
62 -3218.70335 6162.49961
63 -4541.76426 -3218.70335
64 -8809.21626 -4541.76426
65 2601.99081 -8809.21626
66 12063.13539 2601.99081
67 11985.56449 12063.13539
> 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/72nu91258744037.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/8krnd1258744037.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/9t96g1258744037.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/10n73g1258744037.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/11naa41258744037.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/12jy171258744037.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/13cvmp1258744037.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/14cicb1258744037.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/15p9vx1258744037.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/16yca11258744037.tab")
+ }
>
> system("convert tmp/152yf1258744037.ps tmp/152yf1258744037.png")
> system("convert tmp/2v3ta1258744037.ps tmp/2v3ta1258744037.png")
> system("convert tmp/398ny1258744037.ps tmp/398ny1258744037.png")
> system("convert tmp/4wl971258744037.ps tmp/4wl971258744037.png")
> system("convert tmp/5lu4p1258744037.ps tmp/5lu4p1258744037.png")
> system("convert tmp/62t1q1258744037.ps tmp/62t1q1258744037.png")
> system("convert tmp/72nu91258744037.ps tmp/72nu91258744037.png")
> system("convert tmp/8krnd1258744037.ps tmp/8krnd1258744037.png")
> system("convert tmp/9t96g1258744037.ps tmp/9t96g1258744037.png")
> system("convert tmp/10n73g1258744037.ps tmp/10n73g1258744037.png")
>
>
> proc.time()
user system elapsed
2.566 1.575 2.934