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(14.5,14.8,14.3,14.7,15.3,16,14.4,15.4,13.7,15,14.2,15.5,13.5,15.1,11.9,11.7,14.6,16.3,15.6,16.7,14.1,15,14.9,14.9,14.2,14.6,14.6,15.3,17.2,17.9,15.4,16.4,14.3,15.4,17.5,17.9,14.5,15.9,14.4,13.9,16.6,17.8,16.7,17.9,16.6,17.4,16.9,16.7,15.7,16,16.4,16.6,18.4,19.1,16.9,17.8,16.5,17.2,18.3,18.6,15.1,16.3,15.7,15.1,18.1,19.2,16.8,17.7,18.9,19.1,19,18,18.1,17.5,17.8,17.8,21.5,21.1,17.1,17.2,18.7,19.4,19,19.8,16.4,17.6,16.9,16.2,18.6,19.5,19.3,19.9,19.4,20,17.6,17.3,18.6,18.9,18.1,18.6,20.4,21.4,18.1,18.6,19.6,19.8,19.9,20.8,19.2,19.6,17.8,17.7,19.2,19.8,22,22.2,21.1,20.7,19.5,17.9,22.2,20.9,20.9,21.2,22.2,21.4,23.5,23,21.5,21.3,24.3,23.9,22.8,22.4,20.3,18.3,23.7,22.8,23.3,22.3,19.6,17.8,18,16.4,17.3,16,16.8,16.4,18.2,17.7,16.5,16.6,16,16.2,18.4,18.3),dim=c(2,78),dimnames=list(c('Y','X'),1:78))
> y <- array(NA,dim=c(2,78),dimnames=list(c('Y','X'),1:78))
> 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 14.5 14.8 1 0 0 0 0 0 0 0 0 0 0 1
2 14.3 14.7 0 1 0 0 0 0 0 0 0 0 0 2
3 15.3 16.0 0 0 1 0 0 0 0 0 0 0 0 3
4 14.4 15.4 0 0 0 1 0 0 0 0 0 0 0 4
5 13.7 15.0 0 0 0 0 1 0 0 0 0 0 0 5
6 14.2 15.5 0 0 0 0 0 1 0 0 0 0 0 6
7 13.5 15.1 0 0 0 0 0 0 1 0 0 0 0 7
8 11.9 11.7 0 0 0 0 0 0 0 1 0 0 0 8
9 14.6 16.3 0 0 0 0 0 0 0 0 1 0 0 9
10 15.6 16.7 0 0 0 0 0 0 0 0 0 1 0 10
11 14.1 15.0 0 0 0 0 0 0 0 0 0 0 1 11
12 14.9 14.9 0 0 0 0 0 0 0 0 0 0 0 12
13 14.2 14.6 1 0 0 0 0 0 0 0 0 0 0 13
14 14.6 15.3 0 1 0 0 0 0 0 0 0 0 0 14
15 17.2 17.9 0 0 1 0 0 0 0 0 0 0 0 15
16 15.4 16.4 0 0 0 1 0 0 0 0 0 0 0 16
17 14.3 15.4 0 0 0 0 1 0 0 0 0 0 0 17
18 17.5 17.9 0 0 0 0 0 1 0 0 0 0 0 18
19 14.5 15.9 0 0 0 0 0 0 1 0 0 0 0 19
20 14.4 13.9 0 0 0 0 0 0 0 1 0 0 0 20
21 16.6 17.8 0 0 0 0 0 0 0 0 1 0 0 21
22 16.7 17.9 0 0 0 0 0 0 0 0 0 1 0 22
23 16.6 17.4 0 0 0 0 0 0 0 0 0 0 1 23
24 16.9 16.7 0 0 0 0 0 0 0 0 0 0 0 24
25 15.7 16.0 1 0 0 0 0 0 0 0 0 0 0 25
26 16.4 16.6 0 1 0 0 0 0 0 0 0 0 0 26
27 18.4 19.1 0 0 1 0 0 0 0 0 0 0 0 27
28 16.9 17.8 0 0 0 1 0 0 0 0 0 0 0 28
29 16.5 17.2 0 0 0 0 1 0 0 0 0 0 0 29
30 18.3 18.6 0 0 0 0 0 1 0 0 0 0 0 30
31 15.1 16.3 0 0 0 0 0 0 1 0 0 0 0 31
32 15.7 15.1 0 0 0 0 0 0 0 1 0 0 0 32
33 18.1 19.2 0 0 0 0 0 0 0 0 1 0 0 33
34 16.8 17.7 0 0 0 0 0 0 0 0 0 1 0 34
35 18.9 19.1 0 0 0 0 0 0 0 0 0 0 1 35
36 19.0 18.0 0 0 0 0 0 0 0 0 0 0 0 36
37 18.1 17.5 1 0 0 0 0 0 0 0 0 0 0 37
38 17.8 17.8 0 1 0 0 0 0 0 0 0 0 0 38
39 21.5 21.1 0 0 1 0 0 0 0 0 0 0 0 39
40 17.1 17.2 0 0 0 1 0 0 0 0 0 0 0 40
41 18.7 19.4 0 0 0 0 1 0 0 0 0 0 0 41
42 19.0 19.8 0 0 0 0 0 1 0 0 0 0 0 42
43 16.4 17.6 0 0 0 0 0 0 1 0 0 0 0 43
44 16.9 16.2 0 0 0 0 0 0 0 1 0 0 0 44
45 18.6 19.5 0 0 0 0 0 0 0 0 1 0 0 45
46 19.3 19.9 0 0 0 0 0 0 0 0 0 1 0 46
47 19.4 20.0 0 0 0 0 0 0 0 0 0 0 1 47
48 17.6 17.3 0 0 0 0 0 0 0 0 0 0 0 48
49 18.6 18.9 1 0 0 0 0 0 0 0 0 0 0 49
50 18.1 18.6 0 1 0 0 0 0 0 0 0 0 0 50
51 20.4 21.4 0 0 1 0 0 0 0 0 0 0 0 51
52 18.1 18.6 0 0 0 1 0 0 0 0 0 0 0 52
53 19.6 19.8 0 0 0 0 1 0 0 0 0 0 0 53
54 19.9 20.8 0 0 0 0 0 1 0 0 0 0 0 54
55 19.2 19.6 0 0 0 0 0 0 1 0 0 0 0 55
56 17.8 17.7 0 0 0 0 0 0 0 1 0 0 0 56
57 19.2 19.8 0 0 0 0 0 0 0 0 1 0 0 57
58 22.0 22.2 0 0 0 0 0 0 0 0 0 1 0 58
59 21.1 20.7 0 0 0 0 0 0 0 0 0 0 1 59
60 19.5 17.9 0 0 0 0 0 0 0 0 0 0 0 60
61 22.2 20.9 1 0 0 0 0 0 0 0 0 0 0 61
62 20.9 21.2 0 1 0 0 0 0 0 0 0 0 0 62
63 22.2 21.4 0 0 1 0 0 0 0 0 0 0 0 63
64 23.5 23.0 0 0 0 1 0 0 0 0 0 0 0 64
65 21.5 21.3 0 0 0 0 1 0 0 0 0 0 0 65
66 24.3 23.9 0 0 0 0 0 1 0 0 0 0 0 66
67 22.8 22.4 0 0 0 0 0 0 1 0 0 0 0 67
68 20.3 18.3 0 0 0 0 0 0 0 1 0 0 0 68
69 23.7 22.8 0 0 0 0 0 0 0 0 1 0 0 69
70 23.3 22.3 0 0 0 0 0 0 0 0 0 1 0 70
71 19.6 17.8 0 0 0 0 0 0 0 0 0 0 1 71
72 18.0 16.4 0 0 0 0 0 0 0 0 0 0 0 72
73 17.3 16.0 1 0 0 0 0 0 0 0 0 0 0 73
74 16.8 16.4 0 1 0 0 0 0 0 0 0 0 0 74
75 18.2 17.7 0 0 1 0 0 0 0 0 0 0 0 75
76 16.5 16.6 0 0 0 1 0 0 0 0 0 0 0 76
77 16.0 16.2 0 0 0 0 1 0 0 0 0 0 0 77
78 18.4 18.3 0 0 0 0 0 1 0 0 0 0 0 78
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
-0.57952 1.02934 -0.41123 -0.95414 -0.99064 -1.21392
M5 M6 M7 M8 M9 M10
-1.36022 -1.31061 -1.60788 0.02325 -1.55744 -1.31780
M11 t
-0.85570 0.02067
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.1117045 -0.2266274 0.0002810 0.2561121 1.2456656
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.579524 0.528695 -1.096 0.277125
X 1.029341 0.034375 29.945 < 2e-16 ***
M1 -0.411229 0.256803 -1.601 0.114226
M2 -0.954145 0.257203 -3.710 0.000436 ***
M3 -0.990635 0.270918 -3.657 0.000518 ***
M4 -1.213920 0.259064 -4.686 1.50e-05 ***
M5 -1.360223 0.258245 -5.267 1.73e-06 ***
M6 -1.310615 0.268991 -4.872 7.59e-06 ***
M7 -1.607877 0.269670 -5.962 1.17e-07 ***
M8 0.023252 0.268731 0.087 0.931318
M9 -1.557442 0.280128 -5.560 5.65e-07 ***
M10 -1.317798 0.281753 -4.677 1.55e-05 ***
M11 -0.855700 0.270983 -3.158 0.002426 **
t 0.020666 0.003247 6.364 2.39e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4604 on 64 degrees of freedom
Multiple R-squared: 0.9763, Adjusted R-squared: 0.9715
F-statistic: 203.2 on 13 and 64 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,] 3.136885e-02 6.273771e-02 0.9686311
[2,] 8.332862e-02 1.666572e-01 0.9166714
[3,] 3.587128e-02 7.174255e-02 0.9641287
[4,] 1.568129e-02 3.136258e-02 0.9843187
[5,] 7.719978e-03 1.543996e-02 0.9922800
[6,] 4.196931e-03 8.393862e-03 0.9958031
[7,] 2.725941e-03 5.451883e-03 0.9972741
[8,] 9.828045e-04 1.965609e-03 0.9990172
[9,] 3.344379e-04 6.688758e-04 0.9996656
[10,] 1.831744e-04 3.663488e-04 0.9998168
[11,] 1.100412e-04 2.200823e-04 0.9998900
[12,] 3.800911e-05 7.601822e-05 0.9999620
[13,] 2.585720e-05 5.171440e-05 0.9999741
[14,] 3.208706e-05 6.417412e-05 0.9999679
[15,] 1.656601e-05 3.313202e-05 0.9999834
[16,] 7.009362e-06 1.401872e-05 0.9999930
[17,] 2.263551e-06 4.527102e-06 0.9999977
[18,] 9.532531e-07 1.906506e-06 0.9999990
[19,] 4.479208e-07 8.958416e-07 0.9999996
[20,] 1.716051e-06 3.432103e-06 0.9999983
[21,] 5.005080e-06 1.001016e-05 0.9999950
[22,] 5.276404e-06 1.055281e-05 0.9999947
[23,] 2.546402e-05 5.092803e-05 0.9999745
[24,] 3.671334e-04 7.342667e-04 0.9996329
[25,] 4.047455e-04 8.094911e-04 0.9995953
[26,] 1.891892e-03 3.783783e-03 0.9981081
[27,] 1.206123e-03 2.412245e-03 0.9987939
[28,] 2.103629e-03 4.207259e-03 0.9978964
[29,] 1.243716e-03 2.487432e-03 0.9987563
[30,] 9.035468e-04 1.807094e-03 0.9990965
[31,] 1.577123e-03 3.154246e-03 0.9984229
[32,] 1.234041e-03 2.468083e-03 0.9987660
[33,] 4.022080e-03 8.044159e-03 0.9959779
[34,] 6.799867e-03 1.359973e-02 0.9932001
[35,] 4.291747e-02 8.583493e-02 0.9570825
[36,] 4.115650e-02 8.231300e-02 0.9588435
[37,] 8.316063e-02 1.663213e-01 0.9168394
[38,] 7.817611e-02 1.563522e-01 0.9218239
[39,] 1.357594e-01 2.715188e-01 0.8642406
[40,] 2.266724e-01 4.533448e-01 0.7733276
[41,] 1.699726e-01 3.399451e-01 0.8300274
[42,] 1.151253e-01 2.302506e-01 0.8848747
[43,] 5.281357e-01 9.437287e-01 0.4718643
[44,] 4.762439e-01 9.524878e-01 0.5237561
[45,] 3.640771e-01 7.281541e-01 0.6359229
> postscript(file="/var/www/html/rcomp/tmp/1k64k1258795911.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/25dh31258795911.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/38b4p1258795911.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/4qm4i1258795911.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/5r6wt1258795911.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 = 78
Frequency = 1
1 2 3 4 5
0.2358430455 0.6610267678 0.3387081786 0.2589309971 0.0963046629
6 7 8 9 10
0.0113597395 -0.0003080686 0.2476554417 -0.2272842305 0.1006699578
11 12 13 14 15
-0.1322145809 -0.1056469950 -0.1062819095 0.0954291802 0.0349675628
16 17 18 19 20
-0.0184029070 0.0365752334 0.5929487282 -0.0717738144 0.2351125886
21 22 23 24 25
-0.0192885300 -0.2825321044 -0.3506255921 -0.2064535317 -0.2953521299
26 27 28 29 30
0.3092930389 -0.2482344994 -0.2074731273 0.1357686967 0.4244170615
31 32 33 34 35
-0.1315032440 0.0519105264 -0.2083587504 -0.2246570594 -0.0484980497
36 37 38 39 40
0.3074103270 0.3126435707 0.2260909767 0.5450908057 0.3621382340
41 42 43 44 45
-0.1767741564 -0.3587850007 -0.4176393852 -0.1283574567 -0.2651541008
46 47 48 49 50
-0.2371999125 -0.7228978747 -0.6200442326 -0.8764266497 -0.5453747692
51 52 53 54 55
-1.1117045447 -0.3269319864 0.0634964141 -0.7361189048 0.0756859199
56 57 58 59 60
-1.0203617562 -0.2219494513 -0.1526768446 0.0085704585 0.4143581797
61 62 63 64 65
0.4168986554 -0.6696539386 0.4403023421 0.2959754207 0.1714921146
66 67 68 69 70
0.2249315304 0.5455385923 0.6140406561 0.9420350630 0.7963959631
71 72 73 74 75
1.2456656388 0.2103762528 0.3126754174 -0.0768112557 0.0008701550
76 77 78
-0.3642366311 -0.3268629652 -0.1587531540
> postscript(file="/var/www/html/rcomp/tmp/6q1qw1258795911.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 = 78
Frequency = 1
lag(myerror, k = 1) myerror
0 0.2358430455 NA
1 0.6610267678 0.2358430455
2 0.3387081786 0.6610267678
3 0.2589309971 0.3387081786
4 0.0963046629 0.2589309971
5 0.0113597395 0.0963046629
6 -0.0003080686 0.0113597395
7 0.2476554417 -0.0003080686
8 -0.2272842305 0.2476554417
9 0.1006699578 -0.2272842305
10 -0.1322145809 0.1006699578
11 -0.1056469950 -0.1322145809
12 -0.1062819095 -0.1056469950
13 0.0954291802 -0.1062819095
14 0.0349675628 0.0954291802
15 -0.0184029070 0.0349675628
16 0.0365752334 -0.0184029070
17 0.5929487282 0.0365752334
18 -0.0717738144 0.5929487282
19 0.2351125886 -0.0717738144
20 -0.0192885300 0.2351125886
21 -0.2825321044 -0.0192885300
22 -0.3506255921 -0.2825321044
23 -0.2064535317 -0.3506255921
24 -0.2953521299 -0.2064535317
25 0.3092930389 -0.2953521299
26 -0.2482344994 0.3092930389
27 -0.2074731273 -0.2482344994
28 0.1357686967 -0.2074731273
29 0.4244170615 0.1357686967
30 -0.1315032440 0.4244170615
31 0.0519105264 -0.1315032440
32 -0.2083587504 0.0519105264
33 -0.2246570594 -0.2083587504
34 -0.0484980497 -0.2246570594
35 0.3074103270 -0.0484980497
36 0.3126435707 0.3074103270
37 0.2260909767 0.3126435707
38 0.5450908057 0.2260909767
39 0.3621382340 0.5450908057
40 -0.1767741564 0.3621382340
41 -0.3587850007 -0.1767741564
42 -0.4176393852 -0.3587850007
43 -0.1283574567 -0.4176393852
44 -0.2651541008 -0.1283574567
45 -0.2371999125 -0.2651541008
46 -0.7228978747 -0.2371999125
47 -0.6200442326 -0.7228978747
48 -0.8764266497 -0.6200442326
49 -0.5453747692 -0.8764266497
50 -1.1117045447 -0.5453747692
51 -0.3269319864 -1.1117045447
52 0.0634964141 -0.3269319864
53 -0.7361189048 0.0634964141
54 0.0756859199 -0.7361189048
55 -1.0203617562 0.0756859199
56 -0.2219494513 -1.0203617562
57 -0.1526768446 -0.2219494513
58 0.0085704585 -0.1526768446
59 0.4143581797 0.0085704585
60 0.4168986554 0.4143581797
61 -0.6696539386 0.4168986554
62 0.4403023421 -0.6696539386
63 0.2959754207 0.4403023421
64 0.1714921146 0.2959754207
65 0.2249315304 0.1714921146
66 0.5455385923 0.2249315304
67 0.6140406561 0.5455385923
68 0.9420350630 0.6140406561
69 0.7963959631 0.9420350630
70 1.2456656388 0.7963959631
71 0.2103762528 1.2456656388
72 0.3126754174 0.2103762528
73 -0.0768112557 0.3126754174
74 0.0008701550 -0.0768112557
75 -0.3642366311 0.0008701550
76 -0.3268629652 -0.3642366311
77 -0.1587531540 -0.3268629652
78 NA -0.1587531540
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.6610267678 0.2358430455
[2,] 0.3387081786 0.6610267678
[3,] 0.2589309971 0.3387081786
[4,] 0.0963046629 0.2589309971
[5,] 0.0113597395 0.0963046629
[6,] -0.0003080686 0.0113597395
[7,] 0.2476554417 -0.0003080686
[8,] -0.2272842305 0.2476554417
[9,] 0.1006699578 -0.2272842305
[10,] -0.1322145809 0.1006699578
[11,] -0.1056469950 -0.1322145809
[12,] -0.1062819095 -0.1056469950
[13,] 0.0954291802 -0.1062819095
[14,] 0.0349675628 0.0954291802
[15,] -0.0184029070 0.0349675628
[16,] 0.0365752334 -0.0184029070
[17,] 0.5929487282 0.0365752334
[18,] -0.0717738144 0.5929487282
[19,] 0.2351125886 -0.0717738144
[20,] -0.0192885300 0.2351125886
[21,] -0.2825321044 -0.0192885300
[22,] -0.3506255921 -0.2825321044
[23,] -0.2064535317 -0.3506255921
[24,] -0.2953521299 -0.2064535317
[25,] 0.3092930389 -0.2953521299
[26,] -0.2482344994 0.3092930389
[27,] -0.2074731273 -0.2482344994
[28,] 0.1357686967 -0.2074731273
[29,] 0.4244170615 0.1357686967
[30,] -0.1315032440 0.4244170615
[31,] 0.0519105264 -0.1315032440
[32,] -0.2083587504 0.0519105264
[33,] -0.2246570594 -0.2083587504
[34,] -0.0484980497 -0.2246570594
[35,] 0.3074103270 -0.0484980497
[36,] 0.3126435707 0.3074103270
[37,] 0.2260909767 0.3126435707
[38,] 0.5450908057 0.2260909767
[39,] 0.3621382340 0.5450908057
[40,] -0.1767741564 0.3621382340
[41,] -0.3587850007 -0.1767741564
[42,] -0.4176393852 -0.3587850007
[43,] -0.1283574567 -0.4176393852
[44,] -0.2651541008 -0.1283574567
[45,] -0.2371999125 -0.2651541008
[46,] -0.7228978747 -0.2371999125
[47,] -0.6200442326 -0.7228978747
[48,] -0.8764266497 -0.6200442326
[49,] -0.5453747692 -0.8764266497
[50,] -1.1117045447 -0.5453747692
[51,] -0.3269319864 -1.1117045447
[52,] 0.0634964141 -0.3269319864
[53,] -0.7361189048 0.0634964141
[54,] 0.0756859199 -0.7361189048
[55,] -1.0203617562 0.0756859199
[56,] -0.2219494513 -1.0203617562
[57,] -0.1526768446 -0.2219494513
[58,] 0.0085704585 -0.1526768446
[59,] 0.4143581797 0.0085704585
[60,] 0.4168986554 0.4143581797
[61,] -0.6696539386 0.4168986554
[62,] 0.4403023421 -0.6696539386
[63,] 0.2959754207 0.4403023421
[64,] 0.1714921146 0.2959754207
[65,] 0.2249315304 0.1714921146
[66,] 0.5455385923 0.2249315304
[67,] 0.6140406561 0.5455385923
[68,] 0.9420350630 0.6140406561
[69,] 0.7963959631 0.9420350630
[70,] 1.2456656388 0.7963959631
[71,] 0.2103762528 1.2456656388
[72,] 0.3126754174 0.2103762528
[73,] -0.0768112557 0.3126754174
[74,] 0.0008701550 -0.0768112557
[75,] -0.3642366311 0.0008701550
[76,] -0.3268629652 -0.3642366311
[77,] -0.1587531540 -0.3268629652
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.6610267678 0.2358430455
2 0.3387081786 0.6610267678
3 0.2589309971 0.3387081786
4 0.0963046629 0.2589309971
5 0.0113597395 0.0963046629
6 -0.0003080686 0.0113597395
7 0.2476554417 -0.0003080686
8 -0.2272842305 0.2476554417
9 0.1006699578 -0.2272842305
10 -0.1322145809 0.1006699578
11 -0.1056469950 -0.1322145809
12 -0.1062819095 -0.1056469950
13 0.0954291802 -0.1062819095
14 0.0349675628 0.0954291802
15 -0.0184029070 0.0349675628
16 0.0365752334 -0.0184029070
17 0.5929487282 0.0365752334
18 -0.0717738144 0.5929487282
19 0.2351125886 -0.0717738144
20 -0.0192885300 0.2351125886
21 -0.2825321044 -0.0192885300
22 -0.3506255921 -0.2825321044
23 -0.2064535317 -0.3506255921
24 -0.2953521299 -0.2064535317
25 0.3092930389 -0.2953521299
26 -0.2482344994 0.3092930389
27 -0.2074731273 -0.2482344994
28 0.1357686967 -0.2074731273
29 0.4244170615 0.1357686967
30 -0.1315032440 0.4244170615
31 0.0519105264 -0.1315032440
32 -0.2083587504 0.0519105264
33 -0.2246570594 -0.2083587504
34 -0.0484980497 -0.2246570594
35 0.3074103270 -0.0484980497
36 0.3126435707 0.3074103270
37 0.2260909767 0.3126435707
38 0.5450908057 0.2260909767
39 0.3621382340 0.5450908057
40 -0.1767741564 0.3621382340
41 -0.3587850007 -0.1767741564
42 -0.4176393852 -0.3587850007
43 -0.1283574567 -0.4176393852
44 -0.2651541008 -0.1283574567
45 -0.2371999125 -0.2651541008
46 -0.7228978747 -0.2371999125
47 -0.6200442326 -0.7228978747
48 -0.8764266497 -0.6200442326
49 -0.5453747692 -0.8764266497
50 -1.1117045447 -0.5453747692
51 -0.3269319864 -1.1117045447
52 0.0634964141 -0.3269319864
53 -0.7361189048 0.0634964141
54 0.0756859199 -0.7361189048
55 -1.0203617562 0.0756859199
56 -0.2219494513 -1.0203617562
57 -0.1526768446 -0.2219494513
58 0.0085704585 -0.1526768446
59 0.4143581797 0.0085704585
60 0.4168986554 0.4143581797
61 -0.6696539386 0.4168986554
62 0.4403023421 -0.6696539386
63 0.2959754207 0.4403023421
64 0.1714921146 0.2959754207
65 0.2249315304 0.1714921146
66 0.5455385923 0.2249315304
67 0.6140406561 0.5455385923
68 0.9420350630 0.6140406561
69 0.7963959631 0.9420350630
70 1.2456656388 0.7963959631
71 0.2103762528 1.2456656388
72 0.3126754174 0.2103762528
73 -0.0768112557 0.3126754174
74 0.0008701550 -0.0768112557
75 -0.3642366311 0.0008701550
76 -0.3268629652 -0.3642366311
77 -0.1587531540 -0.3268629652
> 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/75wf71258795911.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/8jnyp1258795911.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/98nl31258795911.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/10dc511258795911.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/112wkk1258795911.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/12058r1258795911.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/13lzw61258795911.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/14q3b31258795911.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/15wamo1258795911.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/16b35y1258795912.tab")
+ }
>
> system("convert tmp/1k64k1258795911.ps tmp/1k64k1258795911.png")
> system("convert tmp/25dh31258795911.ps tmp/25dh31258795911.png")
> system("convert tmp/38b4p1258795911.ps tmp/38b4p1258795911.png")
> system("convert tmp/4qm4i1258795911.ps tmp/4qm4i1258795911.png")
> system("convert tmp/5r6wt1258795911.ps tmp/5r6wt1258795911.png")
> system("convert tmp/6q1qw1258795911.ps tmp/6q1qw1258795911.png")
> system("convert tmp/75wf71258795911.ps tmp/75wf71258795911.png")
> system("convert tmp/8jnyp1258795911.ps tmp/8jnyp1258795911.png")
> system("convert tmp/98nl31258795911.ps tmp/98nl31258795911.png")
> system("convert tmp/10dc511258795911.ps tmp/10dc511258795911.png")
>
>
> proc.time()
user system elapsed
2.611 1.586 3.236