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(61.2
+ ,2.08
+ ,83.9
+ ,10554.27
+ ,62
+ ,2.09
+ ,85.6
+ ,10532.54
+ ,65.1
+ ,2.07
+ ,87.5
+ ,10324.31
+ ,63.2
+ ,2.04
+ ,88.5
+ ,10695.25
+ ,66.3
+ ,2.35
+ ,91
+ ,10827.81
+ ,61.9
+ ,2.33
+ ,90.6
+ ,10872.48
+ ,62.1
+ ,2.37
+ ,91.2
+ ,10971.19
+ ,66.3
+ ,2.59
+ ,93.2
+ ,11145.65
+ ,72
+ ,2.62
+ ,90.1
+ ,11234.68
+ ,65.3
+ ,2.6
+ ,95
+ ,11333.88
+ ,67.6
+ ,2.83
+ ,95.4
+ ,10997.97
+ ,70.5
+ ,2.78
+ ,93.7
+ ,11036.89
+ ,74.2
+ ,3.01
+ ,93.9
+ ,11257.35
+ ,77.8
+ ,3.06
+ ,92.5
+ ,11533.59
+ ,78.5
+ ,3.33
+ ,89.2
+ ,11963.12
+ ,77.8
+ ,3.32
+ ,93.3
+ ,12185.15
+ ,81.4
+ ,3.6
+ ,93
+ ,12377.62
+ ,84.5
+ ,3.57
+ ,96.1
+ ,12512.89
+ ,88
+ ,3.57
+ ,96.7
+ ,12631.48
+ ,93.9
+ ,3.83
+ ,97.6
+ ,12268.53
+ ,98.9
+ ,3.84
+ ,102.6
+ ,12754.8
+ ,96.7
+ ,3.8
+ ,107.6
+ ,13407.75
+ ,98.9
+ ,4.07
+ ,103.5
+ ,13480.21
+ ,102.2
+ ,4.05
+ ,100.8
+ ,13673.28
+ ,105.4
+ ,4.272
+ ,94.5
+ ,13239.71
+ ,105.1
+ ,3.858
+ ,100.1
+ ,13557.69
+ ,116.6
+ ,4.067
+ ,97.4
+ ,13901.28
+ ,112
+ ,3.964
+ ,103
+ ,13200.58
+ ,108.8
+ ,3.782
+ ,100.2
+ ,13406.97
+ ,106.9
+ ,4.114
+ ,100.2
+ ,12538.12
+ ,109.5
+ ,4.009
+ ,99
+ ,12419.57
+ ,106.7
+ ,4.025
+ ,102.4
+ ,12193.88
+ ,118.9
+ ,4.082
+ ,99
+ ,12656.63
+ ,117.5
+ ,4.044
+ ,103.7
+ ,12812.48
+ ,113.7
+ ,3.916
+ ,103.4
+ ,12056.67
+ ,119.6
+ ,4.289
+ ,95.3
+ ,11322.38
+ ,120.6
+ ,4.296
+ ,93.6
+ ,11530.75
+ ,117.5
+ ,4.193
+ ,102.4
+ ,11114.08
+ ,120.3
+ ,3.48
+ ,110.5
+ ,9181.73
+ ,119.8
+ ,2.934
+ ,109.1
+ ,8614.55
+ ,108
+ ,2.221
+ ,100.9
+ ,8595.56
+ ,98.8
+ ,1.211
+ ,108.1
+ ,8396.2
+ ,94.6
+ ,1.28
+ ,105
+ ,7690.5
+ ,84.6
+ ,0.96
+ ,111.5
+ ,7235.47
+ ,84.4
+ ,0.5
+ ,109.5
+ ,7992.12
+ ,79.1
+ ,0.687
+ ,110.5
+ ,8398.37
+ ,73.3
+ ,0.344
+ ,114
+ ,8593
+ ,74.3
+ ,0.346
+ ,108.2
+ ,8679.75
+ ,67.8
+ ,0.334
+ ,110.3
+ ,9374.63
+ ,64.8
+ ,0.34
+ ,111.8
+ ,9634.97
+ ,66.5
+ ,0.328
+ ,107.5
+ ,9857.34
+ ,57.7
+ ,0.344
+ ,114.1
+ ,10238.83
+ ,53.8
+ ,0.341
+ ,113.8
+ ,10433.44
+ ,51.8
+ ,0.32
+ ,114.5
+ ,10471.24
+ ,50.9
+ ,0.314
+ ,114.8
+ ,10214.51
+ ,49
+ ,0.325
+ ,117.8
+ ,10677.52
+ ,48.1
+ ,0.339
+ ,116.7
+ ,11052.15
+ ,42.6
+ ,0.329
+ ,122.8
+ ,10500.19
+ ,40.9
+ ,0.48
+ ,122.3
+ ,10159.27
+ ,43.3
+ ,0.399
+ ,115
+ ,10222.24
+ ,43.7
+ ,0.37
+ ,118.5
+ ,10350.4)
+ ,dim=c(4
+ ,61)
+ ,dimnames=list(c('2JAAR'
+ ,'Eonia'
+ ,'deposits'
+ ,'DowJones')
+ ,1:61))
> y <- array(NA,dim=c(4,61),dimnames=list(c('2JAAR','Eonia','deposits','DowJones'),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
2JAAR Eonia deposits DowJones M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 61.2 2.080 83.9 10554.27 1 0 0 0 0 0 0 0 0 0 0 1
2 62.0 2.090 85.6 10532.54 0 1 0 0 0 0 0 0 0 0 0 2
3 65.1 2.070 87.5 10324.31 0 0 1 0 0 0 0 0 0 0 0 3
4 63.2 2.040 88.5 10695.25 0 0 0 1 0 0 0 0 0 0 0 4
5 66.3 2.350 91.0 10827.81 0 0 0 0 1 0 0 0 0 0 0 5
6 61.9 2.330 90.6 10872.48 0 0 0 0 0 1 0 0 0 0 0 6
7 62.1 2.370 91.2 10971.19 0 0 0 0 0 0 1 0 0 0 0 7
8 66.3 2.590 93.2 11145.65 0 0 0 0 0 0 0 1 0 0 0 8
9 72.0 2.620 90.1 11234.68 0 0 0 0 0 0 0 0 1 0 0 9
10 65.3 2.600 95.0 11333.88 0 0 0 0 0 0 0 0 0 1 0 10
11 67.6 2.830 95.4 10997.97 0 0 0 0 0 0 0 0 0 0 1 11
12 70.5 2.780 93.7 11036.89 0 0 0 0 0 0 0 0 0 0 0 12
13 74.2 3.010 93.9 11257.35 1 0 0 0 0 0 0 0 0 0 0 13
14 77.8 3.060 92.5 11533.59 0 1 0 0 0 0 0 0 0 0 0 14
15 78.5 3.330 89.2 11963.12 0 0 1 0 0 0 0 0 0 0 0 15
16 77.8 3.320 93.3 12185.15 0 0 0 1 0 0 0 0 0 0 0 16
17 81.4 3.600 93.0 12377.62 0 0 0 0 1 0 0 0 0 0 0 17
18 84.5 3.570 96.1 12512.89 0 0 0 0 0 1 0 0 0 0 0 18
19 88.0 3.570 96.7 12631.48 0 0 0 0 0 0 1 0 0 0 0 19
20 93.9 3.830 97.6 12268.53 0 0 0 0 0 0 0 1 0 0 0 20
21 98.9 3.840 102.6 12754.80 0 0 0 0 0 0 0 0 1 0 0 21
22 96.7 3.800 107.6 13407.75 0 0 0 0 0 0 0 0 0 1 0 22
23 98.9 4.070 103.5 13480.21 0 0 0 0 0 0 0 0 0 0 1 23
24 102.2 4.050 100.8 13673.28 0 0 0 0 0 0 0 0 0 0 0 24
25 105.4 4.272 94.5 13239.71 1 0 0 0 0 0 0 0 0 0 0 25
26 105.1 3.858 100.1 13557.69 0 1 0 0 0 0 0 0 0 0 0 26
27 116.6 4.067 97.4 13901.28 0 0 1 0 0 0 0 0 0 0 0 27
28 112.0 3.964 103.0 13200.58 0 0 0 1 0 0 0 0 0 0 0 28
29 108.8 3.782 100.2 13406.97 0 0 0 0 1 0 0 0 0 0 0 29
30 106.9 4.114 100.2 12538.12 0 0 0 0 0 1 0 0 0 0 0 30
31 109.5 4.009 99.0 12419.57 0 0 0 0 0 0 1 0 0 0 0 31
32 106.7 4.025 102.4 12193.88 0 0 0 0 0 0 0 1 0 0 0 32
33 118.9 4.082 99.0 12656.63 0 0 0 0 0 0 0 0 1 0 0 33
34 117.5 4.044 103.7 12812.48 0 0 0 0 0 0 0 0 0 1 0 34
35 113.7 3.916 103.4 12056.67 0 0 0 0 0 0 0 0 0 0 1 35
36 119.6 4.289 95.3 11322.38 0 0 0 0 0 0 0 0 0 0 0 36
37 120.6 4.296 93.6 11530.75 1 0 0 0 0 0 0 0 0 0 0 37
38 117.5 4.193 102.4 11114.08 0 1 0 0 0 0 0 0 0 0 0 38
39 120.3 3.480 110.5 9181.73 0 0 1 0 0 0 0 0 0 0 0 39
40 119.8 2.934 109.1 8614.55 0 0 0 1 0 0 0 0 0 0 0 40
41 108.0 2.221 100.9 8595.56 0 0 0 0 1 0 0 0 0 0 0 41
42 98.8 1.211 108.1 8396.20 0 0 0 0 0 1 0 0 0 0 0 42
43 94.6 1.280 105.0 7690.50 0 0 0 0 0 0 1 0 0 0 0 43
44 84.6 0.960 111.5 7235.47 0 0 0 0 0 0 0 1 0 0 0 44
45 84.4 0.500 109.5 7992.12 0 0 0 0 0 0 0 0 1 0 0 45
46 79.1 0.687 110.5 8398.37 0 0 0 0 0 0 0 0 0 1 0 46
47 73.3 0.344 114.0 8593.00 0 0 0 0 0 0 0 0 0 0 1 47
48 74.3 0.346 108.2 8679.75 0 0 0 0 0 0 0 0 0 0 0 48
49 67.8 0.334 110.3 9374.63 1 0 0 0 0 0 0 0 0 0 0 49
50 64.8 0.340 111.8 9634.97 0 1 0 0 0 0 0 0 0 0 0 50
51 66.5 0.328 107.5 9857.34 0 0 1 0 0 0 0 0 0 0 0 51
52 57.7 0.344 114.1 10238.83 0 0 0 1 0 0 0 0 0 0 0 52
53 53.8 0.341 113.8 10433.44 0 0 0 0 1 0 0 0 0 0 0 53
54 51.8 0.320 114.5 10471.24 0 0 0 0 0 1 0 0 0 0 0 54
55 50.9 0.314 114.8 10214.51 0 0 0 0 0 0 1 0 0 0 0 55
56 49.0 0.325 117.8 10677.52 0 0 0 0 0 0 0 1 0 0 0 56
57 48.1 0.339 116.7 11052.15 0 0 0 0 0 0 0 0 1 0 0 57
58 42.6 0.329 122.8 10500.19 0 0 0 0 0 0 0 0 0 1 0 58
59 40.9 0.480 122.3 10159.27 0 0 0 0 0 0 0 0 0 0 1 59
60 43.3 0.399 115.0 10222.24 0 0 0 0 0 0 0 0 0 0 0 60
61 43.7 0.370 118.5 10350.40 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) Eonia deposits DowJones M1 M2
111.020989 23.361214 -0.272891 -0.007766 -0.094675 4.749453
M3 M4 M5 M6 M7 M8
7.355767 6.810722 5.609829 4.685363 2.638864 0.275055
M9 M10 M11 t
8.582649 5.560427 0.692937 0.801501
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-17.04110 -5.24748 0.02093 5.47531 16.11409
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 111.020989 38.525963 2.882 0.006041 **
Eonia 23.361214 1.456574 16.038 < 2e-16 ***
deposits -0.272891 0.451644 -0.604 0.548733
DowJones -0.007766 0.001141 -6.804 1.99e-08 ***
M1 -0.094675 5.605833 -0.017 0.986600
M2 4.749453 5.864748 0.810 0.422299
M3 7.355767 5.852008 1.257 0.215252
M4 6.810722 5.990779 1.137 0.261612
M5 5.609829 5.842139 0.960 0.342068
M6 4.685363 5.909620 0.793 0.432033
M7 2.638864 5.855147 0.451 0.654376
M8 0.275055 6.093598 0.045 0.964197
M9 8.582649 5.907376 1.453 0.153198
M10 5.560427 6.418455 0.866 0.390911
M11 0.692937 6.320002 0.110 0.913181
t 0.801501 0.226036 3.546 0.000927 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.191 on 45 degrees of freedom
Multiple R-squared: 0.8921, Adjusted R-squared: 0.8562
F-statistic: 24.81 on 15 and 45 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.007399337 0.01479867 0.9926007
[2,] 0.414430785 0.82886157 0.5855692
[3,] 0.357628889 0.71525778 0.6423711
[4,] 0.271353500 0.54270700 0.7286465
[5,] 0.233950557 0.46790111 0.7660494
[6,] 0.152931428 0.30586286 0.8470686
[7,] 0.400169823 0.80033965 0.5998302
[8,] 0.327622218 0.65524444 0.6723778
[9,] 0.338061996 0.67612399 0.6619380
[10,] 0.376768744 0.75353749 0.6232313
[11,] 0.284345500 0.56869100 0.7156545
[12,] 0.580097413 0.83980517 0.4199026
[13,] 0.659746575 0.68050685 0.3402534
[14,] 0.849039390 0.30192122 0.1509606
[15,] 0.817977585 0.36404483 0.1820224
[16,] 0.815348991 0.36930202 0.1846510
[17,] 0.737067121 0.52586576 0.2629329
[18,] 0.745586642 0.50882672 0.2544134
[19,] 0.676076424 0.64784715 0.3239236
[20,] 0.761779685 0.47644063 0.2382203
[21,] 0.813812271 0.37237546 0.1861877
[22,] 0.701838245 0.59632351 0.2981618
[23,] 0.640707602 0.71858480 0.3592924
[24,] 0.733139248 0.53372150 0.2668608
> postscript(file="/var/www/html/rcomp/tmp/1q0lj1293384282.ps",horizontal=F,onefile=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/2q0lj1293384282.ps",horizontal=F,onefile=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/309k41293384282.ps",horizontal=F,onefile=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/409k41293384282.ps",horizontal=F,onefile=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/509k41293384282.ps",horizontal=F,onefile=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
5.74429811 0.96021070 0.02093189 1.71904339 -0.31180916 -3.88385504
7 8 9 10 11 12
-2.44295729 0.08057918 -4.18388033 -6.08834838 -7.59506683 -3.79721993
13 14 15 16 17 18
-4.41038375 -5.86075256 -12.44077015 -10.32040231 -11.44923418 -5.62891922
19 20 21 22 23 24
0.20082155 -0.98396961 -0.18569090 7.20495465 6.60730861 11.02860601
25 26 27 28 29 30
3.24913072 10.97277156 16.11408527 9.75016594 12.04009357 -4.24063252
31 32 33 34 35 36
0.80912582 -1.62729530 2.79805439 6.99947199 4.30396605 -6.53148276
37 38 39 40 41 42
-5.24748439 -12.42145615 -9.16955577 -1.95773817 0.91300698 15.84732270
43 44 45 46 47 48
4.95374077 2.23152302 8.99919011 4.97932296 13.72488750 13.66055950
49 50 51 52 53 54
12.70380116 6.34922644 5.47530876 0.80893115 -1.19205721 -2.09391592
55 56 57 58 59 60
-3.52073084 0.29916271 -7.42767327 -13.09540122 -17.04109532 -14.36046282
61
-12.03936186
> postscript(file="/var/www/html/rcomp/tmp/6ti1p1293384282.ps",horizontal=F,onefile=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 5.74429811 NA
1 0.96021070 5.74429811
2 0.02093189 0.96021070
3 1.71904339 0.02093189
4 -0.31180916 1.71904339
5 -3.88385504 -0.31180916
6 -2.44295729 -3.88385504
7 0.08057918 -2.44295729
8 -4.18388033 0.08057918
9 -6.08834838 -4.18388033
10 -7.59506683 -6.08834838
11 -3.79721993 -7.59506683
12 -4.41038375 -3.79721993
13 -5.86075256 -4.41038375
14 -12.44077015 -5.86075256
15 -10.32040231 -12.44077015
16 -11.44923418 -10.32040231
17 -5.62891922 -11.44923418
18 0.20082155 -5.62891922
19 -0.98396961 0.20082155
20 -0.18569090 -0.98396961
21 7.20495465 -0.18569090
22 6.60730861 7.20495465
23 11.02860601 6.60730861
24 3.24913072 11.02860601
25 10.97277156 3.24913072
26 16.11408527 10.97277156
27 9.75016594 16.11408527
28 12.04009357 9.75016594
29 -4.24063252 12.04009357
30 0.80912582 -4.24063252
31 -1.62729530 0.80912582
32 2.79805439 -1.62729530
33 6.99947199 2.79805439
34 4.30396605 6.99947199
35 -6.53148276 4.30396605
36 -5.24748439 -6.53148276
37 -12.42145615 -5.24748439
38 -9.16955577 -12.42145615
39 -1.95773817 -9.16955577
40 0.91300698 -1.95773817
41 15.84732270 0.91300698
42 4.95374077 15.84732270
43 2.23152302 4.95374077
44 8.99919011 2.23152302
45 4.97932296 8.99919011
46 13.72488750 4.97932296
47 13.66055950 13.72488750
48 12.70380116 13.66055950
49 6.34922644 12.70380116
50 5.47530876 6.34922644
51 0.80893115 5.47530876
52 -1.19205721 0.80893115
53 -2.09391592 -1.19205721
54 -3.52073084 -2.09391592
55 0.29916271 -3.52073084
56 -7.42767327 0.29916271
57 -13.09540122 -7.42767327
58 -17.04109532 -13.09540122
59 -14.36046282 -17.04109532
60 -12.03936186 -14.36046282
61 NA -12.03936186
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.96021070 5.74429811
[2,] 0.02093189 0.96021070
[3,] 1.71904339 0.02093189
[4,] -0.31180916 1.71904339
[5,] -3.88385504 -0.31180916
[6,] -2.44295729 -3.88385504
[7,] 0.08057918 -2.44295729
[8,] -4.18388033 0.08057918
[9,] -6.08834838 -4.18388033
[10,] -7.59506683 -6.08834838
[11,] -3.79721993 -7.59506683
[12,] -4.41038375 -3.79721993
[13,] -5.86075256 -4.41038375
[14,] -12.44077015 -5.86075256
[15,] -10.32040231 -12.44077015
[16,] -11.44923418 -10.32040231
[17,] -5.62891922 -11.44923418
[18,] 0.20082155 -5.62891922
[19,] -0.98396961 0.20082155
[20,] -0.18569090 -0.98396961
[21,] 7.20495465 -0.18569090
[22,] 6.60730861 7.20495465
[23,] 11.02860601 6.60730861
[24,] 3.24913072 11.02860601
[25,] 10.97277156 3.24913072
[26,] 16.11408527 10.97277156
[27,] 9.75016594 16.11408527
[28,] 12.04009357 9.75016594
[29,] -4.24063252 12.04009357
[30,] 0.80912582 -4.24063252
[31,] -1.62729530 0.80912582
[32,] 2.79805439 -1.62729530
[33,] 6.99947199 2.79805439
[34,] 4.30396605 6.99947199
[35,] -6.53148276 4.30396605
[36,] -5.24748439 -6.53148276
[37,] -12.42145615 -5.24748439
[38,] -9.16955577 -12.42145615
[39,] -1.95773817 -9.16955577
[40,] 0.91300698 -1.95773817
[41,] 15.84732270 0.91300698
[42,] 4.95374077 15.84732270
[43,] 2.23152302 4.95374077
[44,] 8.99919011 2.23152302
[45,] 4.97932296 8.99919011
[46,] 13.72488750 4.97932296
[47,] 13.66055950 13.72488750
[48,] 12.70380116 13.66055950
[49,] 6.34922644 12.70380116
[50,] 5.47530876 6.34922644
[51,] 0.80893115 5.47530876
[52,] -1.19205721 0.80893115
[53,] -2.09391592 -1.19205721
[54,] -3.52073084 -2.09391592
[55,] 0.29916271 -3.52073084
[56,] -7.42767327 0.29916271
[57,] -13.09540122 -7.42767327
[58,] -17.04109532 -13.09540122
[59,] -14.36046282 -17.04109532
[60,] -12.03936186 -14.36046282
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.96021070 5.74429811
2 0.02093189 0.96021070
3 1.71904339 0.02093189
4 -0.31180916 1.71904339
5 -3.88385504 -0.31180916
6 -2.44295729 -3.88385504
7 0.08057918 -2.44295729
8 -4.18388033 0.08057918
9 -6.08834838 -4.18388033
10 -7.59506683 -6.08834838
11 -3.79721993 -7.59506683
12 -4.41038375 -3.79721993
13 -5.86075256 -4.41038375
14 -12.44077015 -5.86075256
15 -10.32040231 -12.44077015
16 -11.44923418 -10.32040231
17 -5.62891922 -11.44923418
18 0.20082155 -5.62891922
19 -0.98396961 0.20082155
20 -0.18569090 -0.98396961
21 7.20495465 -0.18569090
22 6.60730861 7.20495465
23 11.02860601 6.60730861
24 3.24913072 11.02860601
25 10.97277156 3.24913072
26 16.11408527 10.97277156
27 9.75016594 16.11408527
28 12.04009357 9.75016594
29 -4.24063252 12.04009357
30 0.80912582 -4.24063252
31 -1.62729530 0.80912582
32 2.79805439 -1.62729530
33 6.99947199 2.79805439
34 4.30396605 6.99947199
35 -6.53148276 4.30396605
36 -5.24748439 -6.53148276
37 -12.42145615 -5.24748439
38 -9.16955577 -12.42145615
39 -1.95773817 -9.16955577
40 0.91300698 -1.95773817
41 15.84732270 0.91300698
42 4.95374077 15.84732270
43 2.23152302 4.95374077
44 8.99919011 2.23152302
45 4.97932296 8.99919011
46 13.72488750 4.97932296
47 13.66055950 13.72488750
48 12.70380116 13.66055950
49 6.34922644 12.70380116
50 5.47530876 6.34922644
51 0.80893115 5.47530876
52 -1.19205721 0.80893115
53 -2.09391592 -1.19205721
54 -3.52073084 -2.09391592
55 0.29916271 -3.52073084
56 -7.42767327 0.29916271
57 -13.09540122 -7.42767327
58 -17.04109532 -13.09540122
59 -14.36046282 -17.04109532
60 -12.03936186 -14.36046282
> 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/7491a1293384282.ps",horizontal=F,onefile=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/8491a1293384282.ps",horizontal=F,onefile=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/9491a1293384282.ps",horizontal=F,onefile=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/10x10d1293384282.ps",horizontal=F,onefile=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/11ijz11293384282.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/123kx61293384282.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/130uvx1293384282.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/14d4wg1293384283.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/15hnd41293384283.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/162nts1293384283.tab")
+ }
> try(system("convert tmp/1q0lj1293384282.ps tmp/1q0lj1293384282.png",intern=TRUE))
character(0)
> try(system("convert tmp/2q0lj1293384282.ps tmp/2q0lj1293384282.png",intern=TRUE))
character(0)
> try(system("convert tmp/309k41293384282.ps tmp/309k41293384282.png",intern=TRUE))
character(0)
> try(system("convert tmp/409k41293384282.ps tmp/409k41293384282.png",intern=TRUE))
character(0)
> try(system("convert tmp/509k41293384282.ps tmp/509k41293384282.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ti1p1293384282.ps tmp/6ti1p1293384282.png",intern=TRUE))
character(0)
> try(system("convert tmp/7491a1293384282.ps tmp/7491a1293384282.png",intern=TRUE))
character(0)
> try(system("convert tmp/8491a1293384282.ps tmp/8491a1293384282.png",intern=TRUE))
character(0)
> try(system("convert tmp/9491a1293384282.ps tmp/9491a1293384282.png",intern=TRUE))
character(0)
> try(system("convert tmp/10x10d1293384282.ps tmp/10x10d1293384282.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.417 1.630 5.912