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(130,0,136.7,0,138.1,0,139.5,0,140.4,0,144.6,0,151.4,0,147.9,0,141.5,0,143.8,0,143.6,0,150.5,0,150.1,0,154.9,0,162.1,0,176.7,0,186.6,0,194.8,0,196.3,0,228.8,0,267.2,0,237.2,0,254.7,0,258.2,0,257.9,0,269.6,0,266.9,0,269.6,0,253.9,0,258.6,0,274.2,0,301.5,0,304.5,0,285.1,0,287.7,0,265.5,0,264.1,0,276.1,0,258.9,0,239.1,0,250.1,1,276.8,1,297.6,1,295.4,1,283,1,275.8,1,279.7,1,254.6,1,234.6,1,176.9,1,148.1,1,122.7,1,124.9,1,121.6,1,128.4,1,144.5,1,151.8,1,167.1,1,173.8,1,203.7,1,199.8,1),dim=c(2,61),dimnames=list(c('Y(t)','X(t)'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Y(t)','X(t)'),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(t) X(t) M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 130.0 0 1 0 0 0 0 0 0 0 0 0 0 1
2 136.7 0 0 1 0 0 0 0 0 0 0 0 0 2
3 138.1 0 0 0 1 0 0 0 0 0 0 0 0 3
4 139.5 0 0 0 0 1 0 0 0 0 0 0 0 4
5 140.4 0 0 0 0 0 1 0 0 0 0 0 0 5
6 144.6 0 0 0 0 0 0 1 0 0 0 0 0 6
7 151.4 0 0 0 0 0 0 0 1 0 0 0 0 7
8 147.9 0 0 0 0 0 0 0 0 1 0 0 0 8
9 141.5 0 0 0 0 0 0 0 0 0 1 0 0 9
10 143.8 0 0 0 0 0 0 0 0 0 0 1 0 10
11 143.6 0 0 0 0 0 0 0 0 0 0 0 1 11
12 150.5 0 0 0 0 0 0 0 0 0 0 0 0 12
13 150.1 0 1 0 0 0 0 0 0 0 0 0 0 13
14 154.9 0 0 1 0 0 0 0 0 0 0 0 0 14
15 162.1 0 0 0 1 0 0 0 0 0 0 0 0 15
16 176.7 0 0 0 0 1 0 0 0 0 0 0 0 16
17 186.6 0 0 0 0 0 1 0 0 0 0 0 0 17
18 194.8 0 0 0 0 0 0 1 0 0 0 0 0 18
19 196.3 0 0 0 0 0 0 0 1 0 0 0 0 19
20 228.8 0 0 0 0 0 0 0 0 1 0 0 0 20
21 267.2 0 0 0 0 0 0 0 0 0 1 0 0 21
22 237.2 0 0 0 0 0 0 0 0 0 0 1 0 22
23 254.7 0 0 0 0 0 0 0 0 0 0 0 1 23
24 258.2 0 0 0 0 0 0 0 0 0 0 0 0 24
25 257.9 0 1 0 0 0 0 0 0 0 0 0 0 25
26 269.6 0 0 1 0 0 0 0 0 0 0 0 0 26
27 266.9 0 0 0 1 0 0 0 0 0 0 0 0 27
28 269.6 0 0 0 0 1 0 0 0 0 0 0 0 28
29 253.9 0 0 0 0 0 1 0 0 0 0 0 0 29
30 258.6 0 0 0 0 0 0 1 0 0 0 0 0 30
31 274.2 0 0 0 0 0 0 0 1 0 0 0 0 31
32 301.5 0 0 0 0 0 0 0 0 1 0 0 0 32
33 304.5 0 0 0 0 0 0 0 0 0 1 0 0 33
34 285.1 0 0 0 0 0 0 0 0 0 0 1 0 34
35 287.7 0 0 0 0 0 0 0 0 0 0 0 1 35
36 265.5 0 0 0 0 0 0 0 0 0 0 0 0 36
37 264.1 0 1 0 0 0 0 0 0 0 0 0 0 37
38 276.1 0 0 1 0 0 0 0 0 0 0 0 0 38
39 258.9 0 0 0 1 0 0 0 0 0 0 0 0 39
40 239.1 0 0 0 0 1 0 0 0 0 0 0 0 40
41 250.1 1 0 0 0 0 1 0 0 0 0 0 0 41
42 276.8 1 0 0 0 0 0 1 0 0 0 0 0 42
43 297.6 1 0 0 0 0 0 0 1 0 0 0 0 43
44 295.4 1 0 0 0 0 0 0 0 1 0 0 0 44
45 283.0 1 0 0 0 0 0 0 0 0 1 0 0 45
46 275.8 1 0 0 0 0 0 0 0 0 0 1 0 46
47 279.7 1 0 0 0 0 0 0 0 0 0 0 1 47
48 254.6 1 0 0 0 0 0 0 0 0 0 0 0 48
49 234.6 1 1 0 0 0 0 0 0 0 0 0 0 49
50 176.9 1 0 1 0 0 0 0 0 0 0 0 0 50
51 148.1 1 0 0 1 0 0 0 0 0 0 0 0 51
52 122.7 1 0 0 0 1 0 0 0 0 0 0 0 52
53 124.9 1 0 0 0 0 1 0 0 0 0 0 0 53
54 121.6 1 0 0 0 0 0 1 0 0 0 0 0 54
55 128.4 1 0 0 0 0 0 0 1 0 0 0 0 55
56 144.5 1 0 0 0 0 0 0 0 1 0 0 0 56
57 151.8 1 0 0 0 0 0 0 0 0 1 0 0 57
58 167.1 1 0 0 0 0 0 0 0 0 0 1 0 58
59 173.8 1 0 0 0 0 0 0 0 0 0 0 1 59
60 203.7 1 0 0 0 0 0 0 0 0 0 0 0 60
61 199.8 1 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(t)` M1 M2 M3 M4
157.121 -106.098 -11.960 -13.819 -24.945 -33.351
M5 M6 M7 M8 M9 M10
-13.577 -8.584 -1.390 9.544 12.418 1.512
M11 t
4.506 3.106
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-92.067 -33.687 -3.394 32.434 114.406
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 157.1206 31.4916 4.989 8.72e-06 ***
`X(t)` -106.0984 27.4671 -3.863 0.000342 ***
M1 -11.9595 34.4406 -0.347 0.729953
M2 -13.8189 36.1315 -0.382 0.703843
M3 -24.9450 36.0802 -0.691 0.492730
M4 -33.3511 36.0440 -0.925 0.359544
M5 -13.5775 36.2550 -0.374 0.709718
M6 -8.5835 36.1567 -0.237 0.813379
M7 -1.3896 36.0732 -0.039 0.969435
M8 9.5443 36.0048 0.265 0.792104
M9 12.4182 35.9515 0.345 0.731323
M10 1.5122 35.9134 0.042 0.966593
M11 4.5061 35.8905 0.126 0.900623
t 3.1061 0.7402 4.196 0.000119 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 56.74 on 47 degrees of freedom
Multiple R-squared: 0.3187, Adjusted R-squared: 0.1302
F-statistic: 1.691 on 13 and 47 DF, p-value: 0.09454
> 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.009449635 0.018899270 0.9905504
[2,] 0.004246976 0.008493952 0.9957530
[3,] 0.001364077 0.002728154 0.9986359
[4,] 0.005813181 0.011626363 0.9941868
[5,] 0.057740473 0.115480946 0.9422595
[6,] 0.065640313 0.131280625 0.9343597
[7,] 0.103901894 0.207803787 0.8960981
[8,] 0.152827506 0.305655012 0.8471725
[9,] 0.217631278 0.435262557 0.7823687
[10,] 0.214596558 0.429193116 0.7854034
[11,] 0.181288640 0.362577280 0.8187114
[12,] 0.129820780 0.259641559 0.8701792
[13,] 0.090067781 0.180135562 0.9099322
[14,] 0.063119264 0.126238529 0.9368807
[15,] 0.040813145 0.081626291 0.9591869
[16,] 0.023963760 0.047927520 0.9760362
[17,] 0.012858991 0.025717982 0.9871410
[18,] 0.007786818 0.015573636 0.9922132
[19,] 0.005470389 0.010940778 0.9945296
[20,] 0.017430462 0.034860924 0.9825695
[21,] 0.157762623 0.315525245 0.8422374
[22,] 0.133540245 0.267080489 0.8664598
[23,] 0.130308785 0.260617569 0.8696912
[24,] 0.158522298 0.317044596 0.8414777
[25,] 0.098007904 0.196015808 0.9019921
[26,] 0.083604672 0.167209343 0.9163953
[27,] 0.128689132 0.257378264 0.8713109
[28,] 0.172350487 0.344700973 0.8276495
> postscript(file="/var/www/html/rcomp/tmp/1vr4d1259397191.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/227z91259397191.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/323h71259397191.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/4c9r81259397191.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/5djpk1259397191.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
-18.26718750 -12.81385417 -3.39385417 3.30614583 -18.67354167 -22.57354167
7 8 9 10 11 12
-26.07354167 -43.61354167 -55.99354167 -45.89354167 -52.19354167 -43.89354167
13 14 15 16 17 18
-35.44010417 -31.88677083 -16.66677083 3.23322917 -9.74645833 -9.64645833
19 20 21 22 23 24
-18.44645833 0.01354167 32.43354167 10.23354167 21.63354167 26.53354167
25 26 27 28 29 30
35.08697917 45.54031250 50.86031250 58.86031250 20.28062500 16.88062500
31 32 33 34 35 36
22.18062500 35.44062500 32.46062500 20.86062500 17.36062500 -3.43937500
37 38 39 40 41 42
4.01406250 14.76739583 5.58739583 -8.91260417 85.30614583 103.90614583
43 44 45 46 47 48
114.40614583 98.16614583 79.78614583 80.38614583 78.18614583 54.48614583
49 50 51 52 53 54
43.33958333 -15.60708333 -36.38708333 -56.48708333 -77.16677083 -88.56677083
55 56 57 58 59 60
-92.06677083 -90.00677083 -88.68677083 -65.58677083 -64.98677083 -33.68677083
61
-28.73333333
> postscript(file="/var/www/html/rcomp/tmp/6rtn21259397191.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 -18.26718750 NA
1 -12.81385417 -18.26718750
2 -3.39385417 -12.81385417
3 3.30614583 -3.39385417
4 -18.67354167 3.30614583
5 -22.57354167 -18.67354167
6 -26.07354167 -22.57354167
7 -43.61354167 -26.07354167
8 -55.99354167 -43.61354167
9 -45.89354167 -55.99354167
10 -52.19354167 -45.89354167
11 -43.89354167 -52.19354167
12 -35.44010417 -43.89354167
13 -31.88677083 -35.44010417
14 -16.66677083 -31.88677083
15 3.23322917 -16.66677083
16 -9.74645833 3.23322917
17 -9.64645833 -9.74645833
18 -18.44645833 -9.64645833
19 0.01354167 -18.44645833
20 32.43354167 0.01354167
21 10.23354167 32.43354167
22 21.63354167 10.23354167
23 26.53354167 21.63354167
24 35.08697917 26.53354167
25 45.54031250 35.08697917
26 50.86031250 45.54031250
27 58.86031250 50.86031250
28 20.28062500 58.86031250
29 16.88062500 20.28062500
30 22.18062500 16.88062500
31 35.44062500 22.18062500
32 32.46062500 35.44062500
33 20.86062500 32.46062500
34 17.36062500 20.86062500
35 -3.43937500 17.36062500
36 4.01406250 -3.43937500
37 14.76739583 4.01406250
38 5.58739583 14.76739583
39 -8.91260417 5.58739583
40 85.30614583 -8.91260417
41 103.90614583 85.30614583
42 114.40614583 103.90614583
43 98.16614583 114.40614583
44 79.78614583 98.16614583
45 80.38614583 79.78614583
46 78.18614583 80.38614583
47 54.48614583 78.18614583
48 43.33958333 54.48614583
49 -15.60708333 43.33958333
50 -36.38708333 -15.60708333
51 -56.48708333 -36.38708333
52 -77.16677083 -56.48708333
53 -88.56677083 -77.16677083
54 -92.06677083 -88.56677083
55 -90.00677083 -92.06677083
56 -88.68677083 -90.00677083
57 -65.58677083 -88.68677083
58 -64.98677083 -65.58677083
59 -33.68677083 -64.98677083
60 -28.73333333 -33.68677083
61 NA -28.73333333
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -12.81385417 -18.26718750
[2,] -3.39385417 -12.81385417
[3,] 3.30614583 -3.39385417
[4,] -18.67354167 3.30614583
[5,] -22.57354167 -18.67354167
[6,] -26.07354167 -22.57354167
[7,] -43.61354167 -26.07354167
[8,] -55.99354167 -43.61354167
[9,] -45.89354167 -55.99354167
[10,] -52.19354167 -45.89354167
[11,] -43.89354167 -52.19354167
[12,] -35.44010417 -43.89354167
[13,] -31.88677083 -35.44010417
[14,] -16.66677083 -31.88677083
[15,] 3.23322917 -16.66677083
[16,] -9.74645833 3.23322917
[17,] -9.64645833 -9.74645833
[18,] -18.44645833 -9.64645833
[19,] 0.01354167 -18.44645833
[20,] 32.43354167 0.01354167
[21,] 10.23354167 32.43354167
[22,] 21.63354167 10.23354167
[23,] 26.53354167 21.63354167
[24,] 35.08697917 26.53354167
[25,] 45.54031250 35.08697917
[26,] 50.86031250 45.54031250
[27,] 58.86031250 50.86031250
[28,] 20.28062500 58.86031250
[29,] 16.88062500 20.28062500
[30,] 22.18062500 16.88062500
[31,] 35.44062500 22.18062500
[32,] 32.46062500 35.44062500
[33,] 20.86062500 32.46062500
[34,] 17.36062500 20.86062500
[35,] -3.43937500 17.36062500
[36,] 4.01406250 -3.43937500
[37,] 14.76739583 4.01406250
[38,] 5.58739583 14.76739583
[39,] -8.91260417 5.58739583
[40,] 85.30614583 -8.91260417
[41,] 103.90614583 85.30614583
[42,] 114.40614583 103.90614583
[43,] 98.16614583 114.40614583
[44,] 79.78614583 98.16614583
[45,] 80.38614583 79.78614583
[46,] 78.18614583 80.38614583
[47,] 54.48614583 78.18614583
[48,] 43.33958333 54.48614583
[49,] -15.60708333 43.33958333
[50,] -36.38708333 -15.60708333
[51,] -56.48708333 -36.38708333
[52,] -77.16677083 -56.48708333
[53,] -88.56677083 -77.16677083
[54,] -92.06677083 -88.56677083
[55,] -90.00677083 -92.06677083
[56,] -88.68677083 -90.00677083
[57,] -65.58677083 -88.68677083
[58,] -64.98677083 -65.58677083
[59,] -33.68677083 -64.98677083
[60,] -28.73333333 -33.68677083
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -12.81385417 -18.26718750
2 -3.39385417 -12.81385417
3 3.30614583 -3.39385417
4 -18.67354167 3.30614583
5 -22.57354167 -18.67354167
6 -26.07354167 -22.57354167
7 -43.61354167 -26.07354167
8 -55.99354167 -43.61354167
9 -45.89354167 -55.99354167
10 -52.19354167 -45.89354167
11 -43.89354167 -52.19354167
12 -35.44010417 -43.89354167
13 -31.88677083 -35.44010417
14 -16.66677083 -31.88677083
15 3.23322917 -16.66677083
16 -9.74645833 3.23322917
17 -9.64645833 -9.74645833
18 -18.44645833 -9.64645833
19 0.01354167 -18.44645833
20 32.43354167 0.01354167
21 10.23354167 32.43354167
22 21.63354167 10.23354167
23 26.53354167 21.63354167
24 35.08697917 26.53354167
25 45.54031250 35.08697917
26 50.86031250 45.54031250
27 58.86031250 50.86031250
28 20.28062500 58.86031250
29 16.88062500 20.28062500
30 22.18062500 16.88062500
31 35.44062500 22.18062500
32 32.46062500 35.44062500
33 20.86062500 32.46062500
34 17.36062500 20.86062500
35 -3.43937500 17.36062500
36 4.01406250 -3.43937500
37 14.76739583 4.01406250
38 5.58739583 14.76739583
39 -8.91260417 5.58739583
40 85.30614583 -8.91260417
41 103.90614583 85.30614583
42 114.40614583 103.90614583
43 98.16614583 114.40614583
44 79.78614583 98.16614583
45 80.38614583 79.78614583
46 78.18614583 80.38614583
47 54.48614583 78.18614583
48 43.33958333 54.48614583
49 -15.60708333 43.33958333
50 -36.38708333 -15.60708333
51 -56.48708333 -36.38708333
52 -77.16677083 -56.48708333
53 -88.56677083 -77.16677083
54 -92.06677083 -88.56677083
55 -90.00677083 -92.06677083
56 -88.68677083 -90.00677083
57 -65.58677083 -88.68677083
58 -64.98677083 -65.58677083
59 -33.68677083 -64.98677083
60 -28.73333333 -33.68677083
> 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/70lxs1259397191.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/83y6h1259397191.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/9nime1259397191.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/10gq7n1259397191.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/115a821259397191.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/124kql1259397191.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/13m95h1259397191.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/14i12r1259397192.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/15bycc1259397192.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/16vmg41259397192.tab")
+ }
>
> system("convert tmp/1vr4d1259397191.ps tmp/1vr4d1259397191.png")
> system("convert tmp/227z91259397191.ps tmp/227z91259397191.png")
> system("convert tmp/323h71259397191.ps tmp/323h71259397191.png")
> system("convert tmp/4c9r81259397191.ps tmp/4c9r81259397191.png")
> system("convert tmp/5djpk1259397191.ps tmp/5djpk1259397191.png")
> system("convert tmp/6rtn21259397191.ps tmp/6rtn21259397191.png")
> system("convert tmp/70lxs1259397191.ps tmp/70lxs1259397191.png")
> system("convert tmp/83y6h1259397191.ps tmp/83y6h1259397191.png")
> system("convert tmp/9nime1259397191.ps tmp/9nime1259397191.png")
> system("convert tmp/10gq7n1259397191.ps tmp/10gq7n1259397191.png")
>
>
> proc.time()
user system elapsed
2.381 1.533 4.482