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(0.7461,0.527,0.7775,0.472,0.7790,0,0.7744,0.052,0.7905,0.313,0.7719,0.364,0.7811,0.363,0.7557,-0.155,0.7637,0.052,0.7595,0.568,0.7471,0.668,0.7615,1.378,0.7487,0.252,0.7389,-0.402,0.7337,-0.05,0.7510,0.555,0.7382,0.05,0.7159,0.15,0.7542,0.45,0.7636,0.299,0.7433,0.199,0.7658,0.496,0.7627,0.444,0.7480,-0.393,0.7692,-0.444,0.7850,0.198,0.7913,0.494,0.7720,0.133,0.7880,0.388,0.8070,0.484,0.8268,0.278,0.8244,0.369,0.8487,0.165,0.8572,0.155,0.8214,0.087,0.8827,0.414,0.9216,0.36,0.8865,0.975,0.8816,0.27,0.8884,0.359,0.9466,0.169,0.9180,0.381,0.9337,0.154,0.9559,0.486,0.9626,0.925,0.9434,0.728,0.8639,-0.014,0.7996,0.046,0.6680,-0.819,0.6572,-1.674,0.6928,-0.788,0.6438,0.279,0.6454,0.396,0.6873,-0.141,0.7265,-0.019,0.7912,0.099,0.8114,0.742,0.8281,0.005,0.8393,0.448),dim=c(2,59),dimnames=list(c('USDOLLAR','Amerikaanse_inflatie'),1:59))
> y <- array(NA,dim=c(2,59),dimnames=list(c('USDOLLAR','Amerikaanse_inflatie'),1:59))
> 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
USDOLLAR Amerikaanse_inflatie M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 0.7461 0.527 1 0 0 0 0 0 0 0 0 0 0 1
2 0.7775 0.472 0 1 0 0 0 0 0 0 0 0 0 2
3 0.7790 0.000 0 0 1 0 0 0 0 0 0 0 0 3
4 0.7744 0.052 0 0 0 1 0 0 0 0 0 0 0 4
5 0.7905 0.313 0 0 0 0 1 0 0 0 0 0 0 5
6 0.7719 0.364 0 0 0 0 0 1 0 0 0 0 0 6
7 0.7811 0.363 0 0 0 0 0 0 1 0 0 0 0 7
8 0.7557 -0.155 0 0 0 0 0 0 0 1 0 0 0 8
9 0.7637 0.052 0 0 0 0 0 0 0 0 1 0 0 9
10 0.7595 0.568 0 0 0 0 0 0 0 0 0 1 0 10
11 0.7471 0.668 0 0 0 0 0 0 0 0 0 0 1 11
12 0.7615 1.378 0 0 0 0 0 0 0 0 0 0 0 12
13 0.7487 0.252 1 0 0 0 0 0 0 0 0 0 0 13
14 0.7389 -0.402 0 1 0 0 0 0 0 0 0 0 0 14
15 0.7337 -0.050 0 0 1 0 0 0 0 0 0 0 0 15
16 0.7510 0.555 0 0 0 1 0 0 0 0 0 0 0 16
17 0.7382 0.050 0 0 0 0 1 0 0 0 0 0 0 17
18 0.7159 0.150 0 0 0 0 0 1 0 0 0 0 0 18
19 0.7542 0.450 0 0 0 0 0 0 1 0 0 0 0 19
20 0.7636 0.299 0 0 0 0 0 0 0 1 0 0 0 20
21 0.7433 0.199 0 0 0 0 0 0 0 0 1 0 0 21
22 0.7658 0.496 0 0 0 0 0 0 0 0 0 1 0 22
23 0.7627 0.444 0 0 0 0 0 0 0 0 0 0 1 23
24 0.7480 -0.393 0 0 0 0 0 0 0 0 0 0 0 24
25 0.7692 -0.444 1 0 0 0 0 0 0 0 0 0 0 25
26 0.7850 0.198 0 1 0 0 0 0 0 0 0 0 0 26
27 0.7913 0.494 0 0 1 0 0 0 0 0 0 0 0 27
28 0.7720 0.133 0 0 0 1 0 0 0 0 0 0 0 28
29 0.7880 0.388 0 0 0 0 1 0 0 0 0 0 0 29
30 0.8070 0.484 0 0 0 0 0 1 0 0 0 0 0 30
31 0.8268 0.278 0 0 0 0 0 0 1 0 0 0 0 31
32 0.8244 0.369 0 0 0 0 0 0 0 1 0 0 0 32
33 0.8487 0.165 0 0 0 0 0 0 0 0 1 0 0 33
34 0.8572 0.155 0 0 0 0 0 0 0 0 0 1 0 34
35 0.8214 0.087 0 0 0 0 0 0 0 0 0 0 1 35
36 0.8827 0.414 0 0 0 0 0 0 0 0 0 0 0 36
37 0.9216 0.360 1 0 0 0 0 0 0 0 0 0 0 37
38 0.8865 0.975 0 1 0 0 0 0 0 0 0 0 0 38
39 0.8816 0.270 0 0 1 0 0 0 0 0 0 0 0 39
40 0.8884 0.359 0 0 0 1 0 0 0 0 0 0 0 40
41 0.9466 0.169 0 0 0 0 1 0 0 0 0 0 0 41
42 0.9180 0.381 0 0 0 0 0 1 0 0 0 0 0 42
43 0.9337 0.154 0 0 0 0 0 0 1 0 0 0 0 43
44 0.9559 0.486 0 0 0 0 0 0 0 1 0 0 0 44
45 0.9626 0.925 0 0 0 0 0 0 0 0 1 0 0 45
46 0.9434 0.728 0 0 0 0 0 0 0 0 0 1 0 46
47 0.8639 -0.014 0 0 0 0 0 0 0 0 0 0 1 47
48 0.7996 0.046 0 0 0 0 0 0 0 0 0 0 0 48
49 0.6680 -0.819 1 0 0 0 0 0 0 0 0 0 0 49
50 0.6572 -1.674 0 1 0 0 0 0 0 0 0 0 0 50
51 0.6928 -0.788 0 0 1 0 0 0 0 0 0 0 0 51
52 0.6438 0.279 0 0 0 1 0 0 0 0 0 0 0 52
53 0.6454 0.396 0 0 0 0 1 0 0 0 0 0 0 53
54 0.6873 -0.141 0 0 0 0 0 1 0 0 0 0 0 54
55 0.7265 -0.019 0 0 0 0 0 0 1 0 0 0 0 55
56 0.7912 0.099 0 0 0 0 0 0 0 1 0 0 0 56
57 0.8114 0.742 0 0 0 0 0 0 0 0 1 0 0 57
58 0.8281 0.005 0 0 0 0 0 0 0 0 0 1 0 58
59 0.8393 0.448 0 0 0 0 0 0 0 0 0 0 1 59
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Amerikaanse_inflatie M1
0.728301 0.083321 0.011528
M2 M3 M4
0.013625 0.013018 -0.022257
M5 M6 M7
-0.006722 -0.008461 0.014861
M8 M9 M10
0.029376 0.019423 0.025148
M11 t
0.005225 0.001318
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.179045 -0.042762 -0.002228 0.035367 0.156889
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.7283005 0.0424831 17.143 < 2e-16 ***
Amerikaanse_inflatie 0.0833205 0.0236186 3.528 0.000978 ***
M1 0.0115276 0.0502517 0.229 0.819600
M2 0.0136251 0.0504943 0.270 0.788521
M3 0.0130177 0.0500982 0.260 0.796170
M4 -0.0222569 0.0492608 -0.452 0.653570
M5 -0.0067221 0.0492579 -0.136 0.892061
M6 -0.0084606 0.0492662 -0.172 0.864417
M7 0.0148610 0.0492646 0.302 0.764303
M8 0.0293757 0.0493006 0.596 0.554262
M9 0.0194232 0.0492538 0.394 0.695186
M10 0.0251479 0.0492624 0.510 0.612205
M11 0.0052254 0.0492755 0.106 0.916018
t 0.0013183 0.0005887 2.239 0.030122 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07333 on 45 degrees of freedom
Multiple R-squared: 0.3063, Adjusted R-squared: 0.1059
F-statistic: 1.528 on 13 and 45 DF, p-value: 0.1444
> 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,] 1.214201e-02 2.428402e-02 0.9878580
[2,] 3.239141e-03 6.478282e-03 0.9967609
[3,] 6.035481e-04 1.207096e-03 0.9993965
[4,] 2.541450e-04 5.082901e-04 0.9997459
[5,] 4.767509e-05 9.535019e-05 0.9999523
[6,] 3.252840e-05 6.505680e-05 0.9999675
[7,] 3.585352e-05 7.170704e-05 0.9999641
[8,] 1.141273e-05 2.282545e-05 0.9999886
[9,] 1.551921e-05 3.103841e-05 0.9999845
[10,] 1.418208e-05 2.836417e-05 0.9999858
[11,] 1.221640e-05 2.443280e-05 0.9999878
[12,] 3.447652e-06 6.895303e-06 0.9999966
[13,] 1.329694e-06 2.659388e-06 0.9999987
[14,] 2.542078e-06 5.084157e-06 0.9999975
[15,] 3.608352e-06 7.216703e-06 0.9999964
[16,] 7.279329e-06 1.455866e-05 0.9999927
[17,] 2.564379e-05 5.128759e-05 0.9999744
[18,] 2.218660e-04 4.437319e-04 0.9997781
[19,] 1.322904e-02 2.645808e-02 0.9867710
[20,] 4.527694e-02 9.055389e-02 0.9547231
[21,] 8.322934e-02 1.664587e-01 0.9167707
[22,] 5.041776e-02 1.008355e-01 0.9495822
[23,] 3.025251e-02 6.050503e-02 0.9697475
[24,] 2.677869e-02 5.355738e-02 0.9732213
[25,] 3.277239e-01 6.554477e-01 0.6722761
[26,] 3.186307e-01 6.372614e-01 0.6813693
> postscript(file="/var/www/html/rcomp/tmp/1mq0p1258652008.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/25e0w1258652008.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/3z66e1258652008.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/424rq1258652008.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/5z04o1258652008.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 = 59
Frequency = 1
1 2 3 4 5
-0.0389563091 -0.0063895605 0.0337268177 0.0587504331 0.0362505990
6 7 8 9 10
0.0138214514 -0.0015351973 0.0003918337 -0.0002213716 -0.0544577652
11 12 13 14 15
-0.0565856687 -0.0974361819 -0.0292631430 0.0120126055 -0.0232271364
16 17 18 19 20
-0.0223797764 -0.0099560813 -0.0401679351 -0.0515040646 -0.0453556696
21 22 23 24 25
-0.0486894710 -0.0579786676 -0.0381418499 0.0208045033 0.0334079679
26 27 28 29 30
-0.0076996958 -0.0267734878 0.0179615082 -0.0041384027 0.0072830257
31 32 33 34 35
0.0196070865 -0.0062080878 0.0437234464 0.0460136537 0.0344836000
36 37 38 39 40
0.0724448514 0.1029982777 0.0132402684 0.0663703311 0.0997110866
41 42 43 44 45
0.1568888135 0.1110450599 0.1210188519 0.0997234289 0.0784798596
46 47 48 49 50
0.0686510069 0.0695789931 0.0041868272 -0.0681867934 -0.0111636175
51 52 53 54 55
-0.0500965246 -0.1540432514 -0.1790449286 -0.0919816021 -0.0875866765
56 57 58 59
-0.0485515052 -0.0732924634 -0.0022282278 -0.0093350745
> postscript(file="/var/www/html/rcomp/tmp/6fe2c1258652008.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 = 59
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.0389563091 NA
1 -0.0063895605 -0.0389563091
2 0.0337268177 -0.0063895605
3 0.0587504331 0.0337268177
4 0.0362505990 0.0587504331
5 0.0138214514 0.0362505990
6 -0.0015351973 0.0138214514
7 0.0003918337 -0.0015351973
8 -0.0002213716 0.0003918337
9 -0.0544577652 -0.0002213716
10 -0.0565856687 -0.0544577652
11 -0.0974361819 -0.0565856687
12 -0.0292631430 -0.0974361819
13 0.0120126055 -0.0292631430
14 -0.0232271364 0.0120126055
15 -0.0223797764 -0.0232271364
16 -0.0099560813 -0.0223797764
17 -0.0401679351 -0.0099560813
18 -0.0515040646 -0.0401679351
19 -0.0453556696 -0.0515040646
20 -0.0486894710 -0.0453556696
21 -0.0579786676 -0.0486894710
22 -0.0381418499 -0.0579786676
23 0.0208045033 -0.0381418499
24 0.0334079679 0.0208045033
25 -0.0076996958 0.0334079679
26 -0.0267734878 -0.0076996958
27 0.0179615082 -0.0267734878
28 -0.0041384027 0.0179615082
29 0.0072830257 -0.0041384027
30 0.0196070865 0.0072830257
31 -0.0062080878 0.0196070865
32 0.0437234464 -0.0062080878
33 0.0460136537 0.0437234464
34 0.0344836000 0.0460136537
35 0.0724448514 0.0344836000
36 0.1029982777 0.0724448514
37 0.0132402684 0.1029982777
38 0.0663703311 0.0132402684
39 0.0997110866 0.0663703311
40 0.1568888135 0.0997110866
41 0.1110450599 0.1568888135
42 0.1210188519 0.1110450599
43 0.0997234289 0.1210188519
44 0.0784798596 0.0997234289
45 0.0686510069 0.0784798596
46 0.0695789931 0.0686510069
47 0.0041868272 0.0695789931
48 -0.0681867934 0.0041868272
49 -0.0111636175 -0.0681867934
50 -0.0500965246 -0.0111636175
51 -0.1540432514 -0.0500965246
52 -0.1790449286 -0.1540432514
53 -0.0919816021 -0.1790449286
54 -0.0875866765 -0.0919816021
55 -0.0485515052 -0.0875866765
56 -0.0732924634 -0.0485515052
57 -0.0022282278 -0.0732924634
58 -0.0093350745 -0.0022282278
59 NA -0.0093350745
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.0063895605 -0.0389563091
[2,] 0.0337268177 -0.0063895605
[3,] 0.0587504331 0.0337268177
[4,] 0.0362505990 0.0587504331
[5,] 0.0138214514 0.0362505990
[6,] -0.0015351973 0.0138214514
[7,] 0.0003918337 -0.0015351973
[8,] -0.0002213716 0.0003918337
[9,] -0.0544577652 -0.0002213716
[10,] -0.0565856687 -0.0544577652
[11,] -0.0974361819 -0.0565856687
[12,] -0.0292631430 -0.0974361819
[13,] 0.0120126055 -0.0292631430
[14,] -0.0232271364 0.0120126055
[15,] -0.0223797764 -0.0232271364
[16,] -0.0099560813 -0.0223797764
[17,] -0.0401679351 -0.0099560813
[18,] -0.0515040646 -0.0401679351
[19,] -0.0453556696 -0.0515040646
[20,] -0.0486894710 -0.0453556696
[21,] -0.0579786676 -0.0486894710
[22,] -0.0381418499 -0.0579786676
[23,] 0.0208045033 -0.0381418499
[24,] 0.0334079679 0.0208045033
[25,] -0.0076996958 0.0334079679
[26,] -0.0267734878 -0.0076996958
[27,] 0.0179615082 -0.0267734878
[28,] -0.0041384027 0.0179615082
[29,] 0.0072830257 -0.0041384027
[30,] 0.0196070865 0.0072830257
[31,] -0.0062080878 0.0196070865
[32,] 0.0437234464 -0.0062080878
[33,] 0.0460136537 0.0437234464
[34,] 0.0344836000 0.0460136537
[35,] 0.0724448514 0.0344836000
[36,] 0.1029982777 0.0724448514
[37,] 0.0132402684 0.1029982777
[38,] 0.0663703311 0.0132402684
[39,] 0.0997110866 0.0663703311
[40,] 0.1568888135 0.0997110866
[41,] 0.1110450599 0.1568888135
[42,] 0.1210188519 0.1110450599
[43,] 0.0997234289 0.1210188519
[44,] 0.0784798596 0.0997234289
[45,] 0.0686510069 0.0784798596
[46,] 0.0695789931 0.0686510069
[47,] 0.0041868272 0.0695789931
[48,] -0.0681867934 0.0041868272
[49,] -0.0111636175 -0.0681867934
[50,] -0.0500965246 -0.0111636175
[51,] -0.1540432514 -0.0500965246
[52,] -0.1790449286 -0.1540432514
[53,] -0.0919816021 -0.1790449286
[54,] -0.0875866765 -0.0919816021
[55,] -0.0485515052 -0.0875866765
[56,] -0.0732924634 -0.0485515052
[57,] -0.0022282278 -0.0732924634
[58,] -0.0093350745 -0.0022282278
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.0063895605 -0.0389563091
2 0.0337268177 -0.0063895605
3 0.0587504331 0.0337268177
4 0.0362505990 0.0587504331
5 0.0138214514 0.0362505990
6 -0.0015351973 0.0138214514
7 0.0003918337 -0.0015351973
8 -0.0002213716 0.0003918337
9 -0.0544577652 -0.0002213716
10 -0.0565856687 -0.0544577652
11 -0.0974361819 -0.0565856687
12 -0.0292631430 -0.0974361819
13 0.0120126055 -0.0292631430
14 -0.0232271364 0.0120126055
15 -0.0223797764 -0.0232271364
16 -0.0099560813 -0.0223797764
17 -0.0401679351 -0.0099560813
18 -0.0515040646 -0.0401679351
19 -0.0453556696 -0.0515040646
20 -0.0486894710 -0.0453556696
21 -0.0579786676 -0.0486894710
22 -0.0381418499 -0.0579786676
23 0.0208045033 -0.0381418499
24 0.0334079679 0.0208045033
25 -0.0076996958 0.0334079679
26 -0.0267734878 -0.0076996958
27 0.0179615082 -0.0267734878
28 -0.0041384027 0.0179615082
29 0.0072830257 -0.0041384027
30 0.0196070865 0.0072830257
31 -0.0062080878 0.0196070865
32 0.0437234464 -0.0062080878
33 0.0460136537 0.0437234464
34 0.0344836000 0.0460136537
35 0.0724448514 0.0344836000
36 0.1029982777 0.0724448514
37 0.0132402684 0.1029982777
38 0.0663703311 0.0132402684
39 0.0997110866 0.0663703311
40 0.1568888135 0.0997110866
41 0.1110450599 0.1568888135
42 0.1210188519 0.1110450599
43 0.0997234289 0.1210188519
44 0.0784798596 0.0997234289
45 0.0686510069 0.0784798596
46 0.0695789931 0.0686510069
47 0.0041868272 0.0695789931
48 -0.0681867934 0.0041868272
49 -0.0111636175 -0.0681867934
50 -0.0500965246 -0.0111636175
51 -0.1540432514 -0.0500965246
52 -0.1790449286 -0.1540432514
53 -0.0919816021 -0.1790449286
54 -0.0875866765 -0.0919816021
55 -0.0485515052 -0.0875866765
56 -0.0732924634 -0.0485515052
57 -0.0022282278 -0.0732924634
58 -0.0093350745 -0.0022282278
> 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/7me6w1258652008.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/81u9i1258652008.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/9t68q1258652008.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/107ush1258652008.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/11iocj1258652008.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/12da451258652008.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/13o93e1258652008.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/146jxw1258652008.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/15lhj41258652008.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/1699eg1258652009.tab")
+ }
>
> system("convert tmp/1mq0p1258652008.ps tmp/1mq0p1258652008.png")
> system("convert tmp/25e0w1258652008.ps tmp/25e0w1258652008.png")
> system("convert tmp/3z66e1258652008.ps tmp/3z66e1258652008.png")
> system("convert tmp/424rq1258652008.ps tmp/424rq1258652008.png")
> system("convert tmp/5z04o1258652008.ps tmp/5z04o1258652008.png")
> system("convert tmp/6fe2c1258652008.ps tmp/6fe2c1258652008.png")
> system("convert tmp/7me6w1258652008.ps tmp/7me6w1258652008.png")
> system("convert tmp/81u9i1258652008.ps tmp/81u9i1258652008.png")
> system("convert tmp/9t68q1258652008.ps tmp/9t68q1258652008.png")
> system("convert tmp/107ush1258652008.ps tmp/107ush1258652008.png")
>
>
> proc.time()
user system elapsed
2.427 1.580 3.098