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(467,98.8,460,100.5,448,110.4,443,96.4,436,101.9,431,106.2,484,81,510,94.7,513,101,503,109.4,471,102.3,471,90.7,476,96.2,475,96.1,470,106,461,103.1,455,102,456,104.7,517,86,525,92.1,523,106.9,519,112.6,509,101.7,512,92,519,97.4,517,97,510,105.4,509,102.7,501,98.1,507,104.5,569,87.4,580,89.9,578,109.8,565,111.7,547,98.6,555,96.9,562,95.1,561,97,555,112.7,544,102.9,537,97.4,543,111.4,594,87.4,611,96.8,613,114.1,611,110.3,594,103.9,595,101.6,591,94.6,589,95.9,584,104.7,573,102.8,567,98.1,569,113.9,621,80.9,629,95.7,628,113.2,612,105.9,595,108.8,597,102.3,593,99,590,100.7,580,115.5,574,100.7,573,109.9,573,114.6,620,85.4,626,100.5,620,114.8,588,116.5,566,112.9,557,102),dim=c(2,72),dimnames=list(c('Y','X'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Y','X'),1:72))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 467 98.8 1 0 0 0 0 0 0 0 0 0 0
2 460 100.5 0 1 0 0 0 0 0 0 0 0 0
3 448 110.4 0 0 1 0 0 0 0 0 0 0 0
4 443 96.4 0 0 0 1 0 0 0 0 0 0 0
5 436 101.9 0 0 0 0 1 0 0 0 0 0 0
6 431 106.2 0 0 0 0 0 1 0 0 0 0 0
7 484 81.0 0 0 0 0 0 0 1 0 0 0 0
8 510 94.7 0 0 0 0 0 0 0 1 0 0 0
9 513 101.0 0 0 0 0 0 0 0 0 1 0 0
10 503 109.4 0 0 0 0 0 0 0 0 0 1 0
11 471 102.3 0 0 0 0 0 0 0 0 0 0 1
12 471 90.7 0 0 0 0 0 0 0 0 0 0 0
13 476 96.2 1 0 0 0 0 0 0 0 0 0 0
14 475 96.1 0 1 0 0 0 0 0 0 0 0 0
15 470 106.0 0 0 1 0 0 0 0 0 0 0 0
16 461 103.1 0 0 0 1 0 0 0 0 0 0 0
17 455 102.0 0 0 0 0 1 0 0 0 0 0 0
18 456 104.7 0 0 0 0 0 1 0 0 0 0 0
19 517 86.0 0 0 0 0 0 0 1 0 0 0 0
20 525 92.1 0 0 0 0 0 0 0 1 0 0 0
21 523 106.9 0 0 0 0 0 0 0 0 1 0 0
22 519 112.6 0 0 0 0 0 0 0 0 0 1 0
23 509 101.7 0 0 0 0 0 0 0 0 0 0 1
24 512 92.0 0 0 0 0 0 0 0 0 0 0 0
25 519 97.4 1 0 0 0 0 0 0 0 0 0 0
26 517 97.0 0 1 0 0 0 0 0 0 0 0 0
27 510 105.4 0 0 1 0 0 0 0 0 0 0 0
28 509 102.7 0 0 0 1 0 0 0 0 0 0 0
29 501 98.1 0 0 0 0 1 0 0 0 0 0 0
30 507 104.5 0 0 0 0 0 1 0 0 0 0 0
31 569 87.4 0 0 0 0 0 0 1 0 0 0 0
32 580 89.9 0 0 0 0 0 0 0 1 0 0 0
33 578 109.8 0 0 0 0 0 0 0 0 1 0 0
34 565 111.7 0 0 0 0 0 0 0 0 0 1 0
35 547 98.6 0 0 0 0 0 0 0 0 0 0 1
36 555 96.9 0 0 0 0 0 0 0 0 0 0 0
37 562 95.1 1 0 0 0 0 0 0 0 0 0 0
38 561 97.0 0 1 0 0 0 0 0 0 0 0 0
39 555 112.7 0 0 1 0 0 0 0 0 0 0 0
40 544 102.9 0 0 0 1 0 0 0 0 0 0 0
41 537 97.4 0 0 0 0 1 0 0 0 0 0 0
42 543 111.4 0 0 0 0 0 1 0 0 0 0 0
43 594 87.4 0 0 0 0 0 0 1 0 0 0 0
44 611 96.8 0 0 0 0 0 0 0 1 0 0 0
45 613 114.1 0 0 0 0 0 0 0 0 1 0 0
46 611 110.3 0 0 0 0 0 0 0 0 0 1 0
47 594 103.9 0 0 0 0 0 0 0 0 0 0 1
48 595 101.6 0 0 0 0 0 0 0 0 0 0 0
49 591 94.6 1 0 0 0 0 0 0 0 0 0 0
50 589 95.9 0 1 0 0 0 0 0 0 0 0 0
51 584 104.7 0 0 1 0 0 0 0 0 0 0 0
52 573 102.8 0 0 0 1 0 0 0 0 0 0 0
53 567 98.1 0 0 0 0 1 0 0 0 0 0 0
54 569 113.9 0 0 0 0 0 1 0 0 0 0 0
55 621 80.9 0 0 0 0 0 0 1 0 0 0 0
56 629 95.7 0 0 0 0 0 0 0 1 0 0 0
57 628 113.2 0 0 0 0 0 0 0 0 1 0 0
58 612 105.9 0 0 0 0 0 0 0 0 0 1 0
59 595 108.8 0 0 0 0 0 0 0 0 0 0 1
60 597 102.3 0 0 0 0 0 0 0 0 0 0 0
61 593 99.0 1 0 0 0 0 0 0 0 0 0 0
62 590 100.7 0 1 0 0 0 0 0 0 0 0 0
63 580 115.5 0 0 1 0 0 0 0 0 0 0 0
64 574 100.7 0 0 0 1 0 0 0 0 0 0 0
65 573 109.9 0 0 0 0 1 0 0 0 0 0 0
66 573 114.6 0 0 0 0 0 1 0 0 0 0 0
67 620 85.4 0 0 0 0 0 0 1 0 0 0 0
68 626 100.5 0 0 0 0 0 0 0 1 0 0 0
69 620 114.8 0 0 0 0 0 0 0 0 1 0 0
70 588 116.5 0 0 0 0 0 0 0 0 0 1 0
71 566 112.9 0 0 0 0 0 0 0 0 0 0 1
72 557 102.0 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
4.052 5.572 -9.080 -17.412 -87.603 -51.954
M5 M6 M7 M8 M9 M10
-56.673 -99.493 91.552 47.008 -37.673 -56.636
M11
-40.491
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-86.67 -38.63 10.67 31.57 84.11
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.052 155.550 0.026 0.979305
X 5.572 1.581 3.526 0.000824 ***
M1 -9.080 28.607 -0.317 0.752049
M2 -17.412 28.587 -0.609 0.544794
M3 -87.603 33.901 -2.584 0.012259 *
M4 -51.954 29.224 -1.778 0.080589 .
M5 -56.673 29.159 -1.944 0.056724 .
M6 -99.493 33.986 -2.927 0.004848 **
M7 91.552 35.110 2.608 0.011531 *
M8 47.008 28.885 1.627 0.108975
M9 -37.673 34.642 -1.087 0.281245
M10 -56.636 35.653 -1.589 0.117513
M11 -40.491 30.717 -1.318 0.192530
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 49.51 on 59 degrees of freedom
Multiple R-squared: 0.3328, Adjusted R-squared: 0.1971
F-statistic: 2.453 on 12 and 59 DF, p-value: 0.01147
> 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.05693426 1.138685e-01 9.430657e-01
[2,] 0.03466147 6.932295e-02 9.653385e-01
[3,] 0.02645807 5.291613e-02 9.735419e-01
[4,] 0.04116933 8.233866e-02 9.588307e-01
[5,] 0.02932895 5.865791e-02 9.706710e-01
[6,] 0.01841700 3.683400e-02 9.815830e-01
[7,] 0.01650998 3.301996e-02 9.834900e-01
[8,] 0.03464091 6.928181e-02 9.653591e-01
[9,] 0.05806534 1.161307e-01 9.419347e-01
[10,] 0.17598379 3.519676e-01 8.240162e-01
[11,] 0.37199330 7.439866e-01 6.280067e-01
[12,] 0.56194942 8.761012e-01 4.380506e-01
[13,] 0.77680012 4.463998e-01 2.231999e-01
[14,] 0.89630340 2.073932e-01 1.036966e-01
[15,] 0.95872991 8.254017e-02 4.127009e-02
[16,] 0.98942566 2.114868e-02 1.057434e-02
[17,] 0.99512150 9.757002e-03 4.878501e-03
[18,] 0.99811082 3.778354e-03 1.889177e-03
[19,] 0.99923066 1.538674e-03 7.693369e-04
[20,] 0.99968557 6.288514e-04 3.144257e-04
[21,] 0.99982445 3.510943e-04 1.755472e-04
[22,] 0.99994206 1.158859e-04 5.794293e-05
[23,] 0.99997712 4.576259e-05 2.288130e-05
[24,] 0.99998627 2.746167e-05 1.373084e-05
[25,] 0.99999252 1.495665e-05 7.478327e-06
[26,] 0.99999861 2.772557e-06 1.386279e-06
[27,] 0.99999948 1.035687e-06 5.178435e-07
[28,] 0.99999955 9.004418e-07 4.502209e-07
[29,] 0.99999935 1.308084e-06 6.540418e-07
[30,] 0.99999821 3.573606e-06 1.786803e-06
[31,] 0.99999670 6.597995e-06 3.298997e-06
[32,] 0.99999133 1.734470e-05 8.672352e-06
[33,] 0.99998358 3.284109e-05 1.642055e-05
[34,] 0.99996063 7.874130e-05 3.937065e-05
[35,] 0.99989383 2.123395e-04 1.061698e-04
[36,] 0.99972267 5.546647e-04 2.773324e-04
[37,] 0.99907241 1.855189e-03 9.275947e-04
[38,] 0.99848305 3.033903e-03 1.516951e-03
[39,] 0.99441933 1.116135e-02 5.580673e-03
[40,] 0.98333070 3.333861e-02 1.666930e-02
[41,] 0.94710776 1.057845e-01 5.289224e-02
> postscript(file="/var/www/html/rcomp/tmp/1mri81260633551.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/27pqo1260633551.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/3fbv81260633551.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/48dca1260633551.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/59jwu1260633551.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 = 72
Frequency = 1
1 2 3 4 5 6
-78.5330023 -86.6741968 -83.6513491 -46.2851850 -79.2149865 -65.3563526
7 8 9 10 11 12
-62.9746994 -68.7735467 -16.2000978 -54.0458670 -62.6260485 -38.4760974
13 14 15 16 17 18
-55.0445548 -47.1552857 -37.1324380 -65.6207997 -60.7722345 -31.9976329
19 20 21 22 23 24
-57.8370984 -39.2850993 -39.0777286 -55.8778024 -21.2825606 -4.7203211
25 26 27 28 29 30
-18.7315306 -10.1705175 6.2110499 -15.3918077 6.9604367 20.1168630
31 32 33 34 35 36
-13.6385701 27.9743563 -0.2379200 -4.8625705 33.9921267 10.9745279
37 38 39 40 41 42
37.0851730 33.8294825 10.5319474 18.4936963 46.8611725 17.6667524
43 44 45 46 47 48
11.3614299 20.5242457 10.8004168 48.9389012 51.4579838 24.7838728
49 50 51 52 53 54
68.8714129 67.9592103 84.1117858 48.0509443 72.9604367 29.7355530
55 56 57 58 59 60
74.5825486 44.6539735 30.8156487 74.4578123 25.1528328 22.8831370
61 62 63 64 65 66
46.3525018 42.2113072 19.9290040 60.7531518 13.2051751 29.8348171
67 68 69 70 71 72
48.5063895 14.9060705 13.8996810 -8.6104736 -26.6943343 -15.4451191
> postscript(file="/var/www/html/rcomp/tmp/6uzvx1260633551.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -78.5330023 NA
1 -86.6741968 -78.5330023
2 -83.6513491 -86.6741968
3 -46.2851850 -83.6513491
4 -79.2149865 -46.2851850
5 -65.3563526 -79.2149865
6 -62.9746994 -65.3563526
7 -68.7735467 -62.9746994
8 -16.2000978 -68.7735467
9 -54.0458670 -16.2000978
10 -62.6260485 -54.0458670
11 -38.4760974 -62.6260485
12 -55.0445548 -38.4760974
13 -47.1552857 -55.0445548
14 -37.1324380 -47.1552857
15 -65.6207997 -37.1324380
16 -60.7722345 -65.6207997
17 -31.9976329 -60.7722345
18 -57.8370984 -31.9976329
19 -39.2850993 -57.8370984
20 -39.0777286 -39.2850993
21 -55.8778024 -39.0777286
22 -21.2825606 -55.8778024
23 -4.7203211 -21.2825606
24 -18.7315306 -4.7203211
25 -10.1705175 -18.7315306
26 6.2110499 -10.1705175
27 -15.3918077 6.2110499
28 6.9604367 -15.3918077
29 20.1168630 6.9604367
30 -13.6385701 20.1168630
31 27.9743563 -13.6385701
32 -0.2379200 27.9743563
33 -4.8625705 -0.2379200
34 33.9921267 -4.8625705
35 10.9745279 33.9921267
36 37.0851730 10.9745279
37 33.8294825 37.0851730
38 10.5319474 33.8294825
39 18.4936963 10.5319474
40 46.8611725 18.4936963
41 17.6667524 46.8611725
42 11.3614299 17.6667524
43 20.5242457 11.3614299
44 10.8004168 20.5242457
45 48.9389012 10.8004168
46 51.4579838 48.9389012
47 24.7838728 51.4579838
48 68.8714129 24.7838728
49 67.9592103 68.8714129
50 84.1117858 67.9592103
51 48.0509443 84.1117858
52 72.9604367 48.0509443
53 29.7355530 72.9604367
54 74.5825486 29.7355530
55 44.6539735 74.5825486
56 30.8156487 44.6539735
57 74.4578123 30.8156487
58 25.1528328 74.4578123
59 22.8831370 25.1528328
60 46.3525018 22.8831370
61 42.2113072 46.3525018
62 19.9290040 42.2113072
63 60.7531518 19.9290040
64 13.2051751 60.7531518
65 29.8348171 13.2051751
66 48.5063895 29.8348171
67 14.9060705 48.5063895
68 13.8996810 14.9060705
69 -8.6104736 13.8996810
70 -26.6943343 -8.6104736
71 -15.4451191 -26.6943343
72 NA -15.4451191
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -86.6741968 -78.5330023
[2,] -83.6513491 -86.6741968
[3,] -46.2851850 -83.6513491
[4,] -79.2149865 -46.2851850
[5,] -65.3563526 -79.2149865
[6,] -62.9746994 -65.3563526
[7,] -68.7735467 -62.9746994
[8,] -16.2000978 -68.7735467
[9,] -54.0458670 -16.2000978
[10,] -62.6260485 -54.0458670
[11,] -38.4760974 -62.6260485
[12,] -55.0445548 -38.4760974
[13,] -47.1552857 -55.0445548
[14,] -37.1324380 -47.1552857
[15,] -65.6207997 -37.1324380
[16,] -60.7722345 -65.6207997
[17,] -31.9976329 -60.7722345
[18,] -57.8370984 -31.9976329
[19,] -39.2850993 -57.8370984
[20,] -39.0777286 -39.2850993
[21,] -55.8778024 -39.0777286
[22,] -21.2825606 -55.8778024
[23,] -4.7203211 -21.2825606
[24,] -18.7315306 -4.7203211
[25,] -10.1705175 -18.7315306
[26,] 6.2110499 -10.1705175
[27,] -15.3918077 6.2110499
[28,] 6.9604367 -15.3918077
[29,] 20.1168630 6.9604367
[30,] -13.6385701 20.1168630
[31,] 27.9743563 -13.6385701
[32,] -0.2379200 27.9743563
[33,] -4.8625705 -0.2379200
[34,] 33.9921267 -4.8625705
[35,] 10.9745279 33.9921267
[36,] 37.0851730 10.9745279
[37,] 33.8294825 37.0851730
[38,] 10.5319474 33.8294825
[39,] 18.4936963 10.5319474
[40,] 46.8611725 18.4936963
[41,] 17.6667524 46.8611725
[42,] 11.3614299 17.6667524
[43,] 20.5242457 11.3614299
[44,] 10.8004168 20.5242457
[45,] 48.9389012 10.8004168
[46,] 51.4579838 48.9389012
[47,] 24.7838728 51.4579838
[48,] 68.8714129 24.7838728
[49,] 67.9592103 68.8714129
[50,] 84.1117858 67.9592103
[51,] 48.0509443 84.1117858
[52,] 72.9604367 48.0509443
[53,] 29.7355530 72.9604367
[54,] 74.5825486 29.7355530
[55,] 44.6539735 74.5825486
[56,] 30.8156487 44.6539735
[57,] 74.4578123 30.8156487
[58,] 25.1528328 74.4578123
[59,] 22.8831370 25.1528328
[60,] 46.3525018 22.8831370
[61,] 42.2113072 46.3525018
[62,] 19.9290040 42.2113072
[63,] 60.7531518 19.9290040
[64,] 13.2051751 60.7531518
[65,] 29.8348171 13.2051751
[66,] 48.5063895 29.8348171
[67,] 14.9060705 48.5063895
[68,] 13.8996810 14.9060705
[69,] -8.6104736 13.8996810
[70,] -26.6943343 -8.6104736
[71,] -15.4451191 -26.6943343
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -86.6741968 -78.5330023
2 -83.6513491 -86.6741968
3 -46.2851850 -83.6513491
4 -79.2149865 -46.2851850
5 -65.3563526 -79.2149865
6 -62.9746994 -65.3563526
7 -68.7735467 -62.9746994
8 -16.2000978 -68.7735467
9 -54.0458670 -16.2000978
10 -62.6260485 -54.0458670
11 -38.4760974 -62.6260485
12 -55.0445548 -38.4760974
13 -47.1552857 -55.0445548
14 -37.1324380 -47.1552857
15 -65.6207997 -37.1324380
16 -60.7722345 -65.6207997
17 -31.9976329 -60.7722345
18 -57.8370984 -31.9976329
19 -39.2850993 -57.8370984
20 -39.0777286 -39.2850993
21 -55.8778024 -39.0777286
22 -21.2825606 -55.8778024
23 -4.7203211 -21.2825606
24 -18.7315306 -4.7203211
25 -10.1705175 -18.7315306
26 6.2110499 -10.1705175
27 -15.3918077 6.2110499
28 6.9604367 -15.3918077
29 20.1168630 6.9604367
30 -13.6385701 20.1168630
31 27.9743563 -13.6385701
32 -0.2379200 27.9743563
33 -4.8625705 -0.2379200
34 33.9921267 -4.8625705
35 10.9745279 33.9921267
36 37.0851730 10.9745279
37 33.8294825 37.0851730
38 10.5319474 33.8294825
39 18.4936963 10.5319474
40 46.8611725 18.4936963
41 17.6667524 46.8611725
42 11.3614299 17.6667524
43 20.5242457 11.3614299
44 10.8004168 20.5242457
45 48.9389012 10.8004168
46 51.4579838 48.9389012
47 24.7838728 51.4579838
48 68.8714129 24.7838728
49 67.9592103 68.8714129
50 84.1117858 67.9592103
51 48.0509443 84.1117858
52 72.9604367 48.0509443
53 29.7355530 72.9604367
54 74.5825486 29.7355530
55 44.6539735 74.5825486
56 30.8156487 44.6539735
57 74.4578123 30.8156487
58 25.1528328 74.4578123
59 22.8831370 25.1528328
60 46.3525018 22.8831370
61 42.2113072 46.3525018
62 19.9290040 42.2113072
63 60.7531518 19.9290040
64 13.2051751 60.7531518
65 29.8348171 13.2051751
66 48.5063895 29.8348171
67 14.9060705 48.5063895
68 13.8996810 14.9060705
69 -8.6104736 13.8996810
70 -26.6943343 -8.6104736
71 -15.4451191 -26.6943343
> 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/70s6t1260633551.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/8mf381260633551.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/9fspz1260633551.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/108wo21260633551.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/11qj821260633551.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/12wfte1260633551.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/13l9f21260633551.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/146z1n1260633551.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/150lae1260633551.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/16tkwn1260633551.tab")
+ }
>
> try(system("convert tmp/1mri81260633551.ps tmp/1mri81260633551.png",intern=TRUE))
character(0)
> try(system("convert tmp/27pqo1260633551.ps tmp/27pqo1260633551.png",intern=TRUE))
character(0)
> try(system("convert tmp/3fbv81260633551.ps tmp/3fbv81260633551.png",intern=TRUE))
character(0)
> try(system("convert tmp/48dca1260633551.ps tmp/48dca1260633551.png",intern=TRUE))
character(0)
> try(system("convert tmp/59jwu1260633551.ps tmp/59jwu1260633551.png",intern=TRUE))
character(0)
> try(system("convert tmp/6uzvx1260633551.ps tmp/6uzvx1260633551.png",intern=TRUE))
character(0)
> try(system("convert tmp/70s6t1260633551.ps tmp/70s6t1260633551.png",intern=TRUE))
character(0)
> try(system("convert tmp/8mf381260633551.ps tmp/8mf381260633551.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fspz1260633551.ps tmp/9fspz1260633551.png",intern=TRUE))
character(0)
> try(system("convert tmp/108wo21260633551.ps tmp/108wo21260633551.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.553 1.579 7.010