R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(15044.5,1,14944.2,1,16754.8,1,14254,1,15454.9,1,15644.8,1,14568.3,1,12520.2,1,14803,1,15873.2,1,14755.3,1,12875.1,1,14291.1,1,14205.3,1,15859.4,1,15258.9,1,15498.6,1,15106.5,1,15023.6,1,12083,1,15761.3,1,16943,1,15070.3,1,13659.6,1,14768.9,0,14725.1,0,15998.1,0,15370.6,0,14956.9,0,15469.7,0,15101.8,0,11703.7,0,16283.6,0,16726.5,0,14968.9,0,14861,0,14583.3,0,15305.8,0,17903.9,0,16379.4,0,15420.3,0,17870.5,0,15912.8,0,13866.5,0,17823.2,0,17872,0,17420.4,0,16704.4,0,15991.2,0,16583.6,0,19123.5,0,17838.7,0,17209.4,0,18586.5,0,16258.1,0,15141.6,0,19202.1,0,17746.5,0,19090.1,0,18040.3,0,17515.5,0,17751.8,0,21072.4,0,17170,0,19439.5,0,19795.4,0,17574.9,0,16165.4,0,19464.6,0,19932.1,0,19961.2,0,17343.4,0,18924.2,0,18574.1,0,21350.6,0,18594.6,0,19823.1,0,20844.4,0,19640.2,0,17735.4,0,19813.6,0,22160,0,20664.3,0,17877.4,0,21211.2,0,21423.1,0,21688.7,0,23243.2,0,21490.2,0,22925.8,0,23184.8,0,18562.2,0),dim=c(2,92),dimnames=list(c('Y','X'),1:92))
> y <- array(NA,dim=c(2,92),dimnames=list(c('Y','X'),1:92))
> 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 15044.5 1 1 0 0 0 0 0 0 0 0 0 0 1
2 14944.2 1 0 1 0 0 0 0 0 0 0 0 0 2
3 16754.8 1 0 0 1 0 0 0 0 0 0 0 0 3
4 14254.0 1 0 0 0 1 0 0 0 0 0 0 0 4
5 15454.9 1 0 0 0 0 1 0 0 0 0 0 0 5
6 15644.8 1 0 0 0 0 0 1 0 0 0 0 0 6
7 14568.3 1 0 0 0 0 0 0 1 0 0 0 0 7
8 12520.2 1 0 0 0 0 0 0 0 1 0 0 0 8
9 14803.0 1 0 0 0 0 0 0 0 0 1 0 0 9
10 15873.2 1 0 0 0 0 0 0 0 0 0 1 0 10
11 14755.3 1 0 0 0 0 0 0 0 0 0 0 1 11
12 12875.1 1 0 0 0 0 0 0 0 0 0 0 0 12
13 14291.1 1 1 0 0 0 0 0 0 0 0 0 0 13
14 14205.3 1 0 1 0 0 0 0 0 0 0 0 0 14
15 15859.4 1 0 0 1 0 0 0 0 0 0 0 0 15
16 15258.9 1 0 0 0 1 0 0 0 0 0 0 0 16
17 15498.6 1 0 0 0 0 1 0 0 0 0 0 0 17
18 15106.5 1 0 0 0 0 0 1 0 0 0 0 0 18
19 15023.6 1 0 0 0 0 0 0 1 0 0 0 0 19
20 12083.0 1 0 0 0 0 0 0 0 1 0 0 0 20
21 15761.3 1 0 0 0 0 0 0 0 0 1 0 0 21
22 16943.0 1 0 0 0 0 0 0 0 0 0 1 0 22
23 15070.3 1 0 0 0 0 0 0 0 0 0 0 1 23
24 13659.6 1 0 0 0 0 0 0 0 0 0 0 0 24
25 14768.9 0 1 0 0 0 0 0 0 0 0 0 0 25
26 14725.1 0 0 1 0 0 0 0 0 0 0 0 0 26
27 15998.1 0 0 0 1 0 0 0 0 0 0 0 0 27
28 15370.6 0 0 0 0 1 0 0 0 0 0 0 0 28
29 14956.9 0 0 0 0 0 1 0 0 0 0 0 0 29
30 15469.7 0 0 0 0 0 0 1 0 0 0 0 0 30
31 15101.8 0 0 0 0 0 0 0 1 0 0 0 0 31
32 11703.7 0 0 0 0 0 0 0 0 1 0 0 0 32
33 16283.6 0 0 0 0 0 0 0 0 0 1 0 0 33
34 16726.5 0 0 0 0 0 0 0 0 0 0 1 0 34
35 14968.9 0 0 0 0 0 0 0 0 0 0 0 1 35
36 14861.0 0 0 0 0 0 0 0 0 0 0 0 0 36
37 14583.3 0 1 0 0 0 0 0 0 0 0 0 0 37
38 15305.8 0 0 1 0 0 0 0 0 0 0 0 0 38
39 17903.9 0 0 0 1 0 0 0 0 0 0 0 0 39
40 16379.4 0 0 0 0 1 0 0 0 0 0 0 0 40
41 15420.3 0 0 0 0 0 1 0 0 0 0 0 0 41
42 17870.5 0 0 0 0 0 0 1 0 0 0 0 0 42
43 15912.8 0 0 0 0 0 0 0 1 0 0 0 0 43
44 13866.5 0 0 0 0 0 0 0 0 1 0 0 0 44
45 17823.2 0 0 0 0 0 0 0 0 0 1 0 0 45
46 17872.0 0 0 0 0 0 0 0 0 0 0 1 0 46
47 17420.4 0 0 0 0 0 0 0 0 0 0 0 1 47
48 16704.4 0 0 0 0 0 0 0 0 0 0 0 0 48
49 15991.2 0 1 0 0 0 0 0 0 0 0 0 0 49
50 16583.6 0 0 1 0 0 0 0 0 0 0 0 0 50
51 19123.5 0 0 0 1 0 0 0 0 0 0 0 0 51
52 17838.7 0 0 0 0 1 0 0 0 0 0 0 0 52
53 17209.4 0 0 0 0 0 1 0 0 0 0 0 0 53
54 18586.5 0 0 0 0 0 0 1 0 0 0 0 0 54
55 16258.1 0 0 0 0 0 0 0 1 0 0 0 0 55
56 15141.6 0 0 0 0 0 0 0 0 1 0 0 0 56
57 19202.1 0 0 0 0 0 0 0 0 0 1 0 0 57
58 17746.5 0 0 0 0 0 0 0 0 0 0 1 0 58
59 19090.1 0 0 0 0 0 0 0 0 0 0 0 1 59
60 18040.3 0 0 0 0 0 0 0 0 0 0 0 0 60
61 17515.5 0 1 0 0 0 0 0 0 0 0 0 0 61
62 17751.8 0 0 1 0 0 0 0 0 0 0 0 0 62
63 21072.4 0 0 0 1 0 0 0 0 0 0 0 0 63
64 17170.0 0 0 0 0 1 0 0 0 0 0 0 0 64
65 19439.5 0 0 0 0 0 1 0 0 0 0 0 0 65
66 19795.4 0 0 0 0 0 0 1 0 0 0 0 0 66
67 17574.9 0 0 0 0 0 0 0 1 0 0 0 0 67
68 16165.4 0 0 0 0 0 0 0 0 1 0 0 0 68
69 19464.6 0 0 0 0 0 0 0 0 0 1 0 0 69
70 19932.1 0 0 0 0 0 0 0 0 0 0 1 0 70
71 19961.2 0 0 0 0 0 0 0 0 0 0 0 1 71
72 17343.4 0 0 0 0 0 0 0 0 0 0 0 0 72
73 18924.2 0 1 0 0 0 0 0 0 0 0 0 0 73
74 18574.1 0 0 1 0 0 0 0 0 0 0 0 0 74
75 21350.6 0 0 0 1 0 0 0 0 0 0 0 0 75
76 18594.6 0 0 0 0 1 0 0 0 0 0 0 0 76
77 19823.1 0 0 0 0 0 1 0 0 0 0 0 0 77
78 20844.4 0 0 0 0 0 0 1 0 0 0 0 0 78
79 19640.2 0 0 0 0 0 0 0 1 0 0 0 0 79
80 17735.4 0 0 0 0 0 0 0 0 1 0 0 0 80
81 19813.6 0 0 0 0 0 0 0 0 0 1 0 0 81
82 22160.0 0 0 0 0 0 0 0 0 0 0 1 0 82
83 20664.3 0 0 0 0 0 0 0 0 0 0 0 1 83
84 17877.4 0 0 0 0 0 0 0 0 0 0 0 0 84
85 21211.2 0 1 0 0 0 0 0 0 0 0 0 0 85
86 21423.1 0 0 1 0 0 0 0 0 0 0 0 0 86
87 21688.7 0 0 0 1 0 0 0 0 0 0 0 0 87
88 23243.2 0 0 0 0 1 0 0 0 0 0 0 0 88
89 21490.2 0 0 0 0 0 1 0 0 0 0 0 0 89
90 22925.8 0 0 0 0 0 0 1 0 0 0 0 0 90
91 23184.8 0 0 0 0 0 0 0 1 0 0 0 0 91
92 18562.2 0 0 0 0 0 0 0 0 1 0 0 0 92
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
10563.8 1566.0 1198.6 1244.4 3172.2 1614.9
M5 M6 M7 M8 M9 M10
1660.8 2427.6 1203.2 -1334.6 1990.4 2474.4
M11 t
1611.9 102.0
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1538.74 -494.96 -14.25 460.59 2132.91
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10563.785 439.524 24.035 < 2e-16 ***
X 1566.030 302.055 5.185 1.66e-06 ***
M1 1198.583 426.410 2.811 0.006245 **
M2 1244.438 426.049 2.921 0.004563 **
M3 3172.207 425.746 7.451 1.08e-10 ***
M4 1614.925 425.502 3.795 0.000290 ***
M5 1660.831 425.317 3.905 0.000199 ***
M6 2427.637 425.191 5.710 1.96e-07 ***
M7 1203.218 425.123 2.830 0.005912 **
M8 -1334.627 425.115 -3.139 0.002390 **
M9 1990.409 439.250 4.531 2.08e-05 ***
M10 2474.363 439.108 5.635 2.67e-07 ***
M11 1611.932 439.022 3.672 0.000439 ***
t 102.032 5.005 20.385 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 821.3 on 78 degrees of freedom
Multiple R-squared: 0.9137, Adjusted R-squared: 0.8993
F-statistic: 63.49 on 13 and 78 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.438955227 0.877910454 0.5610448
[2,] 0.284355519 0.568711037 0.7156445
[3,] 0.223503272 0.447006545 0.7764967
[4,] 0.133335399 0.266670799 0.8666646
[5,] 0.159987709 0.319975419 0.8400123
[6,] 0.186867990 0.373735980 0.8131320
[7,] 0.123268961 0.246537922 0.8767310
[8,] 0.096945087 0.193890175 0.9030549
[9,] 0.062793032 0.125586064 0.9372070
[10,] 0.038046202 0.076092404 0.9619538
[11,] 0.024174118 0.048348237 0.9758259
[12,] 0.018167689 0.036335377 0.9818323
[13,] 0.013664226 0.027328453 0.9863358
[14,] 0.007798432 0.015596865 0.9922016
[15,] 0.004498866 0.008997732 0.9955011
[16,] 0.003707004 0.007414009 0.9962930
[17,] 0.004847673 0.009695347 0.9951523
[18,] 0.002789162 0.005578324 0.9972108
[19,] 0.001757328 0.003514656 0.9982427
[20,] 0.005617462 0.011234925 0.9943825
[21,] 0.003692591 0.007385182 0.9963074
[22,] 0.002389993 0.004779987 0.9976100
[23,] 0.005620857 0.011241714 0.9943791
[24,] 0.006017553 0.012035105 0.9939824
[25,] 0.004406841 0.008813683 0.9955932
[26,] 0.018995008 0.037990015 0.9810050
[27,] 0.012595361 0.025190722 0.9874046
[28,] 0.012774741 0.025549483 0.9872253
[29,] 0.020309814 0.040619628 0.9796902
[30,] 0.014946834 0.029893669 0.9850532
[31,] 0.022282098 0.044564197 0.9777179
[32,] 0.060275155 0.120550311 0.9397248
[33,] 0.043642092 0.087284183 0.9563579
[34,] 0.030009040 0.060018079 0.9699910
[35,] 0.026012302 0.052024604 0.9739877
[36,] 0.024646249 0.049292497 0.9753538
[37,] 0.016111824 0.032223647 0.9838882
[38,] 0.012284455 0.024568911 0.9877155
[39,] 0.011103688 0.022207376 0.9888963
[40,] 0.008583965 0.017167931 0.9914160
[41,] 0.013995523 0.027991047 0.9860045
[42,] 0.014405702 0.028811404 0.9855943
[43,] 0.020394508 0.040789017 0.9796055
[44,] 0.091195116 0.182390231 0.9088049
[45,] 0.063486406 0.126972812 0.9365136
[46,] 0.042747228 0.085494456 0.9572528
[47,] 0.096745716 0.193491432 0.9032543
[48,] 0.113569205 0.227138411 0.8864308
[49,] 0.138484379 0.276968758 0.8615156
[50,] 0.110090764 0.220181528 0.8899092
[51,] 0.105484555 0.210969109 0.8945154
[52,] 0.080355597 0.160711194 0.9196444
[53,] 0.086705724 0.173411449 0.9132943
[54,] 0.054790712 0.109581424 0.9452093
[55,] 0.055000696 0.110001392 0.9449993
[56,] 0.068497452 0.136994905 0.9315025
[57,] 0.037445622 0.074891244 0.9625544
[58,] 0.020815401 0.041630801 0.9791846
[59,] 0.039430380 0.078860760 0.9605696
> postscript(file="/var/www/html/rcomp/tmp/19mn51229016348.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/2y04h1229016348.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/3gty81229016348.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/4w5s11229016348.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/59efa1229016348.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 = 92
Frequency = 1
1 2 3 4 5 6
1614.070441 1365.882941 1146.682941 101.132941 1154.095441 475.157941
7 8 9 10 11 12
521.045441 908.757941 -235.509606 248.704680 -108.795320 -479.095320
13 14 15 16 17 18
-363.709636 -597.397136 -973.097136 -118.347136 -26.584636 -1287.522136
19 20 21 22 23 24
-248.034636 -752.822136 -501.589683 94.124603 -1018.175397 -918.975397
25 26 27 28 29 30
455.740058 264.052558 -492.747442 335.002558 -226.634942 -582.672442
31 32 33 34 35 36
171.815058 -790.472442 362.360011 219.274297 -777.925703 624.074297
37 38 39 40 41 42
-954.240019 -379.627519 188.672481 119.422481 -987.615019 593.747481
43 44 45 46 47 48
-241.565019 147.947481 677.579934 140.394220 449.194220 1243.094220
49 50 51 52 53 54
-770.720096 -326.207596 183.892404 354.342404 -422.895096 85.367404
55 56 57 58 59 60
-1120.645096 198.667404 832.099858 -1209.485857 894.514143 1354.614143
61 62 63 64 65 66
-470.800173 -382.387673 908.412327 -1538.737673 582.824827 69.887327
67 68 69 70 71 72
-1028.225173 -1.912673 -129.780219 -248.265933 541.234067 -566.665933
73 74 75 76 77 78
-286.480249 -784.467749 -37.767749 -1338.517749 -257.955249 -105.492749
79 80 81 82 83 84
-187.305249 343.707251 -1005.160296 755.253990 19.953990 -1257.046010
85 86 87 88 89 90
776.139674 840.152174 -924.047826 2085.702174 184.764674 751.527174
91 92
2132.914674 -53.872826
> postscript(file="/var/www/html/rcomp/tmp/6gyci1229016348.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 = 92
Frequency = 1
lag(myerror, k = 1) myerror
0 1614.070441 NA
1 1365.882941 1614.070441
2 1146.682941 1365.882941
3 101.132941 1146.682941
4 1154.095441 101.132941
5 475.157941 1154.095441
6 521.045441 475.157941
7 908.757941 521.045441
8 -235.509606 908.757941
9 248.704680 -235.509606
10 -108.795320 248.704680
11 -479.095320 -108.795320
12 -363.709636 -479.095320
13 -597.397136 -363.709636
14 -973.097136 -597.397136
15 -118.347136 -973.097136
16 -26.584636 -118.347136
17 -1287.522136 -26.584636
18 -248.034636 -1287.522136
19 -752.822136 -248.034636
20 -501.589683 -752.822136
21 94.124603 -501.589683
22 -1018.175397 94.124603
23 -918.975397 -1018.175397
24 455.740058 -918.975397
25 264.052558 455.740058
26 -492.747442 264.052558
27 335.002558 -492.747442
28 -226.634942 335.002558
29 -582.672442 -226.634942
30 171.815058 -582.672442
31 -790.472442 171.815058
32 362.360011 -790.472442
33 219.274297 362.360011
34 -777.925703 219.274297
35 624.074297 -777.925703
36 -954.240019 624.074297
37 -379.627519 -954.240019
38 188.672481 -379.627519
39 119.422481 188.672481
40 -987.615019 119.422481
41 593.747481 -987.615019
42 -241.565019 593.747481
43 147.947481 -241.565019
44 677.579934 147.947481
45 140.394220 677.579934
46 449.194220 140.394220
47 1243.094220 449.194220
48 -770.720096 1243.094220
49 -326.207596 -770.720096
50 183.892404 -326.207596
51 354.342404 183.892404
52 -422.895096 354.342404
53 85.367404 -422.895096
54 -1120.645096 85.367404
55 198.667404 -1120.645096
56 832.099858 198.667404
57 -1209.485857 832.099858
58 894.514143 -1209.485857
59 1354.614143 894.514143
60 -470.800173 1354.614143
61 -382.387673 -470.800173
62 908.412327 -382.387673
63 -1538.737673 908.412327
64 582.824827 -1538.737673
65 69.887327 582.824827
66 -1028.225173 69.887327
67 -1.912673 -1028.225173
68 -129.780219 -1.912673
69 -248.265933 -129.780219
70 541.234067 -248.265933
71 -566.665933 541.234067
72 -286.480249 -566.665933
73 -784.467749 -286.480249
74 -37.767749 -784.467749
75 -1338.517749 -37.767749
76 -257.955249 -1338.517749
77 -105.492749 -257.955249
78 -187.305249 -105.492749
79 343.707251 -187.305249
80 -1005.160296 343.707251
81 755.253990 -1005.160296
82 19.953990 755.253990
83 -1257.046010 19.953990
84 776.139674 -1257.046010
85 840.152174 776.139674
86 -924.047826 840.152174
87 2085.702174 -924.047826
88 184.764674 2085.702174
89 751.527174 184.764674
90 2132.914674 751.527174
91 -53.872826 2132.914674
92 NA -53.872826
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1365.882941 1614.070441
[2,] 1146.682941 1365.882941
[3,] 101.132941 1146.682941
[4,] 1154.095441 101.132941
[5,] 475.157941 1154.095441
[6,] 521.045441 475.157941
[7,] 908.757941 521.045441
[8,] -235.509606 908.757941
[9,] 248.704680 -235.509606
[10,] -108.795320 248.704680
[11,] -479.095320 -108.795320
[12,] -363.709636 -479.095320
[13,] -597.397136 -363.709636
[14,] -973.097136 -597.397136
[15,] -118.347136 -973.097136
[16,] -26.584636 -118.347136
[17,] -1287.522136 -26.584636
[18,] -248.034636 -1287.522136
[19,] -752.822136 -248.034636
[20,] -501.589683 -752.822136
[21,] 94.124603 -501.589683
[22,] -1018.175397 94.124603
[23,] -918.975397 -1018.175397
[24,] 455.740058 -918.975397
[25,] 264.052558 455.740058
[26,] -492.747442 264.052558
[27,] 335.002558 -492.747442
[28,] -226.634942 335.002558
[29,] -582.672442 -226.634942
[30,] 171.815058 -582.672442
[31,] -790.472442 171.815058
[32,] 362.360011 -790.472442
[33,] 219.274297 362.360011
[34,] -777.925703 219.274297
[35,] 624.074297 -777.925703
[36,] -954.240019 624.074297
[37,] -379.627519 -954.240019
[38,] 188.672481 -379.627519
[39,] 119.422481 188.672481
[40,] -987.615019 119.422481
[41,] 593.747481 -987.615019
[42,] -241.565019 593.747481
[43,] 147.947481 -241.565019
[44,] 677.579934 147.947481
[45,] 140.394220 677.579934
[46,] 449.194220 140.394220
[47,] 1243.094220 449.194220
[48,] -770.720096 1243.094220
[49,] -326.207596 -770.720096
[50,] 183.892404 -326.207596
[51,] 354.342404 183.892404
[52,] -422.895096 354.342404
[53,] 85.367404 -422.895096
[54,] -1120.645096 85.367404
[55,] 198.667404 -1120.645096
[56,] 832.099858 198.667404
[57,] -1209.485857 832.099858
[58,] 894.514143 -1209.485857
[59,] 1354.614143 894.514143
[60,] -470.800173 1354.614143
[61,] -382.387673 -470.800173
[62,] 908.412327 -382.387673
[63,] -1538.737673 908.412327
[64,] 582.824827 -1538.737673
[65,] 69.887327 582.824827
[66,] -1028.225173 69.887327
[67,] -1.912673 -1028.225173
[68,] -129.780219 -1.912673
[69,] -248.265933 -129.780219
[70,] 541.234067 -248.265933
[71,] -566.665933 541.234067
[72,] -286.480249 -566.665933
[73,] -784.467749 -286.480249
[74,] -37.767749 -784.467749
[75,] -1338.517749 -37.767749
[76,] -257.955249 -1338.517749
[77,] -105.492749 -257.955249
[78,] -187.305249 -105.492749
[79,] 343.707251 -187.305249
[80,] -1005.160296 343.707251
[81,] 755.253990 -1005.160296
[82,] 19.953990 755.253990
[83,] -1257.046010 19.953990
[84,] 776.139674 -1257.046010
[85,] 840.152174 776.139674
[86,] -924.047826 840.152174
[87,] 2085.702174 -924.047826
[88,] 184.764674 2085.702174
[89,] 751.527174 184.764674
[90,] 2132.914674 751.527174
[91,] -53.872826 2132.914674
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1365.882941 1614.070441
2 1146.682941 1365.882941
3 101.132941 1146.682941
4 1154.095441 101.132941
5 475.157941 1154.095441
6 521.045441 475.157941
7 908.757941 521.045441
8 -235.509606 908.757941
9 248.704680 -235.509606
10 -108.795320 248.704680
11 -479.095320 -108.795320
12 -363.709636 -479.095320
13 -597.397136 -363.709636
14 -973.097136 -597.397136
15 -118.347136 -973.097136
16 -26.584636 -118.347136
17 -1287.522136 -26.584636
18 -248.034636 -1287.522136
19 -752.822136 -248.034636
20 -501.589683 -752.822136
21 94.124603 -501.589683
22 -1018.175397 94.124603
23 -918.975397 -1018.175397
24 455.740058 -918.975397
25 264.052558 455.740058
26 -492.747442 264.052558
27 335.002558 -492.747442
28 -226.634942 335.002558
29 -582.672442 -226.634942
30 171.815058 -582.672442
31 -790.472442 171.815058
32 362.360011 -790.472442
33 219.274297 362.360011
34 -777.925703 219.274297
35 624.074297 -777.925703
36 -954.240019 624.074297
37 -379.627519 -954.240019
38 188.672481 -379.627519
39 119.422481 188.672481
40 -987.615019 119.422481
41 593.747481 -987.615019
42 -241.565019 593.747481
43 147.947481 -241.565019
44 677.579934 147.947481
45 140.394220 677.579934
46 449.194220 140.394220
47 1243.094220 449.194220
48 -770.720096 1243.094220
49 -326.207596 -770.720096
50 183.892404 -326.207596
51 354.342404 183.892404
52 -422.895096 354.342404
53 85.367404 -422.895096
54 -1120.645096 85.367404
55 198.667404 -1120.645096
56 832.099858 198.667404
57 -1209.485857 832.099858
58 894.514143 -1209.485857
59 1354.614143 894.514143
60 -470.800173 1354.614143
61 -382.387673 -470.800173
62 908.412327 -382.387673
63 -1538.737673 908.412327
64 582.824827 -1538.737673
65 69.887327 582.824827
66 -1028.225173 69.887327
67 -1.912673 -1028.225173
68 -129.780219 -1.912673
69 -248.265933 -129.780219
70 541.234067 -248.265933
71 -566.665933 541.234067
72 -286.480249 -566.665933
73 -784.467749 -286.480249
74 -37.767749 -784.467749
75 -1338.517749 -37.767749
76 -257.955249 -1338.517749
77 -105.492749 -257.955249
78 -187.305249 -105.492749
79 343.707251 -187.305249
80 -1005.160296 343.707251
81 755.253990 -1005.160296
82 19.953990 755.253990
83 -1257.046010 19.953990
84 776.139674 -1257.046010
85 840.152174 776.139674
86 -924.047826 840.152174
87 2085.702174 -924.047826
88 184.764674 2085.702174
89 751.527174 184.764674
90 2132.914674 751.527174
91 -53.872826 2132.914674
> 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/76y3n1229016348.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/8esg41229016348.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/9rlhf1229016348.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/10rqzi1229016348.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/11kzlz1229016349.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/12pchb1229016349.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/13qceq1229016349.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/14ebld1229016349.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/15kr3v1229016349.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/16ffvi1229016349.tab")
+ }
>
> system("convert tmp/19mn51229016348.ps tmp/19mn51229016348.png")
> system("convert tmp/2y04h1229016348.ps tmp/2y04h1229016348.png")
> system("convert tmp/3gty81229016348.ps tmp/3gty81229016348.png")
> system("convert tmp/4w5s11229016348.ps tmp/4w5s11229016348.png")
> system("convert tmp/59efa1229016348.ps tmp/59efa1229016348.png")
> system("convert tmp/6gyci1229016348.ps tmp/6gyci1229016348.png")
> system("convert tmp/76y3n1229016348.ps tmp/76y3n1229016348.png")
> system("convert tmp/8esg41229016348.ps tmp/8esg41229016348.png")
> system("convert tmp/9rlhf1229016348.ps tmp/9rlhf1229016348.png")
> system("convert tmp/10rqzi1229016348.ps tmp/10rqzi1229016348.png")
>
>
> proc.time()
user system elapsed
2.832 1.614 3.459