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(119.93
+ ,111.4
+ ,101.21
+ ,108.01
+ ,94.76
+ ,87.4
+ ,119.93
+ ,101.21
+ ,95.26
+ ,96.8
+ ,94.76
+ ,119.93
+ ,117.96
+ ,114.1
+ ,95.26
+ ,94.76
+ ,115.86
+ ,110.3
+ ,117.96
+ ,95.26
+ ,111.44
+ ,103.9
+ ,115.86
+ ,117.96
+ ,108.16
+ ,101.6
+ ,111.44
+ ,115.86
+ ,108.77
+ ,94.6
+ ,108.16
+ ,111.44
+ ,109.45
+ ,95.9
+ ,108.77
+ ,108.16
+ ,124.83
+ ,104.7
+ ,109.45
+ ,108.77
+ ,115.31
+ ,102.8
+ ,124.83
+ ,109.45
+ ,109.49
+ ,98.1
+ ,115.31
+ ,124.83
+ ,124.24
+ ,113.9
+ ,109.49
+ ,115.31
+ ,92.85
+ ,80.9
+ ,124.24
+ ,109.49
+ ,98.42
+ ,95.7
+ ,92.85
+ ,124.24
+ ,120.88
+ ,113.2
+ ,98.42
+ ,92.85
+ ,111.72
+ ,105.9
+ ,120.88
+ ,98.42
+ ,116.1
+ ,108.8
+ ,111.72
+ ,120.88
+ ,109.37
+ ,102.3
+ ,116.1
+ ,111.72
+ ,111.65
+ ,99
+ ,109.37
+ ,116.1
+ ,114.29
+ ,100.7
+ ,111.65
+ ,109.37
+ ,133.68
+ ,115.5
+ ,114.29
+ ,111.65
+ ,114.27
+ ,100.7
+ ,133.68
+ ,114.29
+ ,126.49
+ ,109.9
+ ,114.27
+ ,133.68
+ ,131
+ ,114.6
+ ,126.49
+ ,114.27
+ ,104
+ ,85.4
+ ,131
+ ,126.49
+ ,108.88
+ ,100.5
+ ,104
+ ,131
+ ,128.48
+ ,114.8
+ ,108.88
+ ,104
+ ,132.44
+ ,116.5
+ ,128.48
+ ,108.88
+ ,128.04
+ ,112.9
+ ,132.44
+ ,128.48
+ ,116.35
+ ,102
+ ,128.04
+ ,132.44
+ ,120.93
+ ,106
+ ,116.35
+ ,128.04
+ ,118.59
+ ,105.3
+ ,120.93
+ ,116.35
+ ,133.1
+ ,118.8
+ ,118.59
+ ,120.93
+ ,121.05
+ ,106.1
+ ,133.1
+ ,118.59
+ ,127.62
+ ,109.3
+ ,121.05
+ ,133.1
+ ,135.44
+ ,117.2
+ ,127.62
+ ,121.05
+ ,114.88
+ ,92.5
+ ,135.44
+ ,127.62
+ ,114.34
+ ,104.2
+ ,114.88
+ ,135.44
+ ,128.85
+ ,112.5
+ ,114.34
+ ,114.88
+ ,138.9
+ ,122.4
+ ,128.85
+ ,114.34
+ ,129.44
+ ,113.3
+ ,138.9
+ ,128.85
+ ,114.96
+ ,100
+ ,129.44
+ ,138.9
+ ,127.98
+ ,110.7
+ ,114.96
+ ,129.44
+ ,127.03
+ ,112.8
+ ,127.98
+ ,114.96
+ ,128.75
+ ,109.8
+ ,127.03
+ ,127.98
+ ,137.91
+ ,117.3
+ ,128.75
+ ,127.03
+ ,128.37
+ ,109.1
+ ,137.91
+ ,128.75
+ ,135.9
+ ,115.9
+ ,128.37
+ ,137.91
+ ,122.19
+ ,96
+ ,135.9
+ ,128.37
+ ,113.08
+ ,99.8
+ ,122.19
+ ,135.9
+ ,136.2
+ ,116.8
+ ,113.08
+ ,122.19
+ ,138
+ ,115.7
+ ,136.2
+ ,113.08
+ ,115.24
+ ,99.4
+ ,138
+ ,136.2
+ ,110.95
+ ,94.3
+ ,115.24
+ ,138
+ ,99.23
+ ,91
+ ,110.95
+ ,115.24
+ ,102.39
+ ,93.2
+ ,99.23
+ ,110.95
+ ,112.67
+ ,103.1
+ ,102.39
+ ,99.23)
+ ,dim=c(4
+ ,58)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2')
+ ,1:58))
> y <- array(NA,dim=c(4,58),dimnames=list(c('Y','X','Y1','Y2'),1:58))
> 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 Y1 Y2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 119.93 111.4 101.21 108.01 1 0 0 0 0 0 0 0 0 0 0 1
2 94.76 87.4 119.93 101.21 0 1 0 0 0 0 0 0 0 0 0 2
3 95.26 96.8 94.76 119.93 0 0 1 0 0 0 0 0 0 0 0 3
4 117.96 114.1 95.26 94.76 0 0 0 1 0 0 0 0 0 0 0 4
5 115.86 110.3 117.96 95.26 0 0 0 0 1 0 0 0 0 0 0 5
6 111.44 103.9 115.86 117.96 0 0 0 0 0 1 0 0 0 0 0 6
7 108.16 101.6 111.44 115.86 0 0 0 0 0 0 1 0 0 0 0 7
8 108.77 94.6 108.16 111.44 0 0 0 0 0 0 0 1 0 0 0 8
9 109.45 95.9 108.77 108.16 0 0 0 0 0 0 0 0 1 0 0 9
10 124.83 104.7 109.45 108.77 0 0 0 0 0 0 0 0 0 1 0 10
11 115.31 102.8 124.83 109.45 0 0 0 0 0 0 0 0 0 0 1 11
12 109.49 98.1 115.31 124.83 0 0 0 0 0 0 0 0 0 0 0 12
13 124.24 113.9 109.49 115.31 1 0 0 0 0 0 0 0 0 0 0 13
14 92.85 80.9 124.24 109.49 0 1 0 0 0 0 0 0 0 0 0 14
15 98.42 95.7 92.85 124.24 0 0 1 0 0 0 0 0 0 0 0 15
16 120.88 113.2 98.42 92.85 0 0 0 1 0 0 0 0 0 0 0 16
17 111.72 105.9 120.88 98.42 0 0 0 0 1 0 0 0 0 0 0 17
18 116.10 108.8 111.72 120.88 0 0 0 0 0 1 0 0 0 0 0 18
19 109.37 102.3 116.10 111.72 0 0 0 0 0 0 1 0 0 0 0 19
20 111.65 99.0 109.37 116.10 0 0 0 0 0 0 0 1 0 0 0 20
21 114.29 100.7 111.65 109.37 0 0 0 0 0 0 0 0 1 0 0 21
22 133.68 115.5 114.29 111.65 0 0 0 0 0 0 0 0 0 1 0 22
23 114.27 100.7 133.68 114.29 0 0 0 0 0 0 0 0 0 0 1 23
24 126.49 109.9 114.27 133.68 0 0 0 0 0 0 0 0 0 0 0 24
25 131.00 114.6 126.49 114.27 1 0 0 0 0 0 0 0 0 0 0 25
26 104.00 85.4 131.00 126.49 0 1 0 0 0 0 0 0 0 0 0 26
27 108.88 100.5 104.00 131.00 0 0 1 0 0 0 0 0 0 0 0 27
28 128.48 114.8 108.88 104.00 0 0 0 1 0 0 0 0 0 0 0 28
29 132.44 116.5 128.48 108.88 0 0 0 0 1 0 0 0 0 0 0 29
30 128.04 112.9 132.44 128.48 0 0 0 0 0 1 0 0 0 0 0 30
31 116.35 102.0 128.04 132.44 0 0 0 0 0 0 1 0 0 0 0 31
32 120.93 106.0 116.35 128.04 0 0 0 0 0 0 0 1 0 0 0 32
33 118.59 105.3 120.93 116.35 0 0 0 0 0 0 0 0 1 0 0 33
34 133.10 118.8 118.59 120.93 0 0 0 0 0 0 0 0 0 1 0 34
35 121.05 106.1 133.10 118.59 0 0 0 0 0 0 0 0 0 0 1 35
36 127.62 109.3 121.05 133.10 0 0 0 0 0 0 0 0 0 0 0 36
37 135.44 117.2 127.62 121.05 1 0 0 0 0 0 0 0 0 0 0 37
38 114.88 92.5 135.44 127.62 0 1 0 0 0 0 0 0 0 0 0 38
39 114.34 104.2 114.88 135.44 0 0 1 0 0 0 0 0 0 0 0 39
40 128.85 112.5 114.34 114.88 0 0 0 1 0 0 0 0 0 0 0 40
41 138.90 122.4 128.85 114.34 0 0 0 0 1 0 0 0 0 0 0 41
42 129.44 113.3 138.90 128.85 0 0 0 0 0 1 0 0 0 0 0 42
43 114.96 100.0 129.44 138.90 0 0 0 0 0 0 1 0 0 0 0 43
44 127.98 110.7 114.96 129.44 0 0 0 0 0 0 0 1 0 0 0 44
45 127.03 112.8 127.98 114.96 0 0 0 0 0 0 0 0 1 0 0 45
46 128.75 109.8 127.03 127.98 0 0 0 0 0 0 0 0 0 1 0 46
47 137.91 117.3 128.75 127.03 0 0 0 0 0 0 0 0 0 0 1 47
48 128.37 109.1 137.91 128.75 0 0 0 0 0 0 0 0 0 0 0 48
49 135.90 115.9 128.37 137.91 1 0 0 0 0 0 0 0 0 0 0 49
50 122.19 96.0 135.90 128.37 0 1 0 0 0 0 0 0 0 0 0 50
51 113.08 99.8 122.19 135.90 0 0 1 0 0 0 0 0 0 0 0 51
52 136.20 116.8 113.08 122.19 0 0 0 1 0 0 0 0 0 0 0 52
53 138.00 115.7 136.20 113.08 0 0 0 0 1 0 0 0 0 0 0 53
54 115.24 99.4 138.00 136.20 0 0 0 0 0 1 0 0 0 0 0 54
55 110.95 94.3 115.24 138.00 0 0 0 0 0 0 1 0 0 0 0 55
56 99.23 91.0 110.95 115.24 0 0 0 0 0 0 0 1 0 0 0 56
57 102.39 93.2 99.23 110.95 0 0 0 0 0 0 0 0 1 0 0 57
58 112.67 103.1 102.39 99.23 0 0 0 0 0 0 0 0 0 1 0 58
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 M1 M2
-50.96405 1.07016 0.17399 0.29398 1.59008 4.34910
M3 M4 M5 M6 M7 M8
-6.16809 5.25576 2.64206 -3.92510 -2.86913 2.20032
M9 M10 M11 t
3.48719 5.68633 1.38569 0.01280
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.37639 -1.39009 0.05418 1.10854 6.91513
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -50.96405 9.03349 -5.642 1.3e-06 ***
X 1.07016 0.07956 13.450 < 2e-16 ***
Y1 0.17399 0.06156 2.826 0.007177 **
Y2 0.29398 0.06887 4.269 0.000110 ***
M1 1.59008 2.05716 0.773 0.443882
M2 4.34910 2.46618 1.763 0.085091 .
M3 -6.16809 2.03579 -3.030 0.004176 **
M4 5.25576 2.61804 2.008 0.051158 .
M5 2.64206 2.65565 0.995 0.325491
M6 -3.92510 1.84635 -2.126 0.039436 *
M7 -2.86913 1.84632 -1.554 0.127696
M8 2.20032 1.92968 1.140 0.260642
M9 3.48719 2.12188 1.643 0.107759
M10 5.68633 2.20335 2.581 0.013442 *
M11 1.38569 2.19842 0.630 0.531904
t 0.01280 0.03206 0.399 0.691717
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.659 on 42 degrees of freedom
Multiple R-squared: 0.9621, Adjusted R-squared: 0.9485
F-statistic: 71.02 on 15 and 42 DF, p-value: < 2.2e-16
> 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.4517427 0.9034854 0.5482573
[2,] 0.3871568 0.7743136 0.6128432
[3,] 0.3171503 0.6343007 0.6828497
[4,] 0.3768087 0.7536174 0.6231913
[5,] 0.2851610 0.5703219 0.7148390
[6,] 0.5962869 0.8074262 0.4037131
[7,] 0.7628766 0.4742469 0.2371234
[8,] 0.7007652 0.5984696 0.2992348
[9,] 0.6907704 0.6184593 0.3092296
[10,] 0.6205513 0.7588974 0.3794487
[11,] 0.6574664 0.6850671 0.3425336
[12,] 0.5695940 0.8608120 0.4304060
[13,] 0.5052275 0.9895449 0.4947725
[14,] 0.5965537 0.8068926 0.4034463
[15,] 0.6812065 0.6375869 0.3187935
[16,] 0.8194171 0.3611657 0.1805829
[17,] 0.7197788 0.5604424 0.2802212
[18,] 0.6386672 0.7226657 0.3613328
[19,] 0.6945001 0.6109998 0.3054999
[20,] 0.5785467 0.8429067 0.4214533
[21,] 0.4033528 0.8067057 0.5966472
> postscript(file="/var/www/html/rcomp/tmp/1fq5f1258661756.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/269we1258661756.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/380yp1258661756.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/4g7f51258661756.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/5gryw1258661756.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 = 58
Frequency = 1
1 2 3 4 5 6
0.71358989 -2.80233202 -2.98150396 -2.91951028 -2.44846825 0.22699429
7 8 9 10 11 12
-0.27403664 4.61489633 3.46213056 6.91513366 0.84046944 -1.44191697
13 14 15 16 17 18
-1.39207748 -1.09392375 0.26731703 0.82171206 -3.47039230 -0.64853436
19 20 21 22 23 24
0.43951952 1.05210324 2.15493989 2.36502123 -1.06844335 0.35579877
25 26 27 28 29 30
1.81316353 -1.08703886 1.50968786 1.45808659 1.35497576 0.91096430
31 32 33 34 35 36
-0.58166539 -2.03716440 -2.28797878 -5.37638761 -1.38413112 0.96515396
37 38 39 40 41 42
1.12734662 0.93650026 -0.34175961 -0.01260764 -0.32209319 0.49655496
43 44 45 46 47 48
-2.12764430 -0.34027553 -2.84578311 -3.78954505 1.61210503 0.12096424
49 50 51 52 53 54
-2.26202255 4.04679437 1.54625868 0.65231927 4.88597798 -0.98597919
55 56 57 58
2.54382681 -3.28955963 -0.48330857 -0.11422223
> postscript(file="/var/www/html/rcomp/tmp/63c431258661756.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 = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 0.71358989 NA
1 -2.80233202 0.71358989
2 -2.98150396 -2.80233202
3 -2.91951028 -2.98150396
4 -2.44846825 -2.91951028
5 0.22699429 -2.44846825
6 -0.27403664 0.22699429
7 4.61489633 -0.27403664
8 3.46213056 4.61489633
9 6.91513366 3.46213056
10 0.84046944 6.91513366
11 -1.44191697 0.84046944
12 -1.39207748 -1.44191697
13 -1.09392375 -1.39207748
14 0.26731703 -1.09392375
15 0.82171206 0.26731703
16 -3.47039230 0.82171206
17 -0.64853436 -3.47039230
18 0.43951952 -0.64853436
19 1.05210324 0.43951952
20 2.15493989 1.05210324
21 2.36502123 2.15493989
22 -1.06844335 2.36502123
23 0.35579877 -1.06844335
24 1.81316353 0.35579877
25 -1.08703886 1.81316353
26 1.50968786 -1.08703886
27 1.45808659 1.50968786
28 1.35497576 1.45808659
29 0.91096430 1.35497576
30 -0.58166539 0.91096430
31 -2.03716440 -0.58166539
32 -2.28797878 -2.03716440
33 -5.37638761 -2.28797878
34 -1.38413112 -5.37638761
35 0.96515396 -1.38413112
36 1.12734662 0.96515396
37 0.93650026 1.12734662
38 -0.34175961 0.93650026
39 -0.01260764 -0.34175961
40 -0.32209319 -0.01260764
41 0.49655496 -0.32209319
42 -2.12764430 0.49655496
43 -0.34027553 -2.12764430
44 -2.84578311 -0.34027553
45 -3.78954505 -2.84578311
46 1.61210503 -3.78954505
47 0.12096424 1.61210503
48 -2.26202255 0.12096424
49 4.04679437 -2.26202255
50 1.54625868 4.04679437
51 0.65231927 1.54625868
52 4.88597798 0.65231927
53 -0.98597919 4.88597798
54 2.54382681 -0.98597919
55 -3.28955963 2.54382681
56 -0.48330857 -3.28955963
57 -0.11422223 -0.48330857
58 NA -0.11422223
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.80233202 0.71358989
[2,] -2.98150396 -2.80233202
[3,] -2.91951028 -2.98150396
[4,] -2.44846825 -2.91951028
[5,] 0.22699429 -2.44846825
[6,] -0.27403664 0.22699429
[7,] 4.61489633 -0.27403664
[8,] 3.46213056 4.61489633
[9,] 6.91513366 3.46213056
[10,] 0.84046944 6.91513366
[11,] -1.44191697 0.84046944
[12,] -1.39207748 -1.44191697
[13,] -1.09392375 -1.39207748
[14,] 0.26731703 -1.09392375
[15,] 0.82171206 0.26731703
[16,] -3.47039230 0.82171206
[17,] -0.64853436 -3.47039230
[18,] 0.43951952 -0.64853436
[19,] 1.05210324 0.43951952
[20,] 2.15493989 1.05210324
[21,] 2.36502123 2.15493989
[22,] -1.06844335 2.36502123
[23,] 0.35579877 -1.06844335
[24,] 1.81316353 0.35579877
[25,] -1.08703886 1.81316353
[26,] 1.50968786 -1.08703886
[27,] 1.45808659 1.50968786
[28,] 1.35497576 1.45808659
[29,] 0.91096430 1.35497576
[30,] -0.58166539 0.91096430
[31,] -2.03716440 -0.58166539
[32,] -2.28797878 -2.03716440
[33,] -5.37638761 -2.28797878
[34,] -1.38413112 -5.37638761
[35,] 0.96515396 -1.38413112
[36,] 1.12734662 0.96515396
[37,] 0.93650026 1.12734662
[38,] -0.34175961 0.93650026
[39,] -0.01260764 -0.34175961
[40,] -0.32209319 -0.01260764
[41,] 0.49655496 -0.32209319
[42,] -2.12764430 0.49655496
[43,] -0.34027553 -2.12764430
[44,] -2.84578311 -0.34027553
[45,] -3.78954505 -2.84578311
[46,] 1.61210503 -3.78954505
[47,] 0.12096424 1.61210503
[48,] -2.26202255 0.12096424
[49,] 4.04679437 -2.26202255
[50,] 1.54625868 4.04679437
[51,] 0.65231927 1.54625868
[52,] 4.88597798 0.65231927
[53,] -0.98597919 4.88597798
[54,] 2.54382681 -0.98597919
[55,] -3.28955963 2.54382681
[56,] -0.48330857 -3.28955963
[57,] -0.11422223 -0.48330857
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.80233202 0.71358989
2 -2.98150396 -2.80233202
3 -2.91951028 -2.98150396
4 -2.44846825 -2.91951028
5 0.22699429 -2.44846825
6 -0.27403664 0.22699429
7 4.61489633 -0.27403664
8 3.46213056 4.61489633
9 6.91513366 3.46213056
10 0.84046944 6.91513366
11 -1.44191697 0.84046944
12 -1.39207748 -1.44191697
13 -1.09392375 -1.39207748
14 0.26731703 -1.09392375
15 0.82171206 0.26731703
16 -3.47039230 0.82171206
17 -0.64853436 -3.47039230
18 0.43951952 -0.64853436
19 1.05210324 0.43951952
20 2.15493989 1.05210324
21 2.36502123 2.15493989
22 -1.06844335 2.36502123
23 0.35579877 -1.06844335
24 1.81316353 0.35579877
25 -1.08703886 1.81316353
26 1.50968786 -1.08703886
27 1.45808659 1.50968786
28 1.35497576 1.45808659
29 0.91096430 1.35497576
30 -0.58166539 0.91096430
31 -2.03716440 -0.58166539
32 -2.28797878 -2.03716440
33 -5.37638761 -2.28797878
34 -1.38413112 -5.37638761
35 0.96515396 -1.38413112
36 1.12734662 0.96515396
37 0.93650026 1.12734662
38 -0.34175961 0.93650026
39 -0.01260764 -0.34175961
40 -0.32209319 -0.01260764
41 0.49655496 -0.32209319
42 -2.12764430 0.49655496
43 -0.34027553 -2.12764430
44 -2.84578311 -0.34027553
45 -3.78954505 -2.84578311
46 1.61210503 -3.78954505
47 0.12096424 1.61210503
48 -2.26202255 0.12096424
49 4.04679437 -2.26202255
50 1.54625868 4.04679437
51 0.65231927 1.54625868
52 4.88597798 0.65231927
53 -0.98597919 4.88597798
54 2.54382681 -0.98597919
55 -3.28955963 2.54382681
56 -0.48330857 -3.28955963
57 -0.11422223 -0.48330857
> 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/7y9nr1258661756.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/8zmlu1258661756.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/9bto71258661756.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/10rjvy1258661756.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/11abrr1258661756.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/12imzw1258661756.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/139x881258661756.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/14w86t1258661756.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/15g19g1258661757.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/16pkv21258661757.tab")
+ }
>
> system("convert tmp/1fq5f1258661756.ps tmp/1fq5f1258661756.png")
> system("convert tmp/269we1258661756.ps tmp/269we1258661756.png")
> system("convert tmp/380yp1258661756.ps tmp/380yp1258661756.png")
> system("convert tmp/4g7f51258661756.ps tmp/4g7f51258661756.png")
> system("convert tmp/5gryw1258661756.ps tmp/5gryw1258661756.png")
> system("convert tmp/63c431258661756.ps tmp/63c431258661756.png")
> system("convert tmp/7y9nr1258661756.ps tmp/7y9nr1258661756.png")
> system("convert tmp/8zmlu1258661756.ps tmp/8zmlu1258661756.png")
> system("convert tmp/9bto71258661756.ps tmp/9bto71258661756.png")
> system("convert tmp/10rjvy1258661756.ps tmp/10rjvy1258661756.png")
>
>
> proc.time()
user system elapsed
2.322 1.540 2.793