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(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,561,106,549,105.3,532,118.8,526,106.1,511,109.3,499,117.2,555,92.5,565,104.2,542,112.5,527,122.4,510,113.3,514,100,517,110.7,508,112.8,493,109.8,490,117.3,469,109.1,478,115.9,528,96,534,99.8,518,116.8,506,115.7,502,99.4,516,94.3),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 519 97.4 1 0 0 0 0 0 0 0 0 0 0
2 517 97.0 0 1 0 0 0 0 0 0 0 0 0
3 510 105.4 0 0 1 0 0 0 0 0 0 0 0
4 509 102.7 0 0 0 1 0 0 0 0 0 0 0
5 501 98.1 0 0 0 0 1 0 0 0 0 0 0
6 507 104.5 0 0 0 0 0 1 0 0 0 0 0
7 569 87.4 0 0 0 0 0 0 1 0 0 0 0
8 580 89.9 0 0 0 0 0 0 0 1 0 0 0
9 578 109.8 0 0 0 0 0 0 0 0 1 0 0
10 565 111.7 0 0 0 0 0 0 0 0 0 1 0
11 547 98.6 0 0 0 0 0 0 0 0 0 0 1
12 555 96.9 0 0 0 0 0 0 0 0 0 0 0
13 562 95.1 1 0 0 0 0 0 0 0 0 0 0
14 561 97.0 0 1 0 0 0 0 0 0 0 0 0
15 555 112.7 0 0 1 0 0 0 0 0 0 0 0
16 544 102.9 0 0 0 1 0 0 0 0 0 0 0
17 537 97.4 0 0 0 0 1 0 0 0 0 0 0
18 543 111.4 0 0 0 0 0 1 0 0 0 0 0
19 594 87.4 0 0 0 0 0 0 1 0 0 0 0
20 611 96.8 0 0 0 0 0 0 0 1 0 0 0
21 613 114.1 0 0 0 0 0 0 0 0 1 0 0
22 611 110.3 0 0 0 0 0 0 0 0 0 1 0
23 594 103.9 0 0 0 0 0 0 0 0 0 0 1
24 595 101.6 0 0 0 0 0 0 0 0 0 0 0
25 591 94.6 1 0 0 0 0 0 0 0 0 0 0
26 589 95.9 0 1 0 0 0 0 0 0 0 0 0
27 584 104.7 0 0 1 0 0 0 0 0 0 0 0
28 573 102.8 0 0 0 1 0 0 0 0 0 0 0
29 567 98.1 0 0 0 0 1 0 0 0 0 0 0
30 569 113.9 0 0 0 0 0 1 0 0 0 0 0
31 621 80.9 0 0 0 0 0 0 1 0 0 0 0
32 629 95.7 0 0 0 0 0 0 0 1 0 0 0
33 628 113.2 0 0 0 0 0 0 0 0 1 0 0
34 612 105.9 0 0 0 0 0 0 0 0 0 1 0
35 595 108.8 0 0 0 0 0 0 0 0 0 0 1
36 597 102.3 0 0 0 0 0 0 0 0 0 0 0
37 593 99.0 1 0 0 0 0 0 0 0 0 0 0
38 590 100.7 0 1 0 0 0 0 0 0 0 0 0
39 580 115.5 0 0 1 0 0 0 0 0 0 0 0
40 574 100.7 0 0 0 1 0 0 0 0 0 0 0
41 573 109.9 0 0 0 0 1 0 0 0 0 0 0
42 573 114.6 0 0 0 0 0 1 0 0 0 0 0
43 620 85.4 0 0 0 0 0 0 1 0 0 0 0
44 626 100.5 0 0 0 0 0 0 0 1 0 0 0
45 620 114.8 0 0 0 0 0 0 0 0 1 0 0
46 588 116.5 0 0 0 0 0 0 0 0 0 1 0
47 566 112.9 0 0 0 0 0 0 0 0 0 0 1
48 557 102.0 0 0 0 0 0 0 0 0 0 0 0
49 561 106.0 1 0 0 0 0 0 0 0 0 0 0
50 549 105.3 0 1 0 0 0 0 0 0 0 0 0
51 532 118.8 0 0 1 0 0 0 0 0 0 0 0
52 526 106.1 0 0 0 1 0 0 0 0 0 0 0
53 511 109.3 0 0 0 0 1 0 0 0 0 0 0
54 499 117.2 0 0 0 0 0 1 0 0 0 0 0
55 555 92.5 0 0 0 0 0 0 1 0 0 0 0
56 565 104.2 0 0 0 0 0 0 0 1 0 0 0
57 542 112.5 0 0 0 0 0 0 0 0 1 0 0
58 527 122.4 0 0 0 0 0 0 0 0 0 1 0
59 510 113.3 0 0 0 0 0 0 0 0 0 0 1
60 514 100.0 0 0 0 0 0 0 0 0 0 0 0
61 517 110.7 1 0 0 0 0 0 0 0 0 0 0
62 508 112.8 0 1 0 0 0 0 0 0 0 0 0
63 493 109.8 0 0 1 0 0 0 0 0 0 0 0
64 490 117.3 0 0 0 1 0 0 0 0 0 0 0
65 469 109.1 0 0 0 0 1 0 0 0 0 0 0
66 478 115.9 0 0 0 0 0 1 0 0 0 0 0
67 528 96.0 0 0 0 0 0 0 1 0 0 0 0
68 534 99.8 0 0 0 0 0 0 0 1 0 0 0
69 518 116.8 0 0 0 0 0 0 0 0 1 0 0
70 506 115.7 0 0 0 0 0 0 0 0 0 1 0
71 502 99.4 0 0 0 0 0 0 0 0 0 0 1
72 516 94.3 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
760.9105 -2.0624 3.4593 0.6540 10.6593 -7.4985
M5 M6 M7 M8 M9 M10
-20.8087 0.1362 2.2979 31.6606 56.4081 41.8549
M11
10.3473
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-64.255 -27.267 1.149 29.666 59.557
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 760.9105 88.5940 8.589 5.57e-12 ***
X -2.0624 0.8771 -2.351 0.0221 *
M1 3.4593 21.4323 0.161 0.8723
M2 0.6540 21.4831 0.030 0.9758
M3 10.6593 23.7228 0.449 0.6548
M4 -7.4985 22.0325 -0.340 0.7348
M5 -20.8087 21.7208 -0.958 0.3420
M6 0.1362 24.4294 0.006 0.9956
M7 2.2979 23.5801 0.097 0.9227
M8 31.6606 21.4679 1.475 0.1456
M9 56.4081 24.6942 2.284 0.0260 *
M10 41.8549 24.7894 1.688 0.0966 .
M11 10.3473 22.1924 0.466 0.6428
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 37.09 on 59 degrees of freedom
Multiple R-squared: 0.3167, Adjusted R-squared: 0.1777
F-statistic: 2.279 on 12 and 59 DF, p-value: 0.01866
> 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.43583448 0.87166896 0.56416552
[2,] 0.35928244 0.71856488 0.64071756
[3,] 0.22572776 0.45145552 0.77427224
[4,] 0.15255992 0.30511984 0.84744008
[5,] 0.08490968 0.16981935 0.91509032
[6,] 0.05008883 0.10017767 0.94991117
[7,] 0.06345623 0.12691246 0.93654377
[8,] 0.04165809 0.08331619 0.95834191
[9,] 0.02787147 0.05574293 0.97212853
[10,] 0.05128672 0.10257344 0.94871328
[11,] 0.06741029 0.13482058 0.93258971
[12,] 0.10577345 0.21154690 0.89422655
[13,] 0.10247659 0.20495317 0.89752341
[14,] 0.09617128 0.19234257 0.90382872
[15,] 0.07726899 0.15453797 0.92273101
[16,] 0.08509340 0.17018681 0.91490660
[17,] 0.06753051 0.13506102 0.93246949
[18,] 0.06408603 0.12817206 0.93591397
[19,] 0.05494548 0.10989096 0.94505452
[20,] 0.05095160 0.10190320 0.94904840
[21,] 0.05480638 0.10961275 0.94519362
[22,] 0.04308824 0.08617647 0.95691176
[23,] 0.03506745 0.07013491 0.96493255
[24,] 0.03446393 0.06892787 0.96553607
[25,] 0.03447338 0.06894677 0.96552662
[26,] 0.04554760 0.09109521 0.95445240
[27,] 0.06773523 0.13547047 0.93226477
[28,] 0.09645262 0.19290524 0.90354738
[29,] 0.14390129 0.28780258 0.85609871
[30,] 0.31641673 0.63283347 0.68358327
[31,] 0.52115121 0.95769758 0.47884879
[32,] 0.68105007 0.63789985 0.31894993
[33,] 0.74544866 0.50910269 0.25455134
[34,] 0.77120685 0.45758629 0.22879315
[35,] 0.80177868 0.39644264 0.19822132
[36,] 0.83613117 0.32773767 0.16386883
[37,] 0.85536608 0.28926783 0.14463392
[38,] 0.91325143 0.17349715 0.08674857
[39,] 0.88884257 0.22231486 0.11115743
[40,] 0.89342334 0.21315332 0.10657666
[41,] 0.89505161 0.20989678 0.10494839
> postscript(file="/var/www/html/rcomp/tmp/1lel91258622668.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/2w6pn1258622668.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/32lc51258622668.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/4r2x01258622668.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/5p0n01258622668.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
-44.4913793 -44.5110413 -44.1921695 -32.6028704 -36.7796882 -38.5252529
7 8 9 10 11 12
-13.9540854 -27.1607164 -12.8663168 -7.3945995 -20.9045008 -6.0632964
13 14 15 16 17 18
-6.2349137 -0.5110413 15.8633964 2.8096109 -2.2233726 11.7053505
19 20 21 22 23 24
11.0459146 18.0698869 31.0020302 35.7180317 37.0262525 43.6300131
25 26 27 28 29 30
21.7338831 25.2203118 28.3641461 31.6033702 29.2203118 42.8613662
31 32 33 34 35 36
24.6402737 33.8012400 44.1458646 27.6434440 48.1320433 47.0736975
37 38 39 40 41 42
32.8084708 36.1198620 46.6381340 28.2723170 59.5567059 48.3050506
43 44 45 46 47 48
32.9211020 40.7007902 39.4457146 25.5049506 27.5879091 6.4549756
49 50 51 52 53 54
15.2453148 4.6069309 5.4440747 -8.5906890 -3.6807378 -20.3326931
55 56 57 58 59 60
-17.4358134 -12.6683066 -43.2978198 -23.3268523 -27.5871284 -40.6698370
61 62 63 64 65 66
-19.0613757 -20.9250220 -52.1175818 -21.4917387 -46.0932191 -44.0138213
67 68 69 70 71 72
-37.2173914 -52.7428942 -58.4294728 -58.1449744 -64.2545757 -50.4255528
> postscript(file="/var/www/html/rcomp/tmp/69qj31258622668.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 -44.4913793 NA
1 -44.5110413 -44.4913793
2 -44.1921695 -44.5110413
3 -32.6028704 -44.1921695
4 -36.7796882 -32.6028704
5 -38.5252529 -36.7796882
6 -13.9540854 -38.5252529
7 -27.1607164 -13.9540854
8 -12.8663168 -27.1607164
9 -7.3945995 -12.8663168
10 -20.9045008 -7.3945995
11 -6.0632964 -20.9045008
12 -6.2349137 -6.0632964
13 -0.5110413 -6.2349137
14 15.8633964 -0.5110413
15 2.8096109 15.8633964
16 -2.2233726 2.8096109
17 11.7053505 -2.2233726
18 11.0459146 11.7053505
19 18.0698869 11.0459146
20 31.0020302 18.0698869
21 35.7180317 31.0020302
22 37.0262525 35.7180317
23 43.6300131 37.0262525
24 21.7338831 43.6300131
25 25.2203118 21.7338831
26 28.3641461 25.2203118
27 31.6033702 28.3641461
28 29.2203118 31.6033702
29 42.8613662 29.2203118
30 24.6402737 42.8613662
31 33.8012400 24.6402737
32 44.1458646 33.8012400
33 27.6434440 44.1458646
34 48.1320433 27.6434440
35 47.0736975 48.1320433
36 32.8084708 47.0736975
37 36.1198620 32.8084708
38 46.6381340 36.1198620
39 28.2723170 46.6381340
40 59.5567059 28.2723170
41 48.3050506 59.5567059
42 32.9211020 48.3050506
43 40.7007902 32.9211020
44 39.4457146 40.7007902
45 25.5049506 39.4457146
46 27.5879091 25.5049506
47 6.4549756 27.5879091
48 15.2453148 6.4549756
49 4.6069309 15.2453148
50 5.4440747 4.6069309
51 -8.5906890 5.4440747
52 -3.6807378 -8.5906890
53 -20.3326931 -3.6807378
54 -17.4358134 -20.3326931
55 -12.6683066 -17.4358134
56 -43.2978198 -12.6683066
57 -23.3268523 -43.2978198
58 -27.5871284 -23.3268523
59 -40.6698370 -27.5871284
60 -19.0613757 -40.6698370
61 -20.9250220 -19.0613757
62 -52.1175818 -20.9250220
63 -21.4917387 -52.1175818
64 -46.0932191 -21.4917387
65 -44.0138213 -46.0932191
66 -37.2173914 -44.0138213
67 -52.7428942 -37.2173914
68 -58.4294728 -52.7428942
69 -58.1449744 -58.4294728
70 -64.2545757 -58.1449744
71 -50.4255528 -64.2545757
72 NA -50.4255528
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -44.5110413 -44.4913793
[2,] -44.1921695 -44.5110413
[3,] -32.6028704 -44.1921695
[4,] -36.7796882 -32.6028704
[5,] -38.5252529 -36.7796882
[6,] -13.9540854 -38.5252529
[7,] -27.1607164 -13.9540854
[8,] -12.8663168 -27.1607164
[9,] -7.3945995 -12.8663168
[10,] -20.9045008 -7.3945995
[11,] -6.0632964 -20.9045008
[12,] -6.2349137 -6.0632964
[13,] -0.5110413 -6.2349137
[14,] 15.8633964 -0.5110413
[15,] 2.8096109 15.8633964
[16,] -2.2233726 2.8096109
[17,] 11.7053505 -2.2233726
[18,] 11.0459146 11.7053505
[19,] 18.0698869 11.0459146
[20,] 31.0020302 18.0698869
[21,] 35.7180317 31.0020302
[22,] 37.0262525 35.7180317
[23,] 43.6300131 37.0262525
[24,] 21.7338831 43.6300131
[25,] 25.2203118 21.7338831
[26,] 28.3641461 25.2203118
[27,] 31.6033702 28.3641461
[28,] 29.2203118 31.6033702
[29,] 42.8613662 29.2203118
[30,] 24.6402737 42.8613662
[31,] 33.8012400 24.6402737
[32,] 44.1458646 33.8012400
[33,] 27.6434440 44.1458646
[34,] 48.1320433 27.6434440
[35,] 47.0736975 48.1320433
[36,] 32.8084708 47.0736975
[37,] 36.1198620 32.8084708
[38,] 46.6381340 36.1198620
[39,] 28.2723170 46.6381340
[40,] 59.5567059 28.2723170
[41,] 48.3050506 59.5567059
[42,] 32.9211020 48.3050506
[43,] 40.7007902 32.9211020
[44,] 39.4457146 40.7007902
[45,] 25.5049506 39.4457146
[46,] 27.5879091 25.5049506
[47,] 6.4549756 27.5879091
[48,] 15.2453148 6.4549756
[49,] 4.6069309 15.2453148
[50,] 5.4440747 4.6069309
[51,] -8.5906890 5.4440747
[52,] -3.6807378 -8.5906890
[53,] -20.3326931 -3.6807378
[54,] -17.4358134 -20.3326931
[55,] -12.6683066 -17.4358134
[56,] -43.2978198 -12.6683066
[57,] -23.3268523 -43.2978198
[58,] -27.5871284 -23.3268523
[59,] -40.6698370 -27.5871284
[60,] -19.0613757 -40.6698370
[61,] -20.9250220 -19.0613757
[62,] -52.1175818 -20.9250220
[63,] -21.4917387 -52.1175818
[64,] -46.0932191 -21.4917387
[65,] -44.0138213 -46.0932191
[66,] -37.2173914 -44.0138213
[67,] -52.7428942 -37.2173914
[68,] -58.4294728 -52.7428942
[69,] -58.1449744 -58.4294728
[70,] -64.2545757 -58.1449744
[71,] -50.4255528 -64.2545757
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -44.5110413 -44.4913793
2 -44.1921695 -44.5110413
3 -32.6028704 -44.1921695
4 -36.7796882 -32.6028704
5 -38.5252529 -36.7796882
6 -13.9540854 -38.5252529
7 -27.1607164 -13.9540854
8 -12.8663168 -27.1607164
9 -7.3945995 -12.8663168
10 -20.9045008 -7.3945995
11 -6.0632964 -20.9045008
12 -6.2349137 -6.0632964
13 -0.5110413 -6.2349137
14 15.8633964 -0.5110413
15 2.8096109 15.8633964
16 -2.2233726 2.8096109
17 11.7053505 -2.2233726
18 11.0459146 11.7053505
19 18.0698869 11.0459146
20 31.0020302 18.0698869
21 35.7180317 31.0020302
22 37.0262525 35.7180317
23 43.6300131 37.0262525
24 21.7338831 43.6300131
25 25.2203118 21.7338831
26 28.3641461 25.2203118
27 31.6033702 28.3641461
28 29.2203118 31.6033702
29 42.8613662 29.2203118
30 24.6402737 42.8613662
31 33.8012400 24.6402737
32 44.1458646 33.8012400
33 27.6434440 44.1458646
34 48.1320433 27.6434440
35 47.0736975 48.1320433
36 32.8084708 47.0736975
37 36.1198620 32.8084708
38 46.6381340 36.1198620
39 28.2723170 46.6381340
40 59.5567059 28.2723170
41 48.3050506 59.5567059
42 32.9211020 48.3050506
43 40.7007902 32.9211020
44 39.4457146 40.7007902
45 25.5049506 39.4457146
46 27.5879091 25.5049506
47 6.4549756 27.5879091
48 15.2453148 6.4549756
49 4.6069309 15.2453148
50 5.4440747 4.6069309
51 -8.5906890 5.4440747
52 -3.6807378 -8.5906890
53 -20.3326931 -3.6807378
54 -17.4358134 -20.3326931
55 -12.6683066 -17.4358134
56 -43.2978198 -12.6683066
57 -23.3268523 -43.2978198
58 -27.5871284 -23.3268523
59 -40.6698370 -27.5871284
60 -19.0613757 -40.6698370
61 -20.9250220 -19.0613757
62 -52.1175818 -20.9250220
63 -21.4917387 -52.1175818
64 -46.0932191 -21.4917387
65 -44.0138213 -46.0932191
66 -37.2173914 -44.0138213
67 -52.7428942 -37.2173914
68 -58.4294728 -52.7428942
69 -58.1449744 -58.4294728
70 -64.2545757 -58.1449744
71 -50.4255528 -64.2545757
> 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/7panl1258622668.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/8lqn01258622668.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/9ajhd1258622668.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/10bclj1258622668.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/113unq1258622668.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/12gaqh1258622668.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/13j6yi1258622668.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/14t3ny1258622668.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/15gpal1258622668.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/16203t1258622668.tab")
+ }
>
> system("convert tmp/1lel91258622668.ps tmp/1lel91258622668.png")
> system("convert tmp/2w6pn1258622668.ps tmp/2w6pn1258622668.png")
> system("convert tmp/32lc51258622668.ps tmp/32lc51258622668.png")
> system("convert tmp/4r2x01258622668.ps tmp/4r2x01258622668.png")
> system("convert tmp/5p0n01258622668.ps tmp/5p0n01258622668.png")
> system("convert tmp/69qj31258622668.ps tmp/69qj31258622668.png")
> system("convert tmp/7panl1258622668.ps tmp/7panl1258622668.png")
> system("convert tmp/8lqn01258622668.ps tmp/8lqn01258622668.png")
> system("convert tmp/9ajhd1258622668.ps tmp/9ajhd1258622668.png")
> system("convert tmp/10bclj1258622668.ps tmp/10bclj1258622668.png")
>
>
> proc.time()
user system elapsed
2.534 1.558 2.972