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(562,13.9,561,15.9,555,18.2,544,19.7,537,20.1,543,19.9,594,20,611,22.6,613,20.6,611,20.1,594,20.2,595,21.8,591,22,589,19.5,584,17.5,573,18.2,567,18.8,569,19.7,621,18.8,629,18.5,628,18.7,612,18.5,595,19.3,597,18.9,593,21.4,590,22.5,580,25,574,22.9,573,22.9,573,21.3,620,22.3,626,20.9,620,19.9,588,20.2,566,19.8,557,17.7,561,18.1,549,17.6,532,18.2,526,16,511,16.3,499,17.3,555,19,565,18.6,542,18,527,17.9,510,17.8,514,18.5,517,17.4,508,19,493,17.4,490,20.6,469,18.5,478,20,528,18.8,534,18.8,518,19.7,506,15.3,502,10.6,516,6.1,528,0.9),dim=c(2,61),dimnames=list(c('Y','X'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Y','X'),1:61))
> 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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 562 13.9 1 0 0 0 0 0 0 0 0 0 0 1
2 561 15.9 0 1 0 0 0 0 0 0 0 0 0 2
3 555 18.2 0 0 1 0 0 0 0 0 0 0 0 3
4 544 19.7 0 0 0 1 0 0 0 0 0 0 0 4
5 537 20.1 0 0 0 0 1 0 0 0 0 0 0 5
6 543 19.9 0 0 0 0 0 1 0 0 0 0 0 6
7 594 20.0 0 0 0 0 0 0 1 0 0 0 0 7
8 611 22.6 0 0 0 0 0 0 0 1 0 0 0 8
9 613 20.6 0 0 0 0 0 0 0 0 1 0 0 9
10 611 20.1 0 0 0 0 0 0 0 0 0 1 0 10
11 594 20.2 0 0 0 0 0 0 0 0 0 0 1 11
12 595 21.8 0 0 0 0 0 0 0 0 0 0 0 12
13 591 22.0 1 0 0 0 0 0 0 0 0 0 0 13
14 589 19.5 0 1 0 0 0 0 0 0 0 0 0 14
15 584 17.5 0 0 1 0 0 0 0 0 0 0 0 15
16 573 18.2 0 0 0 1 0 0 0 0 0 0 0 16
17 567 18.8 0 0 0 0 1 0 0 0 0 0 0 17
18 569 19.7 0 0 0 0 0 1 0 0 0 0 0 18
19 621 18.8 0 0 0 0 0 0 1 0 0 0 0 19
20 629 18.5 0 0 0 0 0 0 0 1 0 0 0 20
21 628 18.7 0 0 0 0 0 0 0 0 1 0 0 21
22 612 18.5 0 0 0 0 0 0 0 0 0 1 0 22
23 595 19.3 0 0 0 0 0 0 0 0 0 0 1 23
24 597 18.9 0 0 0 0 0 0 0 0 0 0 0 24
25 593 21.4 1 0 0 0 0 0 0 0 0 0 0 25
26 590 22.5 0 1 0 0 0 0 0 0 0 0 0 26
27 580 25.0 0 0 1 0 0 0 0 0 0 0 0 27
28 574 22.9 0 0 0 1 0 0 0 0 0 0 0 28
29 573 22.9 0 0 0 0 1 0 0 0 0 0 0 29
30 573 21.3 0 0 0 0 0 1 0 0 0 0 0 30
31 620 22.3 0 0 0 0 0 0 1 0 0 0 0 31
32 626 20.9 0 0 0 0 0 0 0 1 0 0 0 32
33 620 19.9 0 0 0 0 0 0 0 0 1 0 0 33
34 588 20.2 0 0 0 0 0 0 0 0 0 1 0 34
35 566 19.8 0 0 0 0 0 0 0 0 0 0 1 35
36 557 17.7 0 0 0 0 0 0 0 0 0 0 0 36
37 561 18.1 1 0 0 0 0 0 0 0 0 0 0 37
38 549 17.6 0 1 0 0 0 0 0 0 0 0 0 38
39 532 18.2 0 0 1 0 0 0 0 0 0 0 0 39
40 526 16.0 0 0 0 1 0 0 0 0 0 0 0 40
41 511 16.3 0 0 0 0 1 0 0 0 0 0 0 41
42 499 17.3 0 0 0 0 0 1 0 0 0 0 0 42
43 555 19.0 0 0 0 0 0 0 1 0 0 0 0 43
44 565 18.6 0 0 0 0 0 0 0 1 0 0 0 44
45 542 18.0 0 0 0 0 0 0 0 0 1 0 0 45
46 527 17.9 0 0 0 0 0 0 0 0 0 1 0 46
47 510 17.8 0 0 0 0 0 0 0 0 0 0 1 47
48 514 18.5 0 0 0 0 0 0 0 0 0 0 0 48
49 517 17.4 1 0 0 0 0 0 0 0 0 0 0 49
50 508 19.0 0 1 0 0 0 0 0 0 0 0 0 50
51 493 17.4 0 0 1 0 0 0 0 0 0 0 0 51
52 490 20.6 0 0 0 1 0 0 0 0 0 0 0 52
53 469 18.5 0 0 0 0 1 0 0 0 0 0 0 53
54 478 20.0 0 0 0 0 0 1 0 0 0 0 0 54
55 528 18.8 0 0 0 0 0 0 1 0 0 0 0 55
56 534 18.8 0 0 0 0 0 0 0 1 0 0 0 56
57 518 19.7 0 0 0 0 0 0 0 0 1 0 0 57
58 506 15.3 0 0 0 0 0 0 0 0 0 1 0 58
59 502 10.6 0 0 0 0 0 0 0 0 0 0 1 59
60 516 6.1 0 0 0 0 0 0 0 0 0 0 0 60
61 528 0.9 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
582.266 1.759 -3.134 -15.905 -25.592 -31.833
M5 M6 M7 M8 M9 M10
-40.006 -38.022 14.477 25.248 18.873 6.742
M11 t
-5.599 -1.546
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-40.0296 -19.9952 -0.7344 18.6025 41.5958
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 582.2656 24.0444 24.216 < 2e-16 ***
X 1.7586 1.0432 1.686 0.0985 .
M1 -3.1344 15.1269 -0.207 0.8367
M2 -15.9054 15.8838 -1.001 0.3218
M3 -25.5924 15.9157 -1.608 0.1145
M4 -31.8332 15.9406 -1.997 0.0516 .
M5 -40.0058 15.9126 -2.514 0.0154 *
M6 -38.0225 15.9634 -2.382 0.0213 *
M7 14.4774 15.9918 0.905 0.3699
M8 25.2476 16.0173 1.576 0.1217
M9 18.8729 15.9310 1.185 0.2421
M10 6.7424 15.7963 0.427 0.6714
M11 -5.5991 15.7257 -0.356 0.7234
t -1.5461 0.2009 -7.697 7.24e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 24.82 on 47 degrees of freedom
Multiple R-squared: 0.724, Adjusted R-squared: 0.6476
F-statistic: 9.483 on 13 and 47 DF, p-value: 3.276e-09
> 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,] 6.292501e-05 1.258500e-04 0.9999370750
[2,] 6.603280e-05 1.320656e-04 0.9999339672
[3,] 8.016894e-06 1.603379e-05 0.9999919831
[4,] 2.494548e-04 4.989096e-04 0.9997505452
[5,] 7.581975e-04 1.516395e-03 0.9992418025
[6,] 1.901897e-02 3.803793e-02 0.9809810348
[7,] 4.417704e-02 8.835409e-02 0.9558229570
[8,] 4.487560e-02 8.975119e-02 0.9551244035
[9,] 5.688126e-02 1.137625e-01 0.9431187439
[10,] 5.113680e-02 1.022736e-01 0.9488632011
[11,] 5.251221e-02 1.050244e-01 0.9474877929
[12,] 3.414729e-02 6.829458e-02 0.9658527109
[13,] 3.631869e-02 7.263739e-02 0.9636813062
[14,] 3.883933e-02 7.767865e-02 0.9611606743
[15,] 4.724820e-02 9.449640e-02 0.9527518007
[16,] 7.251400e-02 1.450280e-01 0.9274860024
[17,] 3.452382e-01 6.904764e-01 0.6547618206
[18,] 8.885139e-01 2.229723e-01 0.1114861436
[19,] 9.939282e-01 1.214360e-02 0.0060718007
[20,] 9.972875e-01 5.425038e-03 0.0027125192
[21,] 9.990750e-01 1.849997e-03 0.0009249984
[22,] 9.991031e-01 1.793783e-03 0.0008968915
[23,] 9.997149e-01 5.702838e-04 0.0002851419
[24,] 9.989824e-01 2.035205e-03 0.0010176026
[25,] 9.995858e-01 8.284414e-04 0.0004142207
[26,] 9.998315e-01 3.370642e-04 0.0001685321
[27,] 9.990866e-01 1.826754e-03 0.0009133770
[28,] 9.989871e-01 2.025769e-03 0.0010128844
> postscript(file="/var/www/html/rcomp/tmp/1xkpi1258757781.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/2i9g11258757781.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/32zkn1258757781.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/4ryyg1258757781.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/5hwzj1258757781.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 = 61
Frequency = 1
1 2 3 4 5 6
-40.0296073 -30.2297138 -29.0413772 -34.8923716 -32.8771816 -26.9627152
7 8 9 10 11 12
-27.0923716 -23.8888451 -10.4509610 2.1049161 -1.1833297 -7.0501467
13 14 15 16 17 18
-6.7214404 9.9921054 19.7423764 15.2982535 17.9617255 17.9417437
19 20 21 22 23 24
20.5706766 19.8741124 26.4430998 24.4714001 19.9521417 18.6025034
25 26 27 28 29 30
14.8864542 24.2690782 21.1056969 26.5856243 35.3042500 37.6807416
31 32 33 34 35 36
31.9683547 31.2062388 34.8855335 16.0345391 8.6255880 -0.7344483
37 38 39 40 41 42
7.2425402 10.4389072 3.6168457 9.2726320 3.4636809 -10.7321599
43 44 45 46 47 48
-8.6755593 -7.1962646 -21.2204057 -22.3679643 -25.3044923 -26.5885789
49 50 51 52 53 54
-16.9737063 -14.4703770 -15.4235418 -16.2641382 -23.8524748 -17.9276103
55 56 57 58 59 60
-16.7711005 -19.9952415 -29.6572667 -20.2428910 -2.0899077 15.7706705
61
41.5957596
> postscript(file="/var/www/html/rcomp/tmp/6a6e71258757781.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -40.0296073 NA
1 -30.2297138 -40.0296073
2 -29.0413772 -30.2297138
3 -34.8923716 -29.0413772
4 -32.8771816 -34.8923716
5 -26.9627152 -32.8771816
6 -27.0923716 -26.9627152
7 -23.8888451 -27.0923716
8 -10.4509610 -23.8888451
9 2.1049161 -10.4509610
10 -1.1833297 2.1049161
11 -7.0501467 -1.1833297
12 -6.7214404 -7.0501467
13 9.9921054 -6.7214404
14 19.7423764 9.9921054
15 15.2982535 19.7423764
16 17.9617255 15.2982535
17 17.9417437 17.9617255
18 20.5706766 17.9417437
19 19.8741124 20.5706766
20 26.4430998 19.8741124
21 24.4714001 26.4430998
22 19.9521417 24.4714001
23 18.6025034 19.9521417
24 14.8864542 18.6025034
25 24.2690782 14.8864542
26 21.1056969 24.2690782
27 26.5856243 21.1056969
28 35.3042500 26.5856243
29 37.6807416 35.3042500
30 31.9683547 37.6807416
31 31.2062388 31.9683547
32 34.8855335 31.2062388
33 16.0345391 34.8855335
34 8.6255880 16.0345391
35 -0.7344483 8.6255880
36 7.2425402 -0.7344483
37 10.4389072 7.2425402
38 3.6168457 10.4389072
39 9.2726320 3.6168457
40 3.4636809 9.2726320
41 -10.7321599 3.4636809
42 -8.6755593 -10.7321599
43 -7.1962646 -8.6755593
44 -21.2204057 -7.1962646
45 -22.3679643 -21.2204057
46 -25.3044923 -22.3679643
47 -26.5885789 -25.3044923
48 -16.9737063 -26.5885789
49 -14.4703770 -16.9737063
50 -15.4235418 -14.4703770
51 -16.2641382 -15.4235418
52 -23.8524748 -16.2641382
53 -17.9276103 -23.8524748
54 -16.7711005 -17.9276103
55 -19.9952415 -16.7711005
56 -29.6572667 -19.9952415
57 -20.2428910 -29.6572667
58 -2.0899077 -20.2428910
59 15.7706705 -2.0899077
60 41.5957596 15.7706705
61 NA 41.5957596
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -30.2297138 -40.0296073
[2,] -29.0413772 -30.2297138
[3,] -34.8923716 -29.0413772
[4,] -32.8771816 -34.8923716
[5,] -26.9627152 -32.8771816
[6,] -27.0923716 -26.9627152
[7,] -23.8888451 -27.0923716
[8,] -10.4509610 -23.8888451
[9,] 2.1049161 -10.4509610
[10,] -1.1833297 2.1049161
[11,] -7.0501467 -1.1833297
[12,] -6.7214404 -7.0501467
[13,] 9.9921054 -6.7214404
[14,] 19.7423764 9.9921054
[15,] 15.2982535 19.7423764
[16,] 17.9617255 15.2982535
[17,] 17.9417437 17.9617255
[18,] 20.5706766 17.9417437
[19,] 19.8741124 20.5706766
[20,] 26.4430998 19.8741124
[21,] 24.4714001 26.4430998
[22,] 19.9521417 24.4714001
[23,] 18.6025034 19.9521417
[24,] 14.8864542 18.6025034
[25,] 24.2690782 14.8864542
[26,] 21.1056969 24.2690782
[27,] 26.5856243 21.1056969
[28,] 35.3042500 26.5856243
[29,] 37.6807416 35.3042500
[30,] 31.9683547 37.6807416
[31,] 31.2062388 31.9683547
[32,] 34.8855335 31.2062388
[33,] 16.0345391 34.8855335
[34,] 8.6255880 16.0345391
[35,] -0.7344483 8.6255880
[36,] 7.2425402 -0.7344483
[37,] 10.4389072 7.2425402
[38,] 3.6168457 10.4389072
[39,] 9.2726320 3.6168457
[40,] 3.4636809 9.2726320
[41,] -10.7321599 3.4636809
[42,] -8.6755593 -10.7321599
[43,] -7.1962646 -8.6755593
[44,] -21.2204057 -7.1962646
[45,] -22.3679643 -21.2204057
[46,] -25.3044923 -22.3679643
[47,] -26.5885789 -25.3044923
[48,] -16.9737063 -26.5885789
[49,] -14.4703770 -16.9737063
[50,] -15.4235418 -14.4703770
[51,] -16.2641382 -15.4235418
[52,] -23.8524748 -16.2641382
[53,] -17.9276103 -23.8524748
[54,] -16.7711005 -17.9276103
[55,] -19.9952415 -16.7711005
[56,] -29.6572667 -19.9952415
[57,] -20.2428910 -29.6572667
[58,] -2.0899077 -20.2428910
[59,] 15.7706705 -2.0899077
[60,] 41.5957596 15.7706705
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -30.2297138 -40.0296073
2 -29.0413772 -30.2297138
3 -34.8923716 -29.0413772
4 -32.8771816 -34.8923716
5 -26.9627152 -32.8771816
6 -27.0923716 -26.9627152
7 -23.8888451 -27.0923716
8 -10.4509610 -23.8888451
9 2.1049161 -10.4509610
10 -1.1833297 2.1049161
11 -7.0501467 -1.1833297
12 -6.7214404 -7.0501467
13 9.9921054 -6.7214404
14 19.7423764 9.9921054
15 15.2982535 19.7423764
16 17.9617255 15.2982535
17 17.9417437 17.9617255
18 20.5706766 17.9417437
19 19.8741124 20.5706766
20 26.4430998 19.8741124
21 24.4714001 26.4430998
22 19.9521417 24.4714001
23 18.6025034 19.9521417
24 14.8864542 18.6025034
25 24.2690782 14.8864542
26 21.1056969 24.2690782
27 26.5856243 21.1056969
28 35.3042500 26.5856243
29 37.6807416 35.3042500
30 31.9683547 37.6807416
31 31.2062388 31.9683547
32 34.8855335 31.2062388
33 16.0345391 34.8855335
34 8.6255880 16.0345391
35 -0.7344483 8.6255880
36 7.2425402 -0.7344483
37 10.4389072 7.2425402
38 3.6168457 10.4389072
39 9.2726320 3.6168457
40 3.4636809 9.2726320
41 -10.7321599 3.4636809
42 -8.6755593 -10.7321599
43 -7.1962646 -8.6755593
44 -21.2204057 -7.1962646
45 -22.3679643 -21.2204057
46 -25.3044923 -22.3679643
47 -26.5885789 -25.3044923
48 -16.9737063 -26.5885789
49 -14.4703770 -16.9737063
50 -15.4235418 -14.4703770
51 -16.2641382 -15.4235418
52 -23.8524748 -16.2641382
53 -17.9276103 -23.8524748
54 -16.7711005 -17.9276103
55 -19.9952415 -16.7711005
56 -29.6572667 -19.9952415
57 -20.2428910 -29.6572667
58 -2.0899077 -20.2428910
59 15.7706705 -2.0899077
60 41.5957596 15.7706705
> 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/7re7u1258757781.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/897881258757781.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/98fao1258757781.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/10yyeo1258757781.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/11p8yl1258757781.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/12sbxz1258757781.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/13skvb1258757781.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/142n5e1258757781.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/152xmz1258757781.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/16h5da1258757781.tab")
+ }
>
> system("convert tmp/1xkpi1258757781.ps tmp/1xkpi1258757781.png")
> system("convert tmp/2i9g11258757781.ps tmp/2i9g11258757781.png")
> system("convert tmp/32zkn1258757781.ps tmp/32zkn1258757781.png")
> system("convert tmp/4ryyg1258757781.ps tmp/4ryyg1258757781.png")
> system("convert tmp/5hwzj1258757781.ps tmp/5hwzj1258757781.png")
> system("convert tmp/6a6e71258757781.ps tmp/6a6e71258757781.png")
> system("convert tmp/7re7u1258757781.ps tmp/7re7u1258757781.png")
> system("convert tmp/897881258757781.ps tmp/897881258757781.png")
> system("convert tmp/98fao1258757781.ps tmp/98fao1258757781.png")
> system("convert tmp/10yyeo1258757781.ps tmp/10yyeo1258757781.png")
>
>
> proc.time()
user system elapsed
2.393 1.556 3.265