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(115.6,37.2,111.9,37.2,107,34.7,107.1,32.5,100.6,33.5,99.2,31.5,108.4,31.2,103,27,99.8,26.7,115,26.5,90.8,26,95.9,27.2,114.4,30.5,108.2,33.7,112.6,34.2,109.1,36.7,105,36.2,105,38.5,118.5,40,103.7,42.5,112.5,43.5,116.6,43.3,96.6,45.5,101.9,44.3,116.5,43,119.3,43.5,115.4,41.5,108.5,42.5,111.5,41.3,108.8,39.5,121.8,38.5,109.6,41,112.2,44.5,119.6,46,104.1,44,105.3,41.5,115,41.3,124.1,38,116.8,38,107.5,36.2,115.6,38.7,116.2,38.7,116.3,39.2,119,35.7,111.9,36.5,118.6,36.7,106.9,34.7,103.2,35,118.6,28.2,118.7,23.7,102.8,15,100.6,8.7,94.9,11,94.5,7.5,102.9,5.7,95.3,9.3,92.5,10.2,102.7,15.7,91.5,18.1,89.5,20.8),dim=c(2,60),dimnames=list(c('Ipzb','Cvn'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Ipzb','Cvn'),1:60))
> 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
Ipzb Cvn M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 115.6 37.2 1 0 0 0 0 0 0 0 0 0 0 1
2 111.9 37.2 0 1 0 0 0 0 0 0 0 0 0 2
3 107.0 34.7 0 0 1 0 0 0 0 0 0 0 0 3
4 107.1 32.5 0 0 0 1 0 0 0 0 0 0 0 4
5 100.6 33.5 0 0 0 0 1 0 0 0 0 0 0 5
6 99.2 31.5 0 0 0 0 0 1 0 0 0 0 0 6
7 108.4 31.2 0 0 0 0 0 0 1 0 0 0 0 7
8 103.0 27.0 0 0 0 0 0 0 0 1 0 0 0 8
9 99.8 26.7 0 0 0 0 0 0 0 0 1 0 0 9
10 115.0 26.5 0 0 0 0 0 0 0 0 0 1 0 10
11 90.8 26.0 0 0 0 0 0 0 0 0 0 0 1 11
12 95.9 27.2 0 0 0 0 0 0 0 0 0 0 0 12
13 114.4 30.5 1 0 0 0 0 0 0 0 0 0 0 13
14 108.2 33.7 0 1 0 0 0 0 0 0 0 0 0 14
15 112.6 34.2 0 0 1 0 0 0 0 0 0 0 0 15
16 109.1 36.7 0 0 0 1 0 0 0 0 0 0 0 16
17 105.0 36.2 0 0 0 0 1 0 0 0 0 0 0 17
18 105.0 38.5 0 0 0 0 0 1 0 0 0 0 0 18
19 118.5 40.0 0 0 0 0 0 0 1 0 0 0 0 19
20 103.7 42.5 0 0 0 0 0 0 0 1 0 0 0 20
21 112.5 43.5 0 0 0 0 0 0 0 0 1 0 0 21
22 116.6 43.3 0 0 0 0 0 0 0 0 0 1 0 22
23 96.6 45.5 0 0 0 0 0 0 0 0 0 0 1 23
24 101.9 44.3 0 0 0 0 0 0 0 0 0 0 0 24
25 116.5 43.0 1 0 0 0 0 0 0 0 0 0 0 25
26 119.3 43.5 0 1 0 0 0 0 0 0 0 0 0 26
27 115.4 41.5 0 0 1 0 0 0 0 0 0 0 0 27
28 108.5 42.5 0 0 0 1 0 0 0 0 0 0 0 28
29 111.5 41.3 0 0 0 0 1 0 0 0 0 0 0 29
30 108.8 39.5 0 0 0 0 0 1 0 0 0 0 0 30
31 121.8 38.5 0 0 0 0 0 0 1 0 0 0 0 31
32 109.6 41.0 0 0 0 0 0 0 0 1 0 0 0 32
33 112.2 44.5 0 0 0 0 0 0 0 0 1 0 0 33
34 119.6 46.0 0 0 0 0 0 0 0 0 0 1 0 34
35 104.1 44.0 0 0 0 0 0 0 0 0 0 0 1 35
36 105.3 41.5 0 0 0 0 0 0 0 0 0 0 0 36
37 115.0 41.3 1 0 0 0 0 0 0 0 0 0 0 37
38 124.1 38.0 0 1 0 0 0 0 0 0 0 0 0 38
39 116.8 38.0 0 0 1 0 0 0 0 0 0 0 0 39
40 107.5 36.2 0 0 0 1 0 0 0 0 0 0 0 40
41 115.6 38.7 0 0 0 0 1 0 0 0 0 0 0 41
42 116.2 38.7 0 0 0 0 0 1 0 0 0 0 0 42
43 116.3 39.2 0 0 0 0 0 0 1 0 0 0 0 43
44 119.0 35.7 0 0 0 0 0 0 0 1 0 0 0 44
45 111.9 36.5 0 0 0 0 0 0 0 0 1 0 0 45
46 118.6 36.7 0 0 0 0 0 0 0 0 0 1 0 46
47 106.9 34.7 0 0 0 0 0 0 0 0 0 0 1 47
48 103.2 35.0 0 0 0 0 0 0 0 0 0 0 0 48
49 118.6 28.2 1 0 0 0 0 0 0 0 0 0 0 49
50 118.7 23.7 0 1 0 0 0 0 0 0 0 0 0 50
51 102.8 15.0 0 0 1 0 0 0 0 0 0 0 0 51
52 100.6 8.7 0 0 0 1 0 0 0 0 0 0 0 52
53 94.9 11.0 0 0 0 0 1 0 0 0 0 0 0 53
54 94.5 7.5 0 0 0 0 0 1 0 0 0 0 0 54
55 102.9 5.7 0 0 0 0 0 0 1 0 0 0 0 55
56 95.3 9.3 0 0 0 0 0 0 0 1 0 0 0 56
57 92.5 10.2 0 0 0 0 0 0 0 0 1 0 0 57
58 102.7 15.7 0 0 0 0 0 0 0 0 0 1 0 58
59 91.5 18.1 0 0 0 0 0 0 0 0 0 0 1 59
60 89.5 20.8 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Cvn M1 M2 M3 M4
75.7410 0.5492 17.0990 17.8337 13.5730 9.8243
M5 M6 M7 M8 M9 M10
8.1985 7.8321 16.6573 8.9629 7.8394 15.6770
M11 t
-0.9895 0.1355
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-7.0538 -2.1824 -0.3271 2.2432 8.7273
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 75.74099 3.20880 23.604 < 2e-16 ***
Cvn 0.54915 0.05486 10.009 3.95e-13 ***
M1 17.09895 2.56235 6.673 2.82e-08 ***
M2 17.83371 2.55925 6.968 1.02e-08 ***
M3 13.57301 2.56253 5.297 3.23e-06 ***
M4 9.82430 2.56473 3.831 0.000386 ***
M5 8.19845 2.55654 3.207 0.002443 **
M6 7.83206 2.55751 3.062 0.003662 **
M7 16.65732 2.55516 6.519 4.81e-08 ***
M8 8.96293 2.55128 3.513 0.001005 **
M9 7.83938 2.54466 3.081 0.003481 **
M10 15.67699 2.54078 6.170 1.61e-07 ***
M11 -0.98954 2.54003 -0.390 0.698647
t 0.13555 0.03439 3.941 0.000274 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.016 on 46 degrees of freedom
Multiple R-squared: 0.8345, Adjusted R-squared: 0.7877
F-statistic: 17.84 on 13 and 46 DF, p-value: 9.158e-14
> 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.24348146 0.48696292 0.7565185
[2,] 0.11635796 0.23271592 0.8836420
[3,] 0.07480164 0.14960328 0.9251984
[4,] 0.27704366 0.55408731 0.7229563
[5,] 0.27981146 0.55962292 0.7201885
[6,] 0.32786340 0.65572680 0.6721366
[7,] 0.26712195 0.53424391 0.7328780
[8,] 0.18560903 0.37121806 0.8143910
[9,] 0.15242089 0.30484179 0.8475791
[10,] 0.16680795 0.33361591 0.8331920
[11,] 0.10861491 0.21722982 0.8913851
[12,] 0.11733818 0.23467635 0.8826618
[13,] 0.10036550 0.20073100 0.8996345
[14,] 0.07592669 0.15185338 0.9240733
[15,] 0.07702871 0.15405741 0.9229713
[16,] 0.06863893 0.13727786 0.9313611
[17,] 0.04113564 0.08227128 0.9588644
[18,] 0.02693971 0.05387943 0.9730603
[19,] 0.02973090 0.05946180 0.9702691
[20,] 0.01623851 0.03247702 0.9837615
[21,] 0.13221576 0.26443153 0.8677842
[22,] 0.23465199 0.46930399 0.7653480
[23,] 0.16639126 0.33278253 0.8336087
[24,] 0.60813129 0.78373742 0.3918687
[25,] 0.52850949 0.94298101 0.4714905
[26,] 0.48265385 0.96530771 0.5173461
[27,] 0.89742672 0.20514657 0.1025733
> postscript(file="/var/www/html/rcomp/tmp/15yn91258728176.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/26q0c1258728176.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/3ijo81258728176.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/4ilb11258728176.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/5sj2y1258728176.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 = 60
Frequency = 1
1 2 3 4 5 6
2.19611274 -2.37419070 -1.77615672 3.14512973 -2.41371736 -2.48456683
7 8 9 10 11 12
-2.08063479 2.38464455 0.33738734 7.67406218 0.27962046 3.59555487
13 14 15 16 17 18
3.04885764 -5.77872751 2.47185487 1.21213381 -1.12298749 -2.15518425
19 20 21 22 23 24
1.56027683 -7.05375241 2.18509468 -1.57823048 -6.25537865 -1.42148295
25 26 27 28 29 30
-3.34208772 -1.68696643 -0.36350771 -4.19950298 0.94978110 -0.53089847
31 32 33 34 35 36
4.05743895 -1.95659029 -0.29061953 -1.68750061 0.44178347 1.88957486
37 38 39 40 41 42
-5.53509549 4.50679783 1.33195548 -3.36641829 4.85100881 5.68185828
43 44 45 46 47 48
-3.45353011 8.72734386 2.17602106 0.79303568 6.72231976 1.73248966
49 50 51 52 53 54
3.63221283 5.33308680 -1.66414591 3.20865773 -2.26408506 -0.51120873
55 56 57 58 59 60
-0.08355088 -2.10164571 -4.40788356 -5.20136677 -1.18834504 -5.79613643
> postscript(file="/var/www/html/rcomp/tmp/62s7y1258728176.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 2.19611274 NA
1 -2.37419070 2.19611274
2 -1.77615672 -2.37419070
3 3.14512973 -1.77615672
4 -2.41371736 3.14512973
5 -2.48456683 -2.41371736
6 -2.08063479 -2.48456683
7 2.38464455 -2.08063479
8 0.33738734 2.38464455
9 7.67406218 0.33738734
10 0.27962046 7.67406218
11 3.59555487 0.27962046
12 3.04885764 3.59555487
13 -5.77872751 3.04885764
14 2.47185487 -5.77872751
15 1.21213381 2.47185487
16 -1.12298749 1.21213381
17 -2.15518425 -1.12298749
18 1.56027683 -2.15518425
19 -7.05375241 1.56027683
20 2.18509468 -7.05375241
21 -1.57823048 2.18509468
22 -6.25537865 -1.57823048
23 -1.42148295 -6.25537865
24 -3.34208772 -1.42148295
25 -1.68696643 -3.34208772
26 -0.36350771 -1.68696643
27 -4.19950298 -0.36350771
28 0.94978110 -4.19950298
29 -0.53089847 0.94978110
30 4.05743895 -0.53089847
31 -1.95659029 4.05743895
32 -0.29061953 -1.95659029
33 -1.68750061 -0.29061953
34 0.44178347 -1.68750061
35 1.88957486 0.44178347
36 -5.53509549 1.88957486
37 4.50679783 -5.53509549
38 1.33195548 4.50679783
39 -3.36641829 1.33195548
40 4.85100881 -3.36641829
41 5.68185828 4.85100881
42 -3.45353011 5.68185828
43 8.72734386 -3.45353011
44 2.17602106 8.72734386
45 0.79303568 2.17602106
46 6.72231976 0.79303568
47 1.73248966 6.72231976
48 3.63221283 1.73248966
49 5.33308680 3.63221283
50 -1.66414591 5.33308680
51 3.20865773 -1.66414591
52 -2.26408506 3.20865773
53 -0.51120873 -2.26408506
54 -0.08355088 -0.51120873
55 -2.10164571 -0.08355088
56 -4.40788356 -2.10164571
57 -5.20136677 -4.40788356
58 -1.18834504 -5.20136677
59 -5.79613643 -1.18834504
60 NA -5.79613643
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.37419070 2.19611274
[2,] -1.77615672 -2.37419070
[3,] 3.14512973 -1.77615672
[4,] -2.41371736 3.14512973
[5,] -2.48456683 -2.41371736
[6,] -2.08063479 -2.48456683
[7,] 2.38464455 -2.08063479
[8,] 0.33738734 2.38464455
[9,] 7.67406218 0.33738734
[10,] 0.27962046 7.67406218
[11,] 3.59555487 0.27962046
[12,] 3.04885764 3.59555487
[13,] -5.77872751 3.04885764
[14,] 2.47185487 -5.77872751
[15,] 1.21213381 2.47185487
[16,] -1.12298749 1.21213381
[17,] -2.15518425 -1.12298749
[18,] 1.56027683 -2.15518425
[19,] -7.05375241 1.56027683
[20,] 2.18509468 -7.05375241
[21,] -1.57823048 2.18509468
[22,] -6.25537865 -1.57823048
[23,] -1.42148295 -6.25537865
[24,] -3.34208772 -1.42148295
[25,] -1.68696643 -3.34208772
[26,] -0.36350771 -1.68696643
[27,] -4.19950298 -0.36350771
[28,] 0.94978110 -4.19950298
[29,] -0.53089847 0.94978110
[30,] 4.05743895 -0.53089847
[31,] -1.95659029 4.05743895
[32,] -0.29061953 -1.95659029
[33,] -1.68750061 -0.29061953
[34,] 0.44178347 -1.68750061
[35,] 1.88957486 0.44178347
[36,] -5.53509549 1.88957486
[37,] 4.50679783 -5.53509549
[38,] 1.33195548 4.50679783
[39,] -3.36641829 1.33195548
[40,] 4.85100881 -3.36641829
[41,] 5.68185828 4.85100881
[42,] -3.45353011 5.68185828
[43,] 8.72734386 -3.45353011
[44,] 2.17602106 8.72734386
[45,] 0.79303568 2.17602106
[46,] 6.72231976 0.79303568
[47,] 1.73248966 6.72231976
[48,] 3.63221283 1.73248966
[49,] 5.33308680 3.63221283
[50,] -1.66414591 5.33308680
[51,] 3.20865773 -1.66414591
[52,] -2.26408506 3.20865773
[53,] -0.51120873 -2.26408506
[54,] -0.08355088 -0.51120873
[55,] -2.10164571 -0.08355088
[56,] -4.40788356 -2.10164571
[57,] -5.20136677 -4.40788356
[58,] -1.18834504 -5.20136677
[59,] -5.79613643 -1.18834504
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.37419070 2.19611274
2 -1.77615672 -2.37419070
3 3.14512973 -1.77615672
4 -2.41371736 3.14512973
5 -2.48456683 -2.41371736
6 -2.08063479 -2.48456683
7 2.38464455 -2.08063479
8 0.33738734 2.38464455
9 7.67406218 0.33738734
10 0.27962046 7.67406218
11 3.59555487 0.27962046
12 3.04885764 3.59555487
13 -5.77872751 3.04885764
14 2.47185487 -5.77872751
15 1.21213381 2.47185487
16 -1.12298749 1.21213381
17 -2.15518425 -1.12298749
18 1.56027683 -2.15518425
19 -7.05375241 1.56027683
20 2.18509468 -7.05375241
21 -1.57823048 2.18509468
22 -6.25537865 -1.57823048
23 -1.42148295 -6.25537865
24 -3.34208772 -1.42148295
25 -1.68696643 -3.34208772
26 -0.36350771 -1.68696643
27 -4.19950298 -0.36350771
28 0.94978110 -4.19950298
29 -0.53089847 0.94978110
30 4.05743895 -0.53089847
31 -1.95659029 4.05743895
32 -0.29061953 -1.95659029
33 -1.68750061 -0.29061953
34 0.44178347 -1.68750061
35 1.88957486 0.44178347
36 -5.53509549 1.88957486
37 4.50679783 -5.53509549
38 1.33195548 4.50679783
39 -3.36641829 1.33195548
40 4.85100881 -3.36641829
41 5.68185828 4.85100881
42 -3.45353011 5.68185828
43 8.72734386 -3.45353011
44 2.17602106 8.72734386
45 0.79303568 2.17602106
46 6.72231976 0.79303568
47 1.73248966 6.72231976
48 3.63221283 1.73248966
49 5.33308680 3.63221283
50 -1.66414591 5.33308680
51 3.20865773 -1.66414591
52 -2.26408506 3.20865773
53 -0.51120873 -2.26408506
54 -0.08355088 -0.51120873
55 -2.10164571 -0.08355088
56 -4.40788356 -2.10164571
57 -5.20136677 -4.40788356
58 -1.18834504 -5.20136677
59 -5.79613643 -1.18834504
> 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/7jrqv1258728176.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/8hli01258728176.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/95nma1258728176.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/10kyfe1258728176.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/110cos1258728176.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/12p51m1258728177.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/13arnd1258728177.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/14f1ce1258728177.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/15ofbj1258728177.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/16x9q81258728177.tab")
+ }
>
> system("convert tmp/15yn91258728176.ps tmp/15yn91258728176.png")
> system("convert tmp/26q0c1258728176.ps tmp/26q0c1258728176.png")
> system("convert tmp/3ijo81258728176.ps tmp/3ijo81258728176.png")
> system("convert tmp/4ilb11258728176.ps tmp/4ilb11258728176.png")
> system("convert tmp/5sj2y1258728176.ps tmp/5sj2y1258728176.png")
> system("convert tmp/62s7y1258728176.ps tmp/62s7y1258728176.png")
> system("convert tmp/7jrqv1258728176.ps tmp/7jrqv1258728176.png")
> system("convert tmp/8hli01258728176.ps tmp/8hli01258728176.png")
> system("convert tmp/95nma1258728176.ps tmp/95nma1258728176.png")
> system("convert tmp/10kyfe1258728176.ps tmp/10kyfe1258728176.png")
>
>
> proc.time()
user system elapsed
2.422 1.575 2.864