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 = '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 108.5 98.71 1 0 0 0 0 0 0 0 0 0 0 1
2 112.3 98.54 0 1 0 0 0 0 0 0 0 0 0 2
3 116.6 98.20 0 0 1 0 0 0 0 0 0 0 0 3
4 115.5 96.92 0 0 0 1 0 0 0 0 0 0 0 4
5 120.1 99.06 0 0 0 0 1 0 0 0 0 0 0 5
6 132.9 99.65 0 0 0 0 0 1 0 0 0 0 0 6
7 128.1 99.82 0 0 0 0 0 0 1 0 0 0 0 7
8 129.3 99.99 0 0 0 0 0 0 0 1 0 0 0 8
9 132.5 100.33 0 0 0 0 0 0 0 0 1 0 0 9
10 131.0 99.31 0 0 0 0 0 0 0 0 0 1 0 10
11 124.9 101.10 0 0 0 0 0 0 0 0 0 0 1 11
12 120.8 101.10 0 0 0 0 0 0 0 0 0 0 0 12
13 122.0 100.93 1 0 0 0 0 0 0 0 0 0 0 13
14 122.1 100.85 0 1 0 0 0 0 0 0 0 0 0 14
15 127.4 100.93 0 0 1 0 0 0 0 0 0 0 0 15
16 135.2 99.60 0 0 0 1 0 0 0 0 0 0 0 16
17 137.3 101.88 0 0 0 0 1 0 0 0 0 0 0 17
18 135.0 101.81 0 0 0 0 0 1 0 0 0 0 0 18
19 136.0 102.38 0 0 0 0 0 0 1 0 0 0 0 19
20 138.4 102.74 0 0 0 0 0 0 0 1 0 0 0 20
21 134.7 102.82 0 0 0 0 0 0 0 0 1 0 0 21
22 138.4 101.72 0 0 0 0 0 0 0 0 0 1 0 22
23 133.9 103.47 0 0 0 0 0 0 0 0 0 0 1 23
24 133.6 102.98 0 0 0 0 0 0 0 0 0 0 0 24
25 141.2 102.68 1 0 0 0 0 0 0 0 0 0 0 25
26 151.8 102.90 0 1 0 0 0 0 0 0 0 0 0 26
27 155.4 103.03 0 0 1 0 0 0 0 0 0 0 0 27
28 156.6 101.29 0 0 0 1 0 0 0 0 0 0 0 28
29 161.6 103.69 0 0 0 0 1 0 0 0 0 0 0 29
30 160.7 103.68 0 0 0 0 0 1 0 0 0 0 0 30
31 156.0 104.20 0 0 0 0 0 0 1 0 0 0 0 31
32 159.5 104.08 0 0 0 0 0 0 0 1 0 0 0 32
33 168.7 104.16 0 0 0 0 0 0 0 0 1 0 0 33
34 169.9 103.05 0 0 0 0 0 0 0 0 0 1 0 34
35 169.9 104.66 0 0 0 0 0 0 0 0 0 0 1 35
36 185.9 104.46 0 0 0 0 0 0 0 0 0 0 0 36
37 190.8 104.95 1 0 0 0 0 0 0 0 0 0 0 37
38 195.8 105.85 0 1 0 0 0 0 0 0 0 0 0 38
39 211.9 106.23 0 0 1 0 0 0 0 0 0 0 0 39
40 227.1 104.86 0 0 0 1 0 0 0 0 0 0 0 40
41 251.3 107.44 0 0 0 0 1 0 0 0 0 0 0 41
42 256.7 108.23 0 0 0 0 0 1 0 0 0 0 0 42
43 251.9 108.45 0 0 0 0 0 0 1 0 0 0 0 43
44 251.2 109.39 0 0 0 0 0 0 0 1 0 0 0 44
45 270.3 110.15 0 0 0 0 0 0 0 0 1 0 0 45
46 267.2 109.13 0 0 0 0 0 0 0 0 0 1 0 46
47 243.0 110.28 0 0 0 0 0 0 0 0 0 0 1 47
48 229.9 110.17 0 0 0 0 0 0 0 0 0 0 0 48
49 187.2 109.99 1 0 0 0 0 0 0 0 0 0 0 49
50 178.2 109.26 0 1 0 0 0 0 0 0 0 0 0 50
51 175.2 109.11 0 0 1 0 0 0 0 0 0 0 0 51
52 192.4 107.06 0 0 0 1 0 0 0 0 0 0 0 52
53 187.0 109.53 0 0 0 0 1 0 0 0 0 0 0 53
54 184.0 108.92 0 0 0 0 0 1 0 0 0 0 0 54
55 194.1 109.24 0 0 0 0 0 0 1 0 0 0 0 55
56 212.7 109.12 0 0 0 0 0 0 0 1 0 0 0 56
57 217.5 109.00 0 0 0 0 0 0 0 0 1 0 0 57
58 200.5 107.23 0 0 0 0 0 0 0 0 0 1 0 58
59 205.9 109.49 0 0 0 0 0 0 0 0 0 0 1 59
60 196.5 109.04 0 0 0 0 0 0 0 0 0 0 0 60
61 206.3 109.02 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
-1581.714 17.139 -1.455 -0.824 5.593 41.788
M5 M6 M7 M8 M9 M10
8.700 10.235 4.925 7.209 11.321 30.117
M11 t
-3.605 -1.500
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-42.2458 -10.3430 0.9688 8.9224 44.3682
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1581.7137 282.0909 -5.607 1.05e-06 ***
X 17.1394 2.8725 5.967 3.02e-07 ***
M1 -1.4549 12.9317 -0.113 0.91090
M2 -0.8240 13.5667 -0.061 0.95182
M3 5.5934 13.5609 0.412 0.68187
M4 41.7883 14.6522 2.852 0.00644 **
M5 8.6997 13.5340 0.643 0.52348
M6 10.2347 13.5147 0.757 0.45265
M7 4.9248 13.5255 0.364 0.71741
M8 7.2087 13.5250 0.533 0.59655
M9 11.3212 13.5232 0.837 0.40674
M10 30.1173 13.8120 2.181 0.03426 *
M11 -3.6051 13.5397 -0.266 0.79120
t -1.5003 0.6128 -2.448 0.01814 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.31 on 47 degrees of freedom
Multiple R-squared: 0.8208, Adjusted R-squared: 0.7712
F-statistic: 16.56 on 13 and 47 DF, p-value: 2.463e-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,] 2.124620e-03 4.249241e-03 0.9978754
[2,] 6.927599e-04 1.385520e-03 0.9993072
[3,] 1.358208e-04 2.716416e-04 0.9998642
[4,] 3.011761e-05 6.023523e-05 0.9999699
[5,] 1.512294e-05 3.024588e-05 0.9999849
[6,] 2.661962e-06 5.323923e-06 0.9999973
[7,] 4.898441e-07 9.796883e-07 0.9999995
[8,] 2.758907e-07 5.517813e-07 0.9999997
[9,] 3.177591e-07 6.355183e-07 0.9999997
[10,] 9.416982e-07 1.883396e-06 0.9999991
[11,] 8.366245e-07 1.673249e-06 0.9999992
[12,] 2.459338e-07 4.918676e-07 0.9999998
[13,] 9.530040e-08 1.906008e-07 0.9999999
[14,] 2.246646e-08 4.493292e-08 1.0000000
[15,] 7.937304e-09 1.587461e-08 1.0000000
[16,] 3.404171e-09 6.808342e-09 1.0000000
[17,] 2.609784e-09 5.219567e-09 1.0000000
[18,] 2.880856e-09 5.761711e-09 1.0000000
[19,] 9.180101e-09 1.836020e-08 1.0000000
[20,] 9.741128e-07 1.948226e-06 0.9999990
[21,] 3.740971e-04 7.481942e-04 0.9996259
[22,] 5.835048e-03 1.167010e-02 0.9941650
[23,] 2.015648e-02 4.031296e-02 0.9798435
[24,] 4.903227e-02 9.806453e-02 0.9509677
[25,] 7.008435e-02 1.401687e-01 0.9299157
[26,] 7.249869e-02 1.449974e-01 0.9275013
[27,] 9.275583e-02 1.855117e-01 0.9072442
[28,] 9.919957e-02 1.983991e-01 0.9008004
> postscript(file="/var/www/html/rcomp/tmp/1hvvw1258722821.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/2chfy1258722821.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/3z20t1258722821.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/4cd8s1258722821.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/5yntj1258722821.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
1.3393209 8.9224380 14.1326194 0.2764258 2.7870437 5.4400382
7 8 9 10 11 12
4.5365228 2.0391167 -3.2004953 -4.5141436 -6.0710162 -12.2758644
13 14 15 16 17 18
-5.2070500 -2.8664781 -3.8548417 -7.9540657 -10.3429628 -11.4779691
19 20 21 22 23 24
-13.4372416 -17.9911323 -25.6745022 -20.4169990 -19.6882960 -13.6948417
25 26 27 28 29 30
2.0020938 9.7008478 6.1555146 2.4834417 0.9378175 0.1744476
31 32 33 34 35 36
-6.6278553 -1.8548375 3.3617926 6.2906897 13.9189078 31.2419381
37 38 39 40 41 42
30.6987534 21.1427203 25.8125389 29.7988906 44.3681757 36.1932916
43 44 45 46 47 48
34.4328066 16.8380680 20.2999110 17.3862627 8.6986015 -4.6209135
49 50 51 52 53 54
-41.2807051 -36.8995280 -42.2458312 -24.6046924 -37.7500741 -30.3298083
55 56 57 58 59 60
-18.9042326 0.9687852 5.2132939 1.2541902 3.1418029 -0.6503185
61
12.4475871
> postscript(file="/var/www/html/rcomp/tmp/6yjr01258722821.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 1.3393209 NA
1 8.9224380 1.3393209
2 14.1326194 8.9224380
3 0.2764258 14.1326194
4 2.7870437 0.2764258
5 5.4400382 2.7870437
6 4.5365228 5.4400382
7 2.0391167 4.5365228
8 -3.2004953 2.0391167
9 -4.5141436 -3.2004953
10 -6.0710162 -4.5141436
11 -12.2758644 -6.0710162
12 -5.2070500 -12.2758644
13 -2.8664781 -5.2070500
14 -3.8548417 -2.8664781
15 -7.9540657 -3.8548417
16 -10.3429628 -7.9540657
17 -11.4779691 -10.3429628
18 -13.4372416 -11.4779691
19 -17.9911323 -13.4372416
20 -25.6745022 -17.9911323
21 -20.4169990 -25.6745022
22 -19.6882960 -20.4169990
23 -13.6948417 -19.6882960
24 2.0020938 -13.6948417
25 9.7008478 2.0020938
26 6.1555146 9.7008478
27 2.4834417 6.1555146
28 0.9378175 2.4834417
29 0.1744476 0.9378175
30 -6.6278553 0.1744476
31 -1.8548375 -6.6278553
32 3.3617926 -1.8548375
33 6.2906897 3.3617926
34 13.9189078 6.2906897
35 31.2419381 13.9189078
36 30.6987534 31.2419381
37 21.1427203 30.6987534
38 25.8125389 21.1427203
39 29.7988906 25.8125389
40 44.3681757 29.7988906
41 36.1932916 44.3681757
42 34.4328066 36.1932916
43 16.8380680 34.4328066
44 20.2999110 16.8380680
45 17.3862627 20.2999110
46 8.6986015 17.3862627
47 -4.6209135 8.6986015
48 -41.2807051 -4.6209135
49 -36.8995280 -41.2807051
50 -42.2458312 -36.8995280
51 -24.6046924 -42.2458312
52 -37.7500741 -24.6046924
53 -30.3298083 -37.7500741
54 -18.9042326 -30.3298083
55 0.9687852 -18.9042326
56 5.2132939 0.9687852
57 1.2541902 5.2132939
58 3.1418029 1.2541902
59 -0.6503185 3.1418029
60 12.4475871 -0.6503185
61 NA 12.4475871
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 8.9224380 1.3393209
[2,] 14.1326194 8.9224380
[3,] 0.2764258 14.1326194
[4,] 2.7870437 0.2764258
[5,] 5.4400382 2.7870437
[6,] 4.5365228 5.4400382
[7,] 2.0391167 4.5365228
[8,] -3.2004953 2.0391167
[9,] -4.5141436 -3.2004953
[10,] -6.0710162 -4.5141436
[11,] -12.2758644 -6.0710162
[12,] -5.2070500 -12.2758644
[13,] -2.8664781 -5.2070500
[14,] -3.8548417 -2.8664781
[15,] -7.9540657 -3.8548417
[16,] -10.3429628 -7.9540657
[17,] -11.4779691 -10.3429628
[18,] -13.4372416 -11.4779691
[19,] -17.9911323 -13.4372416
[20,] -25.6745022 -17.9911323
[21,] -20.4169990 -25.6745022
[22,] -19.6882960 -20.4169990
[23,] -13.6948417 -19.6882960
[24,] 2.0020938 -13.6948417
[25,] 9.7008478 2.0020938
[26,] 6.1555146 9.7008478
[27,] 2.4834417 6.1555146
[28,] 0.9378175 2.4834417
[29,] 0.1744476 0.9378175
[30,] -6.6278553 0.1744476
[31,] -1.8548375 -6.6278553
[32,] 3.3617926 -1.8548375
[33,] 6.2906897 3.3617926
[34,] 13.9189078 6.2906897
[35,] 31.2419381 13.9189078
[36,] 30.6987534 31.2419381
[37,] 21.1427203 30.6987534
[38,] 25.8125389 21.1427203
[39,] 29.7988906 25.8125389
[40,] 44.3681757 29.7988906
[41,] 36.1932916 44.3681757
[42,] 34.4328066 36.1932916
[43,] 16.8380680 34.4328066
[44,] 20.2999110 16.8380680
[45,] 17.3862627 20.2999110
[46,] 8.6986015 17.3862627
[47,] -4.6209135 8.6986015
[48,] -41.2807051 -4.6209135
[49,] -36.8995280 -41.2807051
[50,] -42.2458312 -36.8995280
[51,] -24.6046924 -42.2458312
[52,] -37.7500741 -24.6046924
[53,] -30.3298083 -37.7500741
[54,] -18.9042326 -30.3298083
[55,] 0.9687852 -18.9042326
[56,] 5.2132939 0.9687852
[57,] 1.2541902 5.2132939
[58,] 3.1418029 1.2541902
[59,] -0.6503185 3.1418029
[60,] 12.4475871 -0.6503185
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 8.9224380 1.3393209
2 14.1326194 8.9224380
3 0.2764258 14.1326194
4 2.7870437 0.2764258
5 5.4400382 2.7870437
6 4.5365228 5.4400382
7 2.0391167 4.5365228
8 -3.2004953 2.0391167
9 -4.5141436 -3.2004953
10 -6.0710162 -4.5141436
11 -12.2758644 -6.0710162
12 -5.2070500 -12.2758644
13 -2.8664781 -5.2070500
14 -3.8548417 -2.8664781
15 -7.9540657 -3.8548417
16 -10.3429628 -7.9540657
17 -11.4779691 -10.3429628
18 -13.4372416 -11.4779691
19 -17.9911323 -13.4372416
20 -25.6745022 -17.9911323
21 -20.4169990 -25.6745022
22 -19.6882960 -20.4169990
23 -13.6948417 -19.6882960
24 2.0020938 -13.6948417
25 9.7008478 2.0020938
26 6.1555146 9.7008478
27 2.4834417 6.1555146
28 0.9378175 2.4834417
29 0.1744476 0.9378175
30 -6.6278553 0.1744476
31 -1.8548375 -6.6278553
32 3.3617926 -1.8548375
33 6.2906897 3.3617926
34 13.9189078 6.2906897
35 31.2419381 13.9189078
36 30.6987534 31.2419381
37 21.1427203 30.6987534
38 25.8125389 21.1427203
39 29.7988906 25.8125389
40 44.3681757 29.7988906
41 36.1932916 44.3681757
42 34.4328066 36.1932916
43 16.8380680 34.4328066
44 20.2999110 16.8380680
45 17.3862627 20.2999110
46 8.6986015 17.3862627
47 -4.6209135 8.6986015
48 -41.2807051 -4.6209135
49 -36.8995280 -41.2807051
50 -42.2458312 -36.8995280
51 -24.6046924 -42.2458312
52 -37.7500741 -24.6046924
53 -30.3298083 -37.7500741
54 -18.9042326 -30.3298083
55 0.9687852 -18.9042326
56 5.2132939 0.9687852
57 1.2541902 5.2132939
58 3.1418029 1.2541902
59 -0.6503185 3.1418029
60 12.4475871 -0.6503185
> 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/7g1x61258722821.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/8rd2x1258722821.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/9vw3k1258722821.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/1050qa1258722821.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/11yzwh1258722821.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/12wcm61258722821.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/13o40t1258722821.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/14p69r1258722821.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/15jzsy1258722821.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/161l301258722821.tab")
+ }
>
> system("convert tmp/1hvvw1258722821.ps tmp/1hvvw1258722821.png")
> system("convert tmp/2chfy1258722821.ps tmp/2chfy1258722821.png")
> system("convert tmp/3z20t1258722821.ps tmp/3z20t1258722821.png")
> system("convert tmp/4cd8s1258722821.ps tmp/4cd8s1258722821.png")
> system("convert tmp/5yntj1258722821.ps tmp/5yntj1258722821.png")
> system("convert tmp/6yjr01258722821.ps tmp/6yjr01258722821.png")
> system("convert tmp/7g1x61258722821.ps tmp/7g1x61258722821.png")
> system("convert tmp/8rd2x1258722821.ps tmp/8rd2x1258722821.png")
> system("convert tmp/9vw3k1258722821.ps tmp/9vw3k1258722821.png")
> system("convert tmp/1050qa1258722821.ps tmp/1050qa1258722821.png")
>
>
> proc.time()
user system elapsed
2.345 1.490 2.769