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(2.3
+ ,2.0
+ ,1.9
+ ,2.3
+ ,2.7
+ ,0
+ ,2.8
+ ,2.3
+ ,2.0
+ ,1.9
+ ,2.3
+ ,0
+ ,2.4
+ ,2.8
+ ,2.3
+ ,2.0
+ ,1.9
+ ,0
+ ,2.3
+ ,2.4
+ ,2.8
+ ,2.3
+ ,2.0
+ ,0
+ ,2.7
+ ,2.3
+ ,2.4
+ ,2.8
+ ,2.3
+ ,0
+ ,2.7
+ ,2.7
+ ,2.3
+ ,2.4
+ ,2.8
+ ,0
+ ,2.9
+ ,2.7
+ ,2.7
+ ,2.3
+ ,2.4
+ ,0
+ ,3.0
+ ,2.9
+ ,2.7
+ ,2.7
+ ,2.3
+ ,0
+ ,2.2
+ ,3.0
+ ,2.9
+ ,2.7
+ ,2.7
+ ,0
+ ,2.3
+ ,2.2
+ ,3.0
+ ,2.9
+ ,2.7
+ ,0
+ ,2.8
+ ,2.3
+ ,2.2
+ ,3.0
+ ,2.9
+ ,0
+ ,2.8
+ ,2.8
+ ,2.3
+ ,2.2
+ ,3.0
+ ,0
+ ,2.8
+ ,2.8
+ ,2.8
+ ,2.3
+ ,2.2
+ ,0
+ ,2.2
+ ,2.8
+ ,2.8
+ ,2.8
+ ,2.3
+ ,0
+ ,2.6
+ ,2.2
+ ,2.8
+ ,2.8
+ ,2.8
+ ,0
+ ,2.8
+ ,2.6
+ ,2.2
+ ,2.8
+ ,2.8
+ ,0
+ ,2.5
+ ,2.8
+ ,2.6
+ ,2.2
+ ,2.8
+ ,0
+ ,2.4
+ ,2.5
+ ,2.8
+ ,2.6
+ ,2.2
+ ,0
+ ,2.3
+ ,2.4
+ ,2.5
+ ,2.8
+ ,2.6
+ ,0
+ ,1.9
+ ,2.3
+ ,2.4
+ ,2.5
+ ,2.8
+ ,0
+ ,1.7
+ ,1.9
+ ,2.3
+ ,2.4
+ ,2.5
+ ,0
+ ,2.0
+ ,1.7
+ ,1.9
+ ,2.3
+ ,2.4
+ ,0
+ ,2.1
+ ,2.0
+ ,1.7
+ ,1.9
+ ,2.3
+ ,0
+ ,1.7
+ ,2.1
+ ,2.0
+ ,1.7
+ ,1.9
+ ,0
+ ,1.8
+ ,1.7
+ ,2.1
+ ,2.0
+ ,1.7
+ ,0
+ ,1.8
+ ,1.8
+ ,1.7
+ ,2.1
+ ,2.0
+ ,0
+ ,1.8
+ ,1.8
+ ,1.8
+ ,1.7
+ ,2.1
+ ,0
+ ,1.3
+ ,1.8
+ ,1.8
+ ,1.8
+ ,1.7
+ ,0
+ ,1.3
+ ,1.3
+ ,1.8
+ ,1.8
+ ,1.8
+ ,0
+ ,1.3
+ ,1.3
+ ,1.3
+ ,1.8
+ ,1.8
+ ,0
+ ,1.2
+ ,1.3
+ ,1.3
+ ,1.3
+ ,1.8
+ ,1
+ ,1.4
+ ,1.2
+ ,1.3
+ ,1.3
+ ,1.3
+ ,1
+ ,2.2
+ ,1.4
+ ,1.2
+ ,1.3
+ ,1.3
+ ,1
+ ,2.9
+ ,2.2
+ ,1.4
+ ,1.2
+ ,1.3
+ ,1
+ ,3.1
+ ,2.9
+ ,2.2
+ ,1.4
+ ,1.2
+ ,1
+ ,3.5
+ ,3.1
+ ,2.9
+ ,2.2
+ ,1.4
+ ,1
+ ,3.6
+ ,3.5
+ ,3.1
+ ,2.9
+ ,2.2
+ ,1
+ ,4.4
+ ,3.6
+ ,3.5
+ ,3.1
+ ,2.9
+ ,1
+ ,4.1
+ ,4.4
+ ,3.6
+ ,3.5
+ ,3.1
+ ,1
+ ,5.1
+ ,4.1
+ ,4.4
+ ,3.6
+ ,3.5
+ ,1
+ ,5.8
+ ,5.1
+ ,4.1
+ ,4.4
+ ,3.6
+ ,1
+ ,5.9
+ ,5.8
+ ,5.1
+ ,4.1
+ ,4.4
+ ,1
+ ,5.4
+ ,5.9
+ ,5.8
+ ,5.1
+ ,4.1
+ ,1
+ ,5.5
+ ,5.4
+ ,5.9
+ ,5.8
+ ,5.1
+ ,1
+ ,4.8
+ ,5.5
+ ,5.4
+ ,5.9
+ ,5.8
+ ,1
+ ,3.2
+ ,4.8
+ ,5.5
+ ,5.4
+ ,5.9
+ ,1
+ ,2.7
+ ,3.2
+ ,4.8
+ ,5.5
+ ,5.4
+ ,1
+ ,2.1
+ ,2.7
+ ,3.2
+ ,4.8
+ ,5.5
+ ,1
+ ,1.9
+ ,2.1
+ ,2.7
+ ,3.2
+ ,4.8
+ ,1
+ ,0.6
+ ,1.9
+ ,2.1
+ ,2.7
+ ,3.2
+ ,1
+ ,0.7
+ ,0.6
+ ,1.9
+ ,2.1
+ ,2.7
+ ,1
+ ,-0.2
+ ,0.7
+ ,0.6
+ ,1.9
+ ,2.1
+ ,1
+ ,-1.0
+ ,-0.2
+ ,0.7
+ ,0.6
+ ,1.9
+ ,1
+ ,-1.7
+ ,-1.0
+ ,-0.2
+ ,0.7
+ ,0.6
+ ,1
+ ,-0.7
+ ,-1.7
+ ,-1.0
+ ,-0.2
+ ,0.7
+ ,1
+ ,-1.0
+ ,-0.7
+ ,-1.7
+ ,-1.0
+ ,-0.2
+ ,1)
+ ,dim=c(6
+ ,56)
+ ,dimnames=list(c('Inflatie'
+ ,'yt-1'
+ ,'yt-2'
+ ,'yt-3'
+ ,'yt-4'
+ ,'Kredietcrisis')
+ ,1:56))
> y <- array(NA,dim=c(6,56),dimnames=list(c('Inflatie','yt-1','yt-2','yt-3','yt-4','Kredietcrisis'),1:56))
> 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
Inflatie yt-1 yt-2 yt-3 yt-4 Kredietcrisis M1 M2 M3 M4 M5 M6 M7 M8 M9 M10
1 2.3 2.0 1.9 2.3 2.7 0 1 0 0 0 0 0 0 0 0 0
2 2.8 2.3 2.0 1.9 2.3 0 0 1 0 0 0 0 0 0 0 0
3 2.4 2.8 2.3 2.0 1.9 0 0 0 1 0 0 0 0 0 0 0
4 2.3 2.4 2.8 2.3 2.0 0 0 0 0 1 0 0 0 0 0 0
5 2.7 2.3 2.4 2.8 2.3 0 0 0 0 0 1 0 0 0 0 0
6 2.7 2.7 2.3 2.4 2.8 0 0 0 0 0 0 1 0 0 0 0
7 2.9 2.7 2.7 2.3 2.4 0 0 0 0 0 0 0 1 0 0 0
8 3.0 2.9 2.7 2.7 2.3 0 0 0 0 0 0 0 0 1 0 0
9 2.2 3.0 2.9 2.7 2.7 0 0 0 0 0 0 0 0 0 1 0
10 2.3 2.2 3.0 2.9 2.7 0 0 0 0 0 0 0 0 0 0 1
11 2.8 2.3 2.2 3.0 2.9 0 0 0 0 0 0 0 0 0 0 0
12 2.8 2.8 2.3 2.2 3.0 0 0 0 0 0 0 0 0 0 0 0
13 2.8 2.8 2.8 2.3 2.2 0 1 0 0 0 0 0 0 0 0 0
14 2.2 2.8 2.8 2.8 2.3 0 0 1 0 0 0 0 0 0 0 0
15 2.6 2.2 2.8 2.8 2.8 0 0 0 1 0 0 0 0 0 0 0
16 2.8 2.6 2.2 2.8 2.8 0 0 0 0 1 0 0 0 0 0 0
17 2.5 2.8 2.6 2.2 2.8 0 0 0 0 0 1 0 0 0 0 0
18 2.4 2.5 2.8 2.6 2.2 0 0 0 0 0 0 1 0 0 0 0
19 2.3 2.4 2.5 2.8 2.6 0 0 0 0 0 0 0 1 0 0 0
20 1.9 2.3 2.4 2.5 2.8 0 0 0 0 0 0 0 0 1 0 0
21 1.7 1.9 2.3 2.4 2.5 0 0 0 0 0 0 0 0 0 1 0
22 2.0 1.7 1.9 2.3 2.4 0 0 0 0 0 0 0 0 0 0 1
23 2.1 2.0 1.7 1.9 2.3 0 0 0 0 0 0 0 0 0 0 0
24 1.7 2.1 2.0 1.7 1.9 0 0 0 0 0 0 0 0 0 0 0
25 1.8 1.7 2.1 2.0 1.7 0 1 0 0 0 0 0 0 0 0 0
26 1.8 1.8 1.7 2.1 2.0 0 0 1 0 0 0 0 0 0 0 0
27 1.8 1.8 1.8 1.7 2.1 0 0 0 1 0 0 0 0 0 0 0
28 1.3 1.8 1.8 1.8 1.7 0 0 0 0 1 0 0 0 0 0 0
29 1.3 1.3 1.8 1.8 1.8 0 0 0 0 0 1 0 0 0 0 0
30 1.3 1.3 1.3 1.8 1.8 0 0 0 0 0 0 1 0 0 0 0
31 1.2 1.3 1.3 1.3 1.8 1 0 0 0 0 0 0 1 0 0 0
32 1.4 1.2 1.3 1.3 1.3 1 0 0 0 0 0 0 0 1 0 0
33 2.2 1.4 1.2 1.3 1.3 1 0 0 0 0 0 0 0 0 1 0
34 2.9 2.2 1.4 1.2 1.3 1 0 0 0 0 0 0 0 0 0 1
35 3.1 2.9 2.2 1.4 1.2 1 0 0 0 0 0 0 0 0 0 0
36 3.5 3.1 2.9 2.2 1.4 1 0 0 0 0 0 0 0 0 0 0
37 3.6 3.5 3.1 2.9 2.2 1 1 0 0 0 0 0 0 0 0 0
38 4.4 3.6 3.5 3.1 2.9 1 0 1 0 0 0 0 0 0 0 0
39 4.1 4.4 3.6 3.5 3.1 1 0 0 1 0 0 0 0 0 0 0
40 5.1 4.1 4.4 3.6 3.5 1 0 0 0 1 0 0 0 0 0 0
41 5.8 5.1 4.1 4.4 3.6 1 0 0 0 0 1 0 0 0 0 0
42 5.9 5.8 5.1 4.1 4.4 1 0 0 0 0 0 1 0 0 0 0
43 5.4 5.9 5.8 5.1 4.1 1 0 0 0 0 0 0 1 0 0 0
44 5.5 5.4 5.9 5.8 5.1 1 0 0 0 0 0 0 0 1 0 0
45 4.8 5.5 5.4 5.9 5.8 1 0 0 0 0 0 0 0 0 1 0
46 3.2 4.8 5.5 5.4 5.9 1 0 0 0 0 0 0 0 0 0 1
47 2.7 3.2 4.8 5.5 5.4 1 0 0 0 0 0 0 0 0 0 0
48 2.1 2.7 3.2 4.8 5.5 1 0 0 0 0 0 0 0 0 0 0
49 1.9 2.1 2.7 3.2 4.8 1 1 0 0 0 0 0 0 0 0 0
50 0.6 1.9 2.1 2.7 3.2 1 0 1 0 0 0 0 0 0 0 0
51 0.7 0.6 1.9 2.1 2.7 1 0 0 1 0 0 0 0 0 0 0
52 -0.2 0.7 0.6 1.9 2.1 1 0 0 0 1 0 0 0 0 0 0
53 -1.0 -0.2 0.7 0.6 1.9 1 0 0 0 0 1 0 0 0 0 0
54 -1.7 -1.0 -0.2 0.7 0.6 1 0 0 0 0 0 1 0 0 0 0
55 -0.7 -1.7 -1.0 -0.2 0.7 1 0 0 0 0 0 0 1 0 0 0
56 -1.0 -0.7 -1.7 -1.0 -0.2 1 0 0 0 0 0 0 0 1 0 0
M11 t
1 0 1
2 0 2
3 0 3
4 0 4
5 0 5
6 0 6
7 0 7
8 0 8
9 0 9
10 0 10
11 1 11
12 0 12
13 0 13
14 0 14
15 0 15
16 0 16
17 0 17
18 0 18
19 0 19
20 0 20
21 0 21
22 0 22
23 1 23
24 0 24
25 0 25
26 0 26
27 0 27
28 0 28
29 0 29
30 0 30
31 0 31
32 0 32
33 0 33
34 0 34
35 1 35
36 0 36
37 0 37
38 0 38
39 0 39
40 0 40
41 0 41
42 0 42
43 0 43
44 0 44
45 0 45
46 0 46
47 1 47
48 0 48
49 0 49
50 0 50
51 0 51
52 0 52
53 0 53
54 0 54
55 0 55
56 0 56
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `yt-1` `yt-2` `yt-3` `yt-4`
0.543937 1.031906 0.003116 0.062595 -0.223419
Kredietcrisis M1 M2 M3 M4
0.554773 0.138013 -0.063344 0.041297 0.016058
M5 M6 M7 M8 M9
0.118222 -0.026689 0.120851 -0.036100 -0.136471
M10 M11 t
-0.002260 0.193711 -0.019209
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.00298 -0.26340 0.00827 0.30418 1.06568
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.543937 0.381080 1.427 0.1616
`yt-1` 1.031906 0.160828 6.416 1.53e-07 ***
`yt-2` 0.003116 0.242829 0.013 0.9898
`yt-3` 0.062595 0.256666 0.244 0.8086
`yt-4` -0.223419 0.176370 -1.267 0.2130
Kredietcrisis 0.554773 0.334214 1.660 0.1052
M1 0.138013 0.368967 0.374 0.7104
M2 -0.063344 0.368265 -0.172 0.8643
M3 0.041297 0.370395 0.111 0.9118
M4 0.016058 0.371267 0.043 0.9657
M5 0.118222 0.368993 0.320 0.7504
M6 -0.026689 0.370597 -0.072 0.9430
M7 0.120851 0.373598 0.323 0.7481
M8 -0.036100 0.371068 -0.097 0.9230
M9 -0.136471 0.388445 -0.351 0.7273
M10 -0.002260 0.390185 -0.006 0.9954
M11 0.193711 0.389660 0.497 0.6220
t -0.019209 0.010719 -1.792 0.0811 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5451 on 38 degrees of freedom
Multiple R-squared: 0.9192, Adjusted R-squared: 0.883
F-statistic: 25.41 on 17 and 38 DF, p-value: 1.181e-15
> 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.3864767542 0.7729535084 0.6135232
[2,] 0.2364679532 0.4729359064 0.7635320
[3,] 0.1283087376 0.2566174753 0.8716913
[4,] 0.0707207792 0.1414415583 0.9292792
[5,] 0.0320473324 0.0640946649 0.9679527
[6,] 0.0132101222 0.0264202443 0.9867899
[7,] 0.0049333050 0.0098666100 0.9950667
[8,] 0.0027531405 0.0055062811 0.9972469
[9,] 0.0009502512 0.0019005024 0.9990497
[10,] 0.0002836344 0.0005672688 0.9997164
[11,] 0.0001984477 0.0003968954 0.9998016
[12,] 0.0009441303 0.0018882606 0.9990559
[13,] 0.0115541614 0.0231083229 0.9884458
[14,] 0.0057027747 0.0114055495 0.9942972
[15,] 0.0042461651 0.0084923302 0.9957538
> postscript(file="/var/www/html/rcomp/tmp/1uibo1259252857.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/2vu9h1259252857.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/3il4t1259252857.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/4rz0t1259252857.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/5mkpn1259252857.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 = 56
Frequency = 1
1 2 3 4 5 6
0.02678960 0.37314283 -0.72480486 -0.36558956 0.09162021 -0.01996353
7 8 9 10 11 12
-0.03264924 -0.01025074 -0.70511687 0.09257410 0.35353867 0.12261151
13 14 15 16 17 18
-0.18274586 -0.57113505 0.47428517 0.30783983 -0.24518635 -0.03120761
19 20 21 22 23 24
-0.07856474 -0.13544079 0.13644692 0.51298887 0.12997437 -0.23808006
25 26 27 28 29 30
0.09210380 0.27149213 0.23312752 -0.31805215 0.13728719 0.30296470
31 32 33 34 35 36
-0.44884194 -0.08120139 0.63230904 0.39741745 -0.33903222 0.05993303
37 38 39 40 41 42
-0.23733824 0.82266550 -0.36895775 1.06567687 0.62401632 0.36019953
43 44 45 46 47 48
-0.50312393 0.46827964 -0.06363908 -1.00298042 -0.14448082 0.05553552
49 50 51 52 53 54
0.30119070 -0.89616540 0.38634991 -0.68987499 -0.60773738 -0.61199309
55 56
1.06317984 -0.24138671
> postscript(file="/var/www/html/rcomp/tmp/652i91259252857.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 = 56
Frequency = 1
lag(myerror, k = 1) myerror
0 0.02678960 NA
1 0.37314283 0.02678960
2 -0.72480486 0.37314283
3 -0.36558956 -0.72480486
4 0.09162021 -0.36558956
5 -0.01996353 0.09162021
6 -0.03264924 -0.01996353
7 -0.01025074 -0.03264924
8 -0.70511687 -0.01025074
9 0.09257410 -0.70511687
10 0.35353867 0.09257410
11 0.12261151 0.35353867
12 -0.18274586 0.12261151
13 -0.57113505 -0.18274586
14 0.47428517 -0.57113505
15 0.30783983 0.47428517
16 -0.24518635 0.30783983
17 -0.03120761 -0.24518635
18 -0.07856474 -0.03120761
19 -0.13544079 -0.07856474
20 0.13644692 -0.13544079
21 0.51298887 0.13644692
22 0.12997437 0.51298887
23 -0.23808006 0.12997437
24 0.09210380 -0.23808006
25 0.27149213 0.09210380
26 0.23312752 0.27149213
27 -0.31805215 0.23312752
28 0.13728719 -0.31805215
29 0.30296470 0.13728719
30 -0.44884194 0.30296470
31 -0.08120139 -0.44884194
32 0.63230904 -0.08120139
33 0.39741745 0.63230904
34 -0.33903222 0.39741745
35 0.05993303 -0.33903222
36 -0.23733824 0.05993303
37 0.82266550 -0.23733824
38 -0.36895775 0.82266550
39 1.06567687 -0.36895775
40 0.62401632 1.06567687
41 0.36019953 0.62401632
42 -0.50312393 0.36019953
43 0.46827964 -0.50312393
44 -0.06363908 0.46827964
45 -1.00298042 -0.06363908
46 -0.14448082 -1.00298042
47 0.05553552 -0.14448082
48 0.30119070 0.05553552
49 -0.89616540 0.30119070
50 0.38634991 -0.89616540
51 -0.68987499 0.38634991
52 -0.60773738 -0.68987499
53 -0.61199309 -0.60773738
54 1.06317984 -0.61199309
55 -0.24138671 1.06317984
56 NA -0.24138671
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.37314283 0.02678960
[2,] -0.72480486 0.37314283
[3,] -0.36558956 -0.72480486
[4,] 0.09162021 -0.36558956
[5,] -0.01996353 0.09162021
[6,] -0.03264924 -0.01996353
[7,] -0.01025074 -0.03264924
[8,] -0.70511687 -0.01025074
[9,] 0.09257410 -0.70511687
[10,] 0.35353867 0.09257410
[11,] 0.12261151 0.35353867
[12,] -0.18274586 0.12261151
[13,] -0.57113505 -0.18274586
[14,] 0.47428517 -0.57113505
[15,] 0.30783983 0.47428517
[16,] -0.24518635 0.30783983
[17,] -0.03120761 -0.24518635
[18,] -0.07856474 -0.03120761
[19,] -0.13544079 -0.07856474
[20,] 0.13644692 -0.13544079
[21,] 0.51298887 0.13644692
[22,] 0.12997437 0.51298887
[23,] -0.23808006 0.12997437
[24,] 0.09210380 -0.23808006
[25,] 0.27149213 0.09210380
[26,] 0.23312752 0.27149213
[27,] -0.31805215 0.23312752
[28,] 0.13728719 -0.31805215
[29,] 0.30296470 0.13728719
[30,] -0.44884194 0.30296470
[31,] -0.08120139 -0.44884194
[32,] 0.63230904 -0.08120139
[33,] 0.39741745 0.63230904
[34,] -0.33903222 0.39741745
[35,] 0.05993303 -0.33903222
[36,] -0.23733824 0.05993303
[37,] 0.82266550 -0.23733824
[38,] -0.36895775 0.82266550
[39,] 1.06567687 -0.36895775
[40,] 0.62401632 1.06567687
[41,] 0.36019953 0.62401632
[42,] -0.50312393 0.36019953
[43,] 0.46827964 -0.50312393
[44,] -0.06363908 0.46827964
[45,] -1.00298042 -0.06363908
[46,] -0.14448082 -1.00298042
[47,] 0.05553552 -0.14448082
[48,] 0.30119070 0.05553552
[49,] -0.89616540 0.30119070
[50,] 0.38634991 -0.89616540
[51,] -0.68987499 0.38634991
[52,] -0.60773738 -0.68987499
[53,] -0.61199309 -0.60773738
[54,] 1.06317984 -0.61199309
[55,] -0.24138671 1.06317984
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.37314283 0.02678960
2 -0.72480486 0.37314283
3 -0.36558956 -0.72480486
4 0.09162021 -0.36558956
5 -0.01996353 0.09162021
6 -0.03264924 -0.01996353
7 -0.01025074 -0.03264924
8 -0.70511687 -0.01025074
9 0.09257410 -0.70511687
10 0.35353867 0.09257410
11 0.12261151 0.35353867
12 -0.18274586 0.12261151
13 -0.57113505 -0.18274586
14 0.47428517 -0.57113505
15 0.30783983 0.47428517
16 -0.24518635 0.30783983
17 -0.03120761 -0.24518635
18 -0.07856474 -0.03120761
19 -0.13544079 -0.07856474
20 0.13644692 -0.13544079
21 0.51298887 0.13644692
22 0.12997437 0.51298887
23 -0.23808006 0.12997437
24 0.09210380 -0.23808006
25 0.27149213 0.09210380
26 0.23312752 0.27149213
27 -0.31805215 0.23312752
28 0.13728719 -0.31805215
29 0.30296470 0.13728719
30 -0.44884194 0.30296470
31 -0.08120139 -0.44884194
32 0.63230904 -0.08120139
33 0.39741745 0.63230904
34 -0.33903222 0.39741745
35 0.05993303 -0.33903222
36 -0.23733824 0.05993303
37 0.82266550 -0.23733824
38 -0.36895775 0.82266550
39 1.06567687 -0.36895775
40 0.62401632 1.06567687
41 0.36019953 0.62401632
42 -0.50312393 0.36019953
43 0.46827964 -0.50312393
44 -0.06363908 0.46827964
45 -1.00298042 -0.06363908
46 -0.14448082 -1.00298042
47 0.05553552 -0.14448082
48 0.30119070 0.05553552
49 -0.89616540 0.30119070
50 0.38634991 -0.89616540
51 -0.68987499 0.38634991
52 -0.60773738 -0.68987499
53 -0.61199309 -0.60773738
54 1.06317984 -0.61199309
55 -0.24138671 1.06317984
> 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/7avku1259252857.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/8htdn1259252857.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/9dwed1259252857.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/10eepz1259252857.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/11g5if1259252857.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/1290zj1259252857.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/13x00h1259252857.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/14sm371259252857.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/15wkqp1259252857.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/160dpj1259252857.tab")
+ }
>
> system("convert tmp/1uibo1259252857.ps tmp/1uibo1259252857.png")
> system("convert tmp/2vu9h1259252857.ps tmp/2vu9h1259252857.png")
> system("convert tmp/3il4t1259252857.ps tmp/3il4t1259252857.png")
> system("convert tmp/4rz0t1259252857.ps tmp/4rz0t1259252857.png")
> system("convert tmp/5mkpn1259252857.ps tmp/5mkpn1259252857.png")
> system("convert tmp/652i91259252857.ps tmp/652i91259252857.png")
> system("convert tmp/7avku1259252857.ps tmp/7avku1259252857.png")
> system("convert tmp/8htdn1259252857.ps tmp/8htdn1259252857.png")
> system("convert tmp/9dwed1259252857.ps tmp/9dwed1259252857.png")
> system("convert tmp/10eepz1259252857.ps tmp/10eepz1259252857.png")
>
>
> proc.time()
user system elapsed
2.291 1.534 3.458