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(108.5,98.71,112.3,98.54,116.6,98.2,115.5,96.92,120.1,99.06,132.9,99.65,128.1,99.82,129.3,99.99,132.5,100.33,131,99.31,124.9,101.1,120.8,101.1,122,100.93,122.1,100.85,127.4,100.93,135.2,99.6,137.3,101.88,135,101.81,136,102.38,138.4,102.74,134.7,102.82,138.4,101.72,133.9,103.47,133.6,102.98,141.2,102.68,151.8,102.9,155.4,103.03,156.6,101.29,161.6,103.69,160.7,103.68,156,104.2,159.5,104.08,168.7,104.16,169.9,103.05,169.9,104.66,185.9,104.46,190.8,104.95,195.8,105.85,211.9,106.23,227.1,104.86,251.3,107.44,256.7,108.23,251.9,108.45,251.2,109.39,270.3,110.15,267.2,109.13,243,110.28,229.9,110.17,187.2,109.99,178.2,109.26,175.2,109.11,192.4,107.06,187,109.53,184,108.92,194.1,109.24,212.7,109.12,217.5,109,200.5,107.23,205.9,109.49,196.5,109.04,206.3,109.02),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 = '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 108.5 98.71 1 0 0 0 0 0 0 0 0 0 0
2 112.3 98.54 0 1 0 0 0 0 0 0 0 0 0
3 116.6 98.20 0 0 1 0 0 0 0 0 0 0 0
4 115.5 96.92 0 0 0 1 0 0 0 0 0 0 0
5 120.1 99.06 0 0 0 0 1 0 0 0 0 0 0
6 132.9 99.65 0 0 0 0 0 1 0 0 0 0 0
7 128.1 99.82 0 0 0 0 0 0 1 0 0 0 0
8 129.3 99.99 0 0 0 0 0 0 0 1 0 0 0
9 132.5 100.33 0 0 0 0 0 0 0 0 1 0 0
10 131.0 99.31 0 0 0 0 0 0 0 0 0 1 0
11 124.9 101.10 0 0 0 0 0 0 0 0 0 0 1
12 120.8 101.10 0 0 0 0 0 0 0 0 0 0 0
13 122.0 100.93 1 0 0 0 0 0 0 0 0 0 0
14 122.1 100.85 0 1 0 0 0 0 0 0 0 0 0
15 127.4 100.93 0 0 1 0 0 0 0 0 0 0 0
16 135.2 99.60 0 0 0 1 0 0 0 0 0 0 0
17 137.3 101.88 0 0 0 0 1 0 0 0 0 0 0
18 135.0 101.81 0 0 0 0 0 1 0 0 0 0 0
19 136.0 102.38 0 0 0 0 0 0 1 0 0 0 0
20 138.4 102.74 0 0 0 0 0 0 0 1 0 0 0
21 134.7 102.82 0 0 0 0 0 0 0 0 1 0 0
22 138.4 101.72 0 0 0 0 0 0 0 0 0 1 0
23 133.9 103.47 0 0 0 0 0 0 0 0 0 0 1
24 133.6 102.98 0 0 0 0 0 0 0 0 0 0 0
25 141.2 102.68 1 0 0 0 0 0 0 0 0 0 0
26 151.8 102.90 0 1 0 0 0 0 0 0 0 0 0
27 155.4 103.03 0 0 1 0 0 0 0 0 0 0 0
28 156.6 101.29 0 0 0 1 0 0 0 0 0 0 0
29 161.6 103.69 0 0 0 0 1 0 0 0 0 0 0
30 160.7 103.68 0 0 0 0 0 1 0 0 0 0 0
31 156.0 104.20 0 0 0 0 0 0 1 0 0 0 0
32 159.5 104.08 0 0 0 0 0 0 0 1 0 0 0
33 168.7 104.16 0 0 0 0 0 0 0 0 1 0 0
34 169.9 103.05 0 0 0 0 0 0 0 0 0 1 0
35 169.9 104.66 0 0 0 0 0 0 0 0 0 0 1
36 185.9 104.46 0 0 0 0 0 0 0 0 0 0 0
37 190.8 104.95 1 0 0 0 0 0 0 0 0 0 0
38 195.8 105.85 0 1 0 0 0 0 0 0 0 0 0
39 211.9 106.23 0 0 1 0 0 0 0 0 0 0 0
40 227.1 104.86 0 0 0 1 0 0 0 0 0 0 0
41 251.3 107.44 0 0 0 0 1 0 0 0 0 0 0
42 256.7 108.23 0 0 0 0 0 1 0 0 0 0 0
43 251.9 108.45 0 0 0 0 0 0 1 0 0 0 0
44 251.2 109.39 0 0 0 0 0 0 0 1 0 0 0
45 270.3 110.15 0 0 0 0 0 0 0 0 1 0 0
46 267.2 109.13 0 0 0 0 0 0 0 0 0 1 0
47 243.0 110.28 0 0 0 0 0 0 0 0 0 0 1
48 229.9 110.17 0 0 0 0 0 0 0 0 0 0 0
49 187.2 109.99 1 0 0 0 0 0 0 0 0 0 0
50 178.2 109.26 0 1 0 0 0 0 0 0 0 0 0
51 175.2 109.11 0 0 1 0 0 0 0 0 0 0 0
52 192.4 107.06 0 0 0 1 0 0 0 0 0 0 0
53 187.0 109.53 0 0 0 0 1 0 0 0 0 0 0
54 184.0 108.92 0 0 0 0 0 1 0 0 0 0 0
55 194.1 109.24 0 0 0 0 0 0 1 0 0 0 0
56 212.7 109.12 0 0 0 0 0 0 0 1 0 0 0
57 217.5 109.00 0 0 0 0 0 0 0 0 1 0 0
58 200.5 107.23 0 0 0 0 0 0 0 0 0 1 0
59 205.9 109.49 0 0 0 0 0 0 0 0 0 0 1
60 196.5 109.04 0 0 0 0 0 0 0 0 0 0 0
61 206.3 109.02 1 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
-918.2764 10.3422 -1.9063 0.1083 5.1615 29.2932
M5 M6 M7 M8 M9 M10
10.8409 11.8137 7.4505 9.9063 14.0683 23.1803
M11
-0.4055
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-40.120 -11.474 -1.976 8.765 47.572
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -918.2764 82.3590 -11.150 6.38e-15 ***
X 10.3422 0.7745 13.353 < 2e-16 ***
M1 -1.9063 13.5864 -0.140 0.8890
M2 0.1083 14.2494 0.008 0.9940
M3 5.1615 14.2477 0.362 0.7187
M4 29.2932 14.4314 2.030 0.0479 *
M5 10.8409 14.1909 0.764 0.4486
M6 11.8137 14.1842 0.833 0.4090
M7 7.4505 14.1703 0.526 0.6015
M8 9.9063 14.1639 0.699 0.4877
M9 14.0683 14.1603 0.993 0.3254
M10 23.1803 14.2041 1.632 0.1092
M11 -0.4055 14.1602 -0.029 0.9773
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 22.39 on 48 degrees of freedom
Multiple R-squared: 0.7979, Adjusted R-squared: 0.7474
F-statistic: 15.79 on 12 and 48 DF, p-value: 8.693e-13
> 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.856416e-03 3.712831e-03 0.9981436
[2,] 1.803185e-04 3.606369e-04 0.9998197
[3,] 1.507654e-04 3.015308e-04 0.9998492
[4,] 2.644648e-05 5.289297e-05 0.9999736
[5,] 4.068404e-06 8.136808e-06 0.9999959
[6,] 2.023777e-06 4.047554e-06 0.9999980
[7,] 3.087053e-07 6.174107e-07 0.9999997
[8,] 4.168373e-08 8.336745e-08 1.0000000
[9,] 8.974612e-09 1.794922e-08 1.0000000
[10,] 3.358772e-08 6.717544e-08 1.0000000
[11,] 2.090945e-07 4.181890e-07 0.9999998
[12,] 1.534779e-07 3.069557e-07 0.9999998
[13,] 8.568291e-08 1.713658e-07 0.9999999
[14,] 4.117819e-08 8.235638e-08 1.0000000
[15,] 1.114969e-08 2.229938e-08 1.0000000
[16,] 2.780395e-09 5.560790e-09 1.0000000
[17,] 9.907026e-10 1.981405e-09 1.0000000
[18,] 2.020924e-09 4.041848e-09 1.0000000
[19,] 3.598667e-09 7.197334e-09 1.0000000
[20,] 1.904633e-08 3.809265e-08 1.0000000
[21,] 9.895114e-07 1.979023e-06 0.9999990
[22,] 2.881598e-06 5.763197e-06 0.9999971
[23,] 1.578217e-06 3.156433e-06 0.9999984
[24,] 1.385504e-06 2.771008e-06 0.9999986
[25,] 2.791542e-06 5.583085e-06 0.9999972
[26,] 1.369762e-04 2.739524e-04 0.9998630
[27,] 4.980372e-03 9.960743e-03 0.9950196
[28,] 1.623888e-01 3.247777e-01 0.8376112
[29,] 1.835780e-01 3.671561e-01 0.8164220
[30,] 1.454655e-01 2.909310e-01 0.8545345
> postscript(file="/var/www/html/rcomp/tmp/111xp1258721724.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/2tehz1258721724.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/3nju91258721724.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/4lzrg1258721724.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/5tsz51258721724.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
7.8067892 11.3503361 14.1135184 2.1197630 3.0398315 8.7651692
7 8 9 10 11 12
6.5701821 3.5561873 -0.9221361 -0.9850960 -2.0117855 -6.5173288
13 14 15 16 17 18
-1.6528355 -2.7400842 -3.3206146 -5.8972615 -8.9250972 -11.4739251
19 20 21 22 23 24
-12.0057815 -15.7847893 -24.4741476 -18.5097336 -17.5227362 -13.1606146
25 26 27 28 29 30
-0.5516387 5.7584605 2.9608214 -1.9755343 -3.3444308 -5.1137892
31 32 33 34 35 36
-10.8285369 -8.5433015 -4.3326598 -0.7648241 6.1700776 23.8329689
37 38 39 40 41 42
25.5716279 19.2490493 26.3658669 31.6029070 47.5724193 43.8293223
43 44 45 46 47 48
41.1172266 28.2397583 35.3177221 33.6547623 21.1470637 8.7791594
49 50 51 52 53 54
-30.1529255 -33.6177616 -40.1195922 -25.8498742 -38.3427229 -36.0067772
55 56 57 58 59 60
-24.8530903 -7.4678549 -5.5887786 -13.3951085 -7.7826194 -12.9341848
61
-1.0210174
> postscript(file="/var/www/html/rcomp/tmp/69fxr1258721724.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 7.8067892 NA
1 11.3503361 7.8067892
2 14.1135184 11.3503361
3 2.1197630 14.1135184
4 3.0398315 2.1197630
5 8.7651692 3.0398315
6 6.5701821 8.7651692
7 3.5561873 6.5701821
8 -0.9221361 3.5561873
9 -0.9850960 -0.9221361
10 -2.0117855 -0.9850960
11 -6.5173288 -2.0117855
12 -1.6528355 -6.5173288
13 -2.7400842 -1.6528355
14 -3.3206146 -2.7400842
15 -5.8972615 -3.3206146
16 -8.9250972 -5.8972615
17 -11.4739251 -8.9250972
18 -12.0057815 -11.4739251
19 -15.7847893 -12.0057815
20 -24.4741476 -15.7847893
21 -18.5097336 -24.4741476
22 -17.5227362 -18.5097336
23 -13.1606146 -17.5227362
24 -0.5516387 -13.1606146
25 5.7584605 -0.5516387
26 2.9608214 5.7584605
27 -1.9755343 2.9608214
28 -3.3444308 -1.9755343
29 -5.1137892 -3.3444308
30 -10.8285369 -5.1137892
31 -8.5433015 -10.8285369
32 -4.3326598 -8.5433015
33 -0.7648241 -4.3326598
34 6.1700776 -0.7648241
35 23.8329689 6.1700776
36 25.5716279 23.8329689
37 19.2490493 25.5716279
38 26.3658669 19.2490493
39 31.6029070 26.3658669
40 47.5724193 31.6029070
41 43.8293223 47.5724193
42 41.1172266 43.8293223
43 28.2397583 41.1172266
44 35.3177221 28.2397583
45 33.6547623 35.3177221
46 21.1470637 33.6547623
47 8.7791594 21.1470637
48 -30.1529255 8.7791594
49 -33.6177616 -30.1529255
50 -40.1195922 -33.6177616
51 -25.8498742 -40.1195922
52 -38.3427229 -25.8498742
53 -36.0067772 -38.3427229
54 -24.8530903 -36.0067772
55 -7.4678549 -24.8530903
56 -5.5887786 -7.4678549
57 -13.3951085 -5.5887786
58 -7.7826194 -13.3951085
59 -12.9341848 -7.7826194
60 -1.0210174 -12.9341848
61 NA -1.0210174
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 11.3503361 7.8067892
[2,] 14.1135184 11.3503361
[3,] 2.1197630 14.1135184
[4,] 3.0398315 2.1197630
[5,] 8.7651692 3.0398315
[6,] 6.5701821 8.7651692
[7,] 3.5561873 6.5701821
[8,] -0.9221361 3.5561873
[9,] -0.9850960 -0.9221361
[10,] -2.0117855 -0.9850960
[11,] -6.5173288 -2.0117855
[12,] -1.6528355 -6.5173288
[13,] -2.7400842 -1.6528355
[14,] -3.3206146 -2.7400842
[15,] -5.8972615 -3.3206146
[16,] -8.9250972 -5.8972615
[17,] -11.4739251 -8.9250972
[18,] -12.0057815 -11.4739251
[19,] -15.7847893 -12.0057815
[20,] -24.4741476 -15.7847893
[21,] -18.5097336 -24.4741476
[22,] -17.5227362 -18.5097336
[23,] -13.1606146 -17.5227362
[24,] -0.5516387 -13.1606146
[25,] 5.7584605 -0.5516387
[26,] 2.9608214 5.7584605
[27,] -1.9755343 2.9608214
[28,] -3.3444308 -1.9755343
[29,] -5.1137892 -3.3444308
[30,] -10.8285369 -5.1137892
[31,] -8.5433015 -10.8285369
[32,] -4.3326598 -8.5433015
[33,] -0.7648241 -4.3326598
[34,] 6.1700776 -0.7648241
[35,] 23.8329689 6.1700776
[36,] 25.5716279 23.8329689
[37,] 19.2490493 25.5716279
[38,] 26.3658669 19.2490493
[39,] 31.6029070 26.3658669
[40,] 47.5724193 31.6029070
[41,] 43.8293223 47.5724193
[42,] 41.1172266 43.8293223
[43,] 28.2397583 41.1172266
[44,] 35.3177221 28.2397583
[45,] 33.6547623 35.3177221
[46,] 21.1470637 33.6547623
[47,] 8.7791594 21.1470637
[48,] -30.1529255 8.7791594
[49,] -33.6177616 -30.1529255
[50,] -40.1195922 -33.6177616
[51,] -25.8498742 -40.1195922
[52,] -38.3427229 -25.8498742
[53,] -36.0067772 -38.3427229
[54,] -24.8530903 -36.0067772
[55,] -7.4678549 -24.8530903
[56,] -5.5887786 -7.4678549
[57,] -13.3951085 -5.5887786
[58,] -7.7826194 -13.3951085
[59,] -12.9341848 -7.7826194
[60,] -1.0210174 -12.9341848
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 11.3503361 7.8067892
2 14.1135184 11.3503361
3 2.1197630 14.1135184
4 3.0398315 2.1197630
5 8.7651692 3.0398315
6 6.5701821 8.7651692
7 3.5561873 6.5701821
8 -0.9221361 3.5561873
9 -0.9850960 -0.9221361
10 -2.0117855 -0.9850960
11 -6.5173288 -2.0117855
12 -1.6528355 -6.5173288
13 -2.7400842 -1.6528355
14 -3.3206146 -2.7400842
15 -5.8972615 -3.3206146
16 -8.9250972 -5.8972615
17 -11.4739251 -8.9250972
18 -12.0057815 -11.4739251
19 -15.7847893 -12.0057815
20 -24.4741476 -15.7847893
21 -18.5097336 -24.4741476
22 -17.5227362 -18.5097336
23 -13.1606146 -17.5227362
24 -0.5516387 -13.1606146
25 5.7584605 -0.5516387
26 2.9608214 5.7584605
27 -1.9755343 2.9608214
28 -3.3444308 -1.9755343
29 -5.1137892 -3.3444308
30 -10.8285369 -5.1137892
31 -8.5433015 -10.8285369
32 -4.3326598 -8.5433015
33 -0.7648241 -4.3326598
34 6.1700776 -0.7648241
35 23.8329689 6.1700776
36 25.5716279 23.8329689
37 19.2490493 25.5716279
38 26.3658669 19.2490493
39 31.6029070 26.3658669
40 47.5724193 31.6029070
41 43.8293223 47.5724193
42 41.1172266 43.8293223
43 28.2397583 41.1172266
44 35.3177221 28.2397583
45 33.6547623 35.3177221
46 21.1470637 33.6547623
47 8.7791594 21.1470637
48 -30.1529255 8.7791594
49 -33.6177616 -30.1529255
50 -40.1195922 -33.6177616
51 -25.8498742 -40.1195922
52 -38.3427229 -25.8498742
53 -36.0067772 -38.3427229
54 -24.8530903 -36.0067772
55 -7.4678549 -24.8530903
56 -5.5887786 -7.4678549
57 -13.3951085 -5.5887786
58 -7.7826194 -13.3951085
59 -12.9341848 -7.7826194
60 -1.0210174 -12.9341848
> 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/7w27f1258721724.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/83ep81258721724.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/9pxsn1258721724.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/10c5kz1258721724.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/119ivp1258721724.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/12l1aa1258721724.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/137s2c1258721724.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/14b48r1258721724.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/152be11258721724.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/161amx1258721724.tab")
+ }
>
> system("convert tmp/111xp1258721724.ps tmp/111xp1258721724.png")
> system("convert tmp/2tehz1258721724.ps tmp/2tehz1258721724.png")
> system("convert tmp/3nju91258721724.ps tmp/3nju91258721724.png")
> system("convert tmp/4lzrg1258721724.ps tmp/4lzrg1258721724.png")
> system("convert tmp/5tsz51258721724.ps tmp/5tsz51258721724.png")
> system("convert tmp/69fxr1258721724.ps tmp/69fxr1258721724.png")
> system("convert tmp/7w27f1258721724.ps tmp/7w27f1258721724.png")
> system("convert tmp/83ep81258721724.ps tmp/83ep81258721724.png")
> system("convert tmp/9pxsn1258721724.ps tmp/9pxsn1258721724.png")
> system("convert tmp/10c5kz1258721724.ps tmp/10c5kz1258721724.png")
>
>
> proc.time()
user system elapsed
2.391 1.536 4.782