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(105.7
+ ,102.9
+ ,105.7
+ ,105.6
+ ,105.4
+ ,105.4
+ ,105.8
+ ,103.1
+ ,105.8
+ ,105.7
+ ,105.6
+ ,105.4
+ ,105.8
+ ,103
+ ,105.8
+ ,105.8
+ ,105.7
+ ,105.6
+ ,105.8
+ ,102.8
+ ,105.8
+ ,105.8
+ ,105.8
+ ,105.7
+ ,105.9
+ ,102.5
+ ,105.9
+ ,105.8
+ ,105.8
+ ,105.8
+ ,106.1
+ ,101.9
+ ,106.1
+ ,105.9
+ ,105.8
+ ,105.8
+ ,106.4
+ ,101.9
+ ,106.4
+ ,106.1
+ ,105.9
+ ,105.8
+ ,106.4
+ ,101.8
+ ,106.4
+ ,106.4
+ ,106.1
+ ,105.9
+ ,106.3
+ ,102
+ ,106.3
+ ,106.4
+ ,106.4
+ ,106.1
+ ,106.2
+ ,102.6
+ ,106.2
+ ,106.3
+ ,106.4
+ ,106.4
+ ,106.2
+ ,102.5
+ ,106.2
+ ,106.2
+ ,106.3
+ ,106.4
+ ,106.3
+ ,102.5
+ ,106.3
+ ,106.2
+ ,106.2
+ ,106.3
+ ,106.4
+ ,101.6
+ ,106.4
+ ,106.3
+ ,106.2
+ ,106.2
+ ,106.5
+ ,101.4
+ ,106.5
+ ,106.4
+ ,106.3
+ ,106.2
+ ,106.6
+ ,100.8
+ ,106.6
+ ,106.5
+ ,106.4
+ ,106.3
+ ,106.6
+ ,101.1
+ ,106.6
+ ,106.6
+ ,106.5
+ ,106.4
+ ,106.6
+ ,101.3
+ ,106.6
+ ,106.6
+ ,106.6
+ ,106.5
+ ,106.8
+ ,101.2
+ ,106.8
+ ,106.6
+ ,106.6
+ ,106.6
+ ,107
+ ,101.3
+ ,107
+ ,106.8
+ ,106.6
+ ,106.6
+ ,107.2
+ ,101.1
+ ,107.2
+ ,107
+ ,106.8
+ ,106.6
+ ,107.3
+ ,101.3
+ ,107.3
+ ,107.2
+ ,107
+ ,106.8
+ ,107.5
+ ,101.2
+ ,107.5
+ ,107.3
+ ,107.2
+ ,107
+ ,107.6
+ ,101.6
+ ,107.6
+ ,107.5
+ ,107.3
+ ,107.2
+ ,107.6
+ ,101.7
+ ,107.6
+ ,107.6
+ ,107.5
+ ,107.3
+ ,107.7
+ ,101.5
+ ,107.7
+ ,107.6
+ ,107.6
+ ,107.5
+ ,107.7
+ ,100.9
+ ,107.7
+ ,107.7
+ ,107.6
+ ,107.6
+ ,107.7
+ ,101.5
+ ,107.7
+ ,107.7
+ ,107.7
+ ,107.6
+ ,107.7
+ ,101.4
+ ,107.7
+ ,107.7
+ ,107.7
+ ,107.7
+ ,107.6
+ ,101.6
+ ,107.6
+ ,107.7
+ ,107.7
+ ,107.7
+ ,107.7
+ ,101.7
+ ,107.7
+ ,107.6
+ ,107.7
+ ,107.7
+ ,107.9
+ ,101.4
+ ,107.9
+ ,107.7
+ ,107.6
+ ,107.7
+ ,107.9
+ ,101.8
+ ,107.9
+ ,107.9
+ ,107.7
+ ,107.6
+ ,107.9
+ ,101.7
+ ,107.9
+ ,107.9
+ ,107.9
+ ,107.7
+ ,107.8
+ ,101.4
+ ,107.8
+ ,107.9
+ ,107.9
+ ,107.9
+ ,107.6
+ ,101.2
+ ,107.6
+ ,107.8
+ ,107.9
+ ,107.9
+ ,107.4
+ ,101
+ ,107.4
+ ,107.6
+ ,107.8
+ ,107.9
+ ,107
+ ,101.7
+ ,107
+ ,107.4
+ ,107.6
+ ,107.8
+ ,107
+ ,102.4
+ ,107
+ ,107
+ ,107.4
+ ,107.6
+ ,107.2
+ ,102
+ ,107.2
+ ,107
+ ,107
+ ,107.4
+ ,107.5
+ ,102.1
+ ,107.5
+ ,107.2
+ ,107
+ ,107
+ ,107.8
+ ,102
+ ,107.8
+ ,107.5
+ ,107.2
+ ,107
+ ,107.8
+ ,101.8
+ ,107.8
+ ,107.8
+ ,107.5
+ ,107.2
+ ,107.7
+ ,102.7
+ ,107.7
+ ,107.8
+ ,107.8
+ ,107.5
+ ,107.6
+ ,102.3
+ ,107.6
+ ,107.7
+ ,107.8
+ ,107.8
+ ,107.6
+ ,101.9
+ ,107.6
+ ,107.6
+ ,107.7
+ ,107.8
+ ,107.5
+ ,102
+ ,107.5
+ ,107.6
+ ,107.6
+ ,107.7
+ ,107.5
+ ,102.3
+ ,107.5
+ ,107.5
+ ,107.6
+ ,107.6
+ ,107.6
+ ,102.8
+ ,107.6
+ ,107.5
+ ,107.5
+ ,107.6
+ ,107.6
+ ,102.4
+ ,107.6
+ ,107.6
+ ,107.5
+ ,107.5
+ ,107.9
+ ,102.3
+ ,107.9
+ ,107.6
+ ,107.6
+ ,107.5
+ ,107.6
+ ,102.7
+ ,107.6
+ ,107.9
+ ,107.6
+ ,107.6
+ ,107.5
+ ,102.7
+ ,107.5
+ ,107.6
+ ,107.9
+ ,107.6
+ ,107.5
+ ,102.9
+ ,107.5
+ ,107.5
+ ,107.6
+ ,107.9
+ ,107.6
+ ,103
+ ,107.6
+ ,107.5
+ ,107.5
+ ,107.6
+ ,107.7
+ ,102.2
+ ,107.7
+ ,107.6
+ ,107.5
+ ,107.5
+ ,107.8
+ ,102.3
+ ,107.8
+ ,107.7
+ ,107.6
+ ,107.5
+ ,107.9
+ ,102.8
+ ,107.9
+ ,107.8
+ ,107.7
+ ,107.6
+ ,107.9
+ ,102.8
+ ,107.9
+ ,107.9
+ ,107.8
+ ,107.7)
+ ,dim=c(6
+ ,58)
+ ,dimnames=list(c('Werkl'
+ ,'Infl'
+ ,'Yt-1'
+ ,'Yt-2'
+ ,'Yt-3'
+ ,'Yt-4')
+ ,1:58))
> y <- array(NA,dim=c(6,58),dimnames=list(c('Werkl','Infl','Yt-1','Yt-2','Yt-3','Yt-4'),1:58))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Do not include Seasonal 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
Werkl Infl Yt-1 Yt-2 Yt-3 Yt-4 t
1 105.7 102.9 105.7 105.6 105.4 105.4 1
2 105.8 103.1 105.8 105.7 105.6 105.4 2
3 105.8 103.0 105.8 105.8 105.7 105.6 3
4 105.8 102.8 105.8 105.8 105.8 105.7 4
5 105.9 102.5 105.9 105.8 105.8 105.8 5
6 106.1 101.9 106.1 105.9 105.8 105.8 6
7 106.4 101.9 106.4 106.1 105.9 105.8 7
8 106.4 101.8 106.4 106.4 106.1 105.9 8
9 106.3 102.0 106.3 106.4 106.4 106.1 9
10 106.2 102.6 106.2 106.3 106.4 106.4 10
11 106.2 102.5 106.2 106.2 106.3 106.4 11
12 106.3 102.5 106.3 106.2 106.2 106.3 12
13 106.4 101.6 106.4 106.3 106.2 106.2 13
14 106.5 101.4 106.5 106.4 106.3 106.2 14
15 106.6 100.8 106.6 106.5 106.4 106.3 15
16 106.6 101.1 106.6 106.6 106.5 106.4 16
17 106.6 101.3 106.6 106.6 106.6 106.5 17
18 106.8 101.2 106.8 106.6 106.6 106.6 18
19 107.0 101.3 107.0 106.8 106.6 106.6 19
20 107.2 101.1 107.2 107.0 106.8 106.6 20
21 107.3 101.3 107.3 107.2 107.0 106.8 21
22 107.5 101.2 107.5 107.3 107.2 107.0 22
23 107.6 101.6 107.6 107.5 107.3 107.2 23
24 107.6 101.7 107.6 107.6 107.5 107.3 24
25 107.7 101.5 107.7 107.6 107.6 107.5 25
26 107.7 100.9 107.7 107.7 107.6 107.6 26
27 107.7 101.5 107.7 107.7 107.7 107.6 27
28 107.7 101.4 107.7 107.7 107.7 107.7 28
29 107.6 101.6 107.6 107.7 107.7 107.7 29
30 107.7 101.7 107.7 107.6 107.7 107.7 30
31 107.9 101.4 107.9 107.7 107.6 107.7 31
32 107.9 101.8 107.9 107.9 107.7 107.6 32
33 107.9 101.7 107.9 107.9 107.9 107.7 33
34 107.8 101.4 107.8 107.9 107.9 107.9 34
35 107.6 101.2 107.6 107.8 107.9 107.9 35
36 107.4 101.0 107.4 107.6 107.8 107.9 36
37 107.0 101.7 107.0 107.4 107.6 107.8 37
38 107.0 102.4 107.0 107.0 107.4 107.6 38
39 107.2 102.0 107.2 107.0 107.0 107.4 39
40 107.5 102.1 107.5 107.2 107.0 107.0 40
41 107.8 102.0 107.8 107.5 107.2 107.0 41
42 107.8 101.8 107.8 107.8 107.5 107.2 42
43 107.7 102.7 107.7 107.8 107.8 107.5 43
44 107.6 102.3 107.6 107.7 107.8 107.8 44
45 107.6 101.9 107.6 107.6 107.7 107.8 45
46 107.5 102.0 107.5 107.6 107.6 107.7 46
47 107.5 102.3 107.5 107.5 107.6 107.6 47
48 107.6 102.8 107.6 107.5 107.5 107.6 48
49 107.6 102.4 107.6 107.6 107.5 107.5 49
50 107.9 102.3 107.9 107.6 107.6 107.5 50
51 107.6 102.7 107.6 107.9 107.6 107.6 51
52 107.5 102.7 107.5 107.6 107.9 107.6 52
53 107.5 102.9 107.5 107.5 107.6 107.9 53
54 107.6 103.0 107.6 107.5 107.5 107.6 54
55 107.7 102.2 107.7 107.6 107.5 107.5 55
56 107.8 102.3 107.8 107.7 107.6 107.5 56
57 107.9 102.8 107.9 107.8 107.7 107.6 57
58 107.9 102.8 107.9 107.9 107.8 107.7 58
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Infl `Yt-1` `Yt-2` `Yt-3` `Yt-4`
-2.986e-14 0.000e+00 1.000e+00 -2.068e-16 1.759e-16 -1.028e-16
t
6.246e-18
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.352e-15 -2.696e-17 3.253e-18 6.244e-17 3.945e-16
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.986e-14 1.686e-14 -1.771e+00 0.0825 .
Infl 0.000e+00 6.491e-17 0.000e+00 1.0000
`Yt-1` 1.000e+00 2.395e-16 4.176e+15 <2e-16 ***
`Yt-2` -2.068e-16 3.907e-16 -5.290e-01 0.5988
`Yt-3` 1.759e-16 3.895e-16 4.520e-01 0.6535
`Yt-4` -1.028e-16 2.323e-16 -4.420e-01 0.6601
t 6.246e-18 4.568e-18 1.367e+00 0.1775
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.127e-16 on 51 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 9.833e+31 on 6 and 51 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,] 4.645919e-01 9.291838e-01 5.354081e-01
[2,] 2.934680e-04 5.869361e-04 9.997065e-01
[3,] 3.710533e-02 7.421067e-02 9.628947e-01
[4,] 6.103936e-02 1.220787e-01 9.389606e-01
[5,] 5.897307e-02 1.179461e-01 9.410269e-01
[6,] 7.391816e-03 1.478363e-02 9.926082e-01
[7,] 5.630856e-01 8.738288e-01 4.369144e-01
[8,] 2.910526e-04 5.821052e-04 9.997089e-01
[9,] 1.623024e-02 3.246048e-02 9.837698e-01
[10,] 3.267210e-01 6.534420e-01 6.732790e-01
[11,] 3.680126e-02 7.360252e-02 9.631987e-01
[12,] 3.137959e-01 6.275917e-01 6.862041e-01
[13,] 1.000000e+00 2.479878e-17 1.239939e-17
[14,] 1.452540e-08 2.905079e-08 1.000000e+00
[15,] 1.000000e+00 6.751851e-09 3.375925e-09
[16,] 7.654566e-01 4.690867e-01 2.345434e-01
[17,] 1.956626e-04 3.913253e-04 9.998043e-01
[18,] 3.935603e-01 7.871205e-01 6.064397e-01
[19,] 4.289401e-07 8.578802e-07 9.999996e-01
[20,] 1.000000e+00 8.289437e-22 4.144718e-22
[21,] 9.342286e-01 1.315428e-01 6.577140e-02
[22,] 1.238581e-03 2.477162e-03 9.987614e-01
[23,] 9.857998e-01 2.840035e-02 1.420018e-02
[24,] 1.000000e+00 3.916183e-08 1.958091e-08
[25,] 1.000000e+00 2.501190e-17 1.250595e-17
[26,] 9.686107e-01 6.277859e-02 3.138929e-02
[27,] 8.895976e-01 2.208049e-01 1.104024e-01
[28,] 1.000000e+00 3.890862e-16 1.945431e-16
[29,] 1.000000e+00 2.387273e-08 1.193637e-08
[30,] 1.000000e+00 7.584749e-08 3.792374e-08
[31,] 1.000000e+00 4.325906e-10 2.162953e-10
[32,] 9.999999e-01 2.304032e-07 1.152016e-07
[33,] 1.000000e+00 1.665858e-10 8.329288e-11
[34,] 1.000000e+00 2.427511e-08 1.213756e-08
[35,] 9.999976e-01 4.829615e-06 2.414808e-06
[36,] 2.371820e-01 4.743640e-01 7.628180e-01
[37,] 9.957387e-01 8.522549e-03 4.261275e-03
[38,] 9.999459e-01 1.082593e-04 5.412964e-05
[39,] 6.979439e-01 6.041122e-01 3.020561e-01
> postscript(file="/var/www/html/rcomp/tmp/1q2981259333359.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/2gwz91259333359.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/3b1jq1259333359.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/48sx61259333359.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/55c8p1259333359.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 58
Frequency = 1
1 2 3 4 5
2.461622e-16 3.945390e-16 -1.351695e-15 1.882506e-16 1.522706e-16
6 7 8 9 10
-2.433888e-17 6.898320e-17 7.812881e-17 5.792602e-17 5.609736e-17
11 12 13 14 15
1.360203e-16 5.132567e-17 1.069652e-16 -2.546267e-17 -1.260164e-17
16 17 18 19 20
5.980052e-17 -2.126153e-17 1.153925e-17 7.922977e-18 -4.063581e-17
21 22 23 24 25
-7.300544e-18 -6.037763e-18 5.570569e-18 1.686661e-17 5.084593e-17
26 27 28 29 30
5.211995e-19 5.172600e-17 -2.151158e-17 -9.504481e-18 -4.492323e-17
31 32 33 34 35
-2.764578e-17 -4.438112e-17 1.373677e-17 -6.249183e-17 -1.693246e-17
36 37 38 39 40
-8.202029e-17 -2.919231e-17 -6.752440e-18 -2.746258e-17 8.581714e-18
41 42 43 44 45
-1.388050e-16 6.331619e-17 9.467956e-17 9.730991e-17 9.349053e-19
46 47 48 49 50
-1.506105e-18 2.656127e-17 -5.008436e-17 -2.193791e-17 -1.027634e-16
51 52 53 54 55
8.558452e-17 -6.074251e-17 3.944747e-17 6.697657e-17 -7.766992e-17
56 57 58
-1.448221e-16 1.073831e-16 1.145095e-16
> postscript(file="/var/www/html/rcomp/tmp/6a1ef1259333359.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 2.461622e-16 NA
1 3.945390e-16 2.461622e-16
2 -1.351695e-15 3.945390e-16
3 1.882506e-16 -1.351695e-15
4 1.522706e-16 1.882506e-16
5 -2.433888e-17 1.522706e-16
6 6.898320e-17 -2.433888e-17
7 7.812881e-17 6.898320e-17
8 5.792602e-17 7.812881e-17
9 5.609736e-17 5.792602e-17
10 1.360203e-16 5.609736e-17
11 5.132567e-17 1.360203e-16
12 1.069652e-16 5.132567e-17
13 -2.546267e-17 1.069652e-16
14 -1.260164e-17 -2.546267e-17
15 5.980052e-17 -1.260164e-17
16 -2.126153e-17 5.980052e-17
17 1.153925e-17 -2.126153e-17
18 7.922977e-18 1.153925e-17
19 -4.063581e-17 7.922977e-18
20 -7.300544e-18 -4.063581e-17
21 -6.037763e-18 -7.300544e-18
22 5.570569e-18 -6.037763e-18
23 1.686661e-17 5.570569e-18
24 5.084593e-17 1.686661e-17
25 5.211995e-19 5.084593e-17
26 5.172600e-17 5.211995e-19
27 -2.151158e-17 5.172600e-17
28 -9.504481e-18 -2.151158e-17
29 -4.492323e-17 -9.504481e-18
30 -2.764578e-17 -4.492323e-17
31 -4.438112e-17 -2.764578e-17
32 1.373677e-17 -4.438112e-17
33 -6.249183e-17 1.373677e-17
34 -1.693246e-17 -6.249183e-17
35 -8.202029e-17 -1.693246e-17
36 -2.919231e-17 -8.202029e-17
37 -6.752440e-18 -2.919231e-17
38 -2.746258e-17 -6.752440e-18
39 8.581714e-18 -2.746258e-17
40 -1.388050e-16 8.581714e-18
41 6.331619e-17 -1.388050e-16
42 9.467956e-17 6.331619e-17
43 9.730991e-17 9.467956e-17
44 9.349053e-19 9.730991e-17
45 -1.506105e-18 9.349053e-19
46 2.656127e-17 -1.506105e-18
47 -5.008436e-17 2.656127e-17
48 -2.193791e-17 -5.008436e-17
49 -1.027634e-16 -2.193791e-17
50 8.558452e-17 -1.027634e-16
51 -6.074251e-17 8.558452e-17
52 3.944747e-17 -6.074251e-17
53 6.697657e-17 3.944747e-17
54 -7.766992e-17 6.697657e-17
55 -1.448221e-16 -7.766992e-17
56 1.073831e-16 -1.448221e-16
57 1.145095e-16 1.073831e-16
58 NA 1.145095e-16
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.945390e-16 2.461622e-16
[2,] -1.351695e-15 3.945390e-16
[3,] 1.882506e-16 -1.351695e-15
[4,] 1.522706e-16 1.882506e-16
[5,] -2.433888e-17 1.522706e-16
[6,] 6.898320e-17 -2.433888e-17
[7,] 7.812881e-17 6.898320e-17
[8,] 5.792602e-17 7.812881e-17
[9,] 5.609736e-17 5.792602e-17
[10,] 1.360203e-16 5.609736e-17
[11,] 5.132567e-17 1.360203e-16
[12,] 1.069652e-16 5.132567e-17
[13,] -2.546267e-17 1.069652e-16
[14,] -1.260164e-17 -2.546267e-17
[15,] 5.980052e-17 -1.260164e-17
[16,] -2.126153e-17 5.980052e-17
[17,] 1.153925e-17 -2.126153e-17
[18,] 7.922977e-18 1.153925e-17
[19,] -4.063581e-17 7.922977e-18
[20,] -7.300544e-18 -4.063581e-17
[21,] -6.037763e-18 -7.300544e-18
[22,] 5.570569e-18 -6.037763e-18
[23,] 1.686661e-17 5.570569e-18
[24,] 5.084593e-17 1.686661e-17
[25,] 5.211995e-19 5.084593e-17
[26,] 5.172600e-17 5.211995e-19
[27,] -2.151158e-17 5.172600e-17
[28,] -9.504481e-18 -2.151158e-17
[29,] -4.492323e-17 -9.504481e-18
[30,] -2.764578e-17 -4.492323e-17
[31,] -4.438112e-17 -2.764578e-17
[32,] 1.373677e-17 -4.438112e-17
[33,] -6.249183e-17 1.373677e-17
[34,] -1.693246e-17 -6.249183e-17
[35,] -8.202029e-17 -1.693246e-17
[36,] -2.919231e-17 -8.202029e-17
[37,] -6.752440e-18 -2.919231e-17
[38,] -2.746258e-17 -6.752440e-18
[39,] 8.581714e-18 -2.746258e-17
[40,] -1.388050e-16 8.581714e-18
[41,] 6.331619e-17 -1.388050e-16
[42,] 9.467956e-17 6.331619e-17
[43,] 9.730991e-17 9.467956e-17
[44,] 9.349053e-19 9.730991e-17
[45,] -1.506105e-18 9.349053e-19
[46,] 2.656127e-17 -1.506105e-18
[47,] -5.008436e-17 2.656127e-17
[48,] -2.193791e-17 -5.008436e-17
[49,] -1.027634e-16 -2.193791e-17
[50,] 8.558452e-17 -1.027634e-16
[51,] -6.074251e-17 8.558452e-17
[52,] 3.944747e-17 -6.074251e-17
[53,] 6.697657e-17 3.944747e-17
[54,] -7.766992e-17 6.697657e-17
[55,] -1.448221e-16 -7.766992e-17
[56,] 1.073831e-16 -1.448221e-16
[57,] 1.145095e-16 1.073831e-16
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.945390e-16 2.461622e-16
2 -1.351695e-15 3.945390e-16
3 1.882506e-16 -1.351695e-15
4 1.522706e-16 1.882506e-16
5 -2.433888e-17 1.522706e-16
6 6.898320e-17 -2.433888e-17
7 7.812881e-17 6.898320e-17
8 5.792602e-17 7.812881e-17
9 5.609736e-17 5.792602e-17
10 1.360203e-16 5.609736e-17
11 5.132567e-17 1.360203e-16
12 1.069652e-16 5.132567e-17
13 -2.546267e-17 1.069652e-16
14 -1.260164e-17 -2.546267e-17
15 5.980052e-17 -1.260164e-17
16 -2.126153e-17 5.980052e-17
17 1.153925e-17 -2.126153e-17
18 7.922977e-18 1.153925e-17
19 -4.063581e-17 7.922977e-18
20 -7.300544e-18 -4.063581e-17
21 -6.037763e-18 -7.300544e-18
22 5.570569e-18 -6.037763e-18
23 1.686661e-17 5.570569e-18
24 5.084593e-17 1.686661e-17
25 5.211995e-19 5.084593e-17
26 5.172600e-17 5.211995e-19
27 -2.151158e-17 5.172600e-17
28 -9.504481e-18 -2.151158e-17
29 -4.492323e-17 -9.504481e-18
30 -2.764578e-17 -4.492323e-17
31 -4.438112e-17 -2.764578e-17
32 1.373677e-17 -4.438112e-17
33 -6.249183e-17 1.373677e-17
34 -1.693246e-17 -6.249183e-17
35 -8.202029e-17 -1.693246e-17
36 -2.919231e-17 -8.202029e-17
37 -6.752440e-18 -2.919231e-17
38 -2.746258e-17 -6.752440e-18
39 8.581714e-18 -2.746258e-17
40 -1.388050e-16 8.581714e-18
41 6.331619e-17 -1.388050e-16
42 9.467956e-17 6.331619e-17
43 9.730991e-17 9.467956e-17
44 9.349053e-19 9.730991e-17
45 -1.506105e-18 9.349053e-19
46 2.656127e-17 -1.506105e-18
47 -5.008436e-17 2.656127e-17
48 -2.193791e-17 -5.008436e-17
49 -1.027634e-16 -2.193791e-17
50 8.558452e-17 -1.027634e-16
51 -6.074251e-17 8.558452e-17
52 3.944747e-17 -6.074251e-17
53 6.697657e-17 3.944747e-17
54 -7.766992e-17 6.697657e-17
55 -1.448221e-16 -7.766992e-17
56 1.073831e-16 -1.448221e-16
57 1.145095e-16 1.073831e-16
> 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/743vd1259333359.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/8n1pp1259333359.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/9os8j1259333359.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/10m86a1259333359.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/11k65x1259333359.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/1267ox1259333359.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/135l9i1259333359.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/140vfv1259333359.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/15jakv1259333359.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/16fsvz1259333359.tab")
+ }
>
> system("convert tmp/1q2981259333359.ps tmp/1q2981259333359.png")
> system("convert tmp/2gwz91259333359.ps tmp/2gwz91259333359.png")
> system("convert tmp/3b1jq1259333359.ps tmp/3b1jq1259333359.png")
> system("convert tmp/48sx61259333359.ps tmp/48sx61259333359.png")
> system("convert tmp/55c8p1259333359.ps tmp/55c8p1259333359.png")
> system("convert tmp/6a1ef1259333359.ps tmp/6a1ef1259333359.png")
> system("convert tmp/743vd1259333359.ps tmp/743vd1259333359.png")
> system("convert tmp/8n1pp1259333359.ps tmp/8n1pp1259333359.png")
> system("convert tmp/9os8j1259333359.ps tmp/9os8j1259333359.png")
> system("convert tmp/10m86a1259333359.ps tmp/10m86a1259333359.png")
>
>
> proc.time()
user system elapsed
2.534 1.613 3.586