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.62
+ ,125.03
+ ,105.57
+ ,105.24
+ ,105.15
+ ,104.89
+ ,106.17
+ ,130.09
+ ,105.62
+ ,105.57
+ ,105.24
+ ,105.15
+ ,106.27
+ ,126.65
+ ,106.17
+ ,105.62
+ ,105.57
+ ,105.24
+ ,106.41
+ ,121.7
+ ,106.27
+ ,106.17
+ ,105.62
+ ,105.57
+ ,106.94
+ ,119.24
+ ,106.41
+ ,106.27
+ ,106.17
+ ,105.62
+ ,107.16
+ ,122.63
+ ,106.94
+ ,106.41
+ ,106.27
+ ,106.17
+ ,107.32
+ ,116.66
+ ,107.16
+ ,106.94
+ ,106.41
+ ,106.27
+ ,107.32
+ ,114.12
+ ,107.32
+ ,107.16
+ ,106.94
+ ,106.41
+ ,107.35
+ ,113.11
+ ,107.32
+ ,107.32
+ ,107.16
+ ,106.94
+ ,107.55
+ ,112.61
+ ,107.35
+ ,107.32
+ ,107.32
+ ,107.16
+ ,107.87
+ ,113.4
+ ,107.55
+ ,107.35
+ ,107.32
+ ,107.32
+ ,108.37
+ ,115.18
+ ,107.87
+ ,107.55
+ ,107.35
+ ,107.32
+ ,108.38
+ ,121.01
+ ,108.37
+ ,107.87
+ ,107.55
+ ,107.35
+ ,107.92
+ ,119.44
+ ,108.38
+ ,108.37
+ ,107.87
+ ,107.55
+ ,108.03
+ ,116.68
+ ,107.92
+ ,108.38
+ ,108.37
+ ,107.87
+ ,108.14
+ ,117.07
+ ,108.03
+ ,107.92
+ ,108.38
+ ,108.37
+ ,108.3
+ ,117.41
+ ,108.14
+ ,108.03
+ ,107.92
+ ,108.38
+ ,108.64
+ ,119.58
+ ,108.3
+ ,108.14
+ ,108.03
+ ,107.92
+ ,108.66
+ ,120.92
+ ,108.64
+ ,108.3
+ ,108.14
+ ,108.03
+ ,109.04
+ ,117.09
+ ,108.66
+ ,108.64
+ ,108.3
+ ,108.14
+ ,109.03
+ ,116.77
+ ,109.04
+ ,108.66
+ ,108.64
+ ,108.3
+ ,109.03
+ ,119.39
+ ,109.03
+ ,109.04
+ ,108.66
+ ,108.64
+ ,109.54
+ ,122.49
+ ,109.03
+ ,109.03
+ ,109.04
+ ,108.66
+ ,109.75
+ ,124.08
+ ,109.54
+ ,109.03
+ ,109.03
+ ,109.04
+ ,109.83
+ ,118.29
+ ,109.75
+ ,109.54
+ ,109.03
+ ,109.03
+ ,109.65
+ ,112.94
+ ,109.83
+ ,109.75
+ ,109.54
+ ,109.03
+ ,109.82
+ ,113.79
+ ,109.65
+ ,109.83
+ ,109.75
+ ,109.54
+ ,109.95
+ ,114.43
+ ,109.82
+ ,109.65
+ ,109.83
+ ,109.75
+ ,110.12
+ ,118.7
+ ,109.95
+ ,109.82
+ ,109.65
+ ,109.83
+ ,110.15
+ ,120.36
+ ,110.12
+ ,109.95
+ ,109.82
+ ,109.65
+ ,110.21
+ ,118.27
+ ,110.15
+ ,110.12
+ ,109.95
+ ,109.82
+ ,109.99
+ ,118.34
+ ,110.21
+ ,110.15
+ ,110.12
+ ,109.95
+ ,110.14
+ ,117.82
+ ,109.99
+ ,110.21
+ ,110.15
+ ,110.12
+ ,110.14
+ ,117.65
+ ,110.14
+ ,109.99
+ ,110.21
+ ,110.15
+ ,110.81
+ ,118.18
+ ,110.14
+ ,110.14
+ ,109.99
+ ,110.21
+ ,110.97
+ ,121.02
+ ,110.81
+ ,110.14
+ ,110.14
+ ,109.99
+ ,110.99
+ ,124.78
+ ,110.97
+ ,110.81
+ ,110.14
+ ,110.14
+ ,109.73
+ ,131.16
+ ,110.99
+ ,110.97
+ ,110.81
+ ,110.14
+ ,109.81
+ ,130.14
+ ,109.73
+ ,110.99
+ ,110.97
+ ,110.81
+ ,110.02
+ ,131.75
+ ,109.81
+ ,109.73
+ ,110.99
+ ,110.97
+ ,110.18
+ ,134.73
+ ,110.02
+ ,109.81
+ ,109.73
+ ,110.99
+ ,110.21
+ ,135.35
+ ,110.18
+ ,110.02
+ ,109.81
+ ,109.73
+ ,110.25
+ ,140.32
+ ,110.21
+ ,110.18
+ ,110.02
+ ,109.81
+ ,110.36
+ ,136.35
+ ,110.25
+ ,110.21
+ ,110.18
+ ,110.02
+ ,110.51
+ ,131.6
+ ,110.36
+ ,110.25
+ ,110.21
+ ,110.18
+ ,110.6
+ ,128.9
+ ,110.51
+ ,110.36
+ ,110.25
+ ,110.21
+ ,110.95
+ ,133.89
+ ,110.6
+ ,110.51
+ ,110.36
+ ,110.25
+ ,111.18
+ ,138.25
+ ,110.95
+ ,110.6
+ ,110.51
+ ,110.36
+ ,111.19
+ ,146.23
+ ,111.18
+ ,110.95
+ ,110.6
+ ,110.51
+ ,111.69
+ ,144.76
+ ,111.19
+ ,111.18
+ ,110.95
+ ,110.6
+ ,111.7
+ ,149.3
+ ,111.69
+ ,111.19
+ ,111.18
+ ,110.95
+ ,111.83
+ ,156.8
+ ,111.7
+ ,111.69
+ ,111.19
+ ,111.18
+ ,111.77
+ ,159.08
+ ,111.83
+ ,111.7
+ ,111.69
+ ,111.19
+ ,111.73
+ ,165.12
+ ,111.77
+ ,111.83
+ ,111.7
+ ,111.69
+ ,112.01
+ ,163.14
+ ,111.73
+ ,111.77
+ ,111.83
+ ,111.7
+ ,111.86
+ ,153.43
+ ,112.01
+ ,111.73
+ ,111.77
+ ,111.83
+ ,112.04
+ ,151.01
+ ,111.86
+ ,112.01
+ ,111.73
+ ,111.77)
+ ,dim=c(6
+ ,57)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y(t-1)'
+ ,'Y(t-2)'
+ ,'Y(t-3)'
+ ,'Y(t-4)')
+ ,1:57))
> y <- array(NA,dim=c(6,57),dimnames=list(c('Y','X','Y(t-1)','Y(t-2)','Y(t-3)','Y(t-4)'),1:57))
> 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
Y X Y(t-1) Y(t-2) Y(t-3) Y(t-4) t
1 105.62 125.03 105.57 105.24 105.15 104.89 1
2 106.17 130.09 105.62 105.57 105.24 105.15 2
3 106.27 126.65 106.17 105.62 105.57 105.24 3
4 106.41 121.70 106.27 106.17 105.62 105.57 4
5 106.94 119.24 106.41 106.27 106.17 105.62 5
6 107.16 122.63 106.94 106.41 106.27 106.17 6
7 107.32 116.66 107.16 106.94 106.41 106.27 7
8 107.32 114.12 107.32 107.16 106.94 106.41 8
9 107.35 113.11 107.32 107.32 107.16 106.94 9
10 107.55 112.61 107.35 107.32 107.32 107.16 10
11 107.87 113.40 107.55 107.35 107.32 107.32 11
12 108.37 115.18 107.87 107.55 107.35 107.32 12
13 108.38 121.01 108.37 107.87 107.55 107.35 13
14 107.92 119.44 108.38 108.37 107.87 107.55 14
15 108.03 116.68 107.92 108.38 108.37 107.87 15
16 108.14 117.07 108.03 107.92 108.38 108.37 16
17 108.30 117.41 108.14 108.03 107.92 108.38 17
18 108.64 119.58 108.30 108.14 108.03 107.92 18
19 108.66 120.92 108.64 108.30 108.14 108.03 19
20 109.04 117.09 108.66 108.64 108.30 108.14 20
21 109.03 116.77 109.04 108.66 108.64 108.30 21
22 109.03 119.39 109.03 109.04 108.66 108.64 22
23 109.54 122.49 109.03 109.03 109.04 108.66 23
24 109.75 124.08 109.54 109.03 109.03 109.04 24
25 109.83 118.29 109.75 109.54 109.03 109.03 25
26 109.65 112.94 109.83 109.75 109.54 109.03 26
27 109.82 113.79 109.65 109.83 109.75 109.54 27
28 109.95 114.43 109.82 109.65 109.83 109.75 28
29 110.12 118.70 109.95 109.82 109.65 109.83 29
30 110.15 120.36 110.12 109.95 109.82 109.65 30
31 110.21 118.27 110.15 110.12 109.95 109.82 31
32 109.99 118.34 110.21 110.15 110.12 109.95 32
33 110.14 117.82 109.99 110.21 110.15 110.12 33
34 110.14 117.65 110.14 109.99 110.21 110.15 34
35 110.81 118.18 110.14 110.14 109.99 110.21 35
36 110.97 121.02 110.81 110.14 110.14 109.99 36
37 110.99 124.78 110.97 110.81 110.14 110.14 37
38 109.73 131.16 110.99 110.97 110.81 110.14 38
39 109.81 130.14 109.73 110.99 110.97 110.81 39
40 110.02 131.75 109.81 109.73 110.99 110.97 40
41 110.18 134.73 110.02 109.81 109.73 110.99 41
42 110.21 135.35 110.18 110.02 109.81 109.73 42
43 110.25 140.32 110.21 110.18 110.02 109.81 43
44 110.36 136.35 110.25 110.21 110.18 110.02 44
45 110.51 131.60 110.36 110.25 110.21 110.18 45
46 110.60 128.90 110.51 110.36 110.25 110.21 46
47 110.95 133.89 110.60 110.51 110.36 110.25 47
48 111.18 138.25 110.95 110.60 110.51 110.36 48
49 111.19 146.23 111.18 110.95 110.60 110.51 49
50 111.69 144.76 111.19 111.18 110.95 110.60 50
51 111.70 149.30 111.69 111.19 111.18 110.95 51
52 111.83 156.80 111.70 111.69 111.19 111.18 52
53 111.77 159.08 111.83 111.70 111.69 111.19 53
54 111.73 165.12 111.77 111.83 111.70 111.69 54
55 112.01 163.14 111.73 111.77 111.83 111.70 55
56 111.86 153.43 112.01 111.73 111.77 111.83 56
57 112.04 151.01 111.86 112.01 111.73 111.77 57
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X `Y(t-1)` `Y(t-2)` `Y(t-3)` `Y(t-4)`
24.991374 -0.004501 0.832984 -0.057174 -0.232787 0.228951
t
0.022653
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.061699 -0.103502 0.002305 0.101695 0.497613
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.991374 10.294385 2.428 0.0188 *
X -0.004501 0.004947 -0.910 0.3673
`Y(t-1)` 0.832984 0.135743 6.136 1.34e-07 ***
`Y(t-2)` -0.057174 0.178307 -0.321 0.7498
`Y(t-3)` -0.232787 0.180137 -1.292 0.2022
`Y(t-4)` 0.228951 0.134340 1.704 0.0945 .
t 0.022653 0.012086 1.874 0.0667 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2534 on 50 degrees of freedom
Multiple R-squared: 0.9795, Adjusted R-squared: 0.977
F-statistic: 398 on 6 and 50 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.34388705 0.68777410 0.65611295
[2,] 0.19624285 0.39248570 0.80375715
[3,] 0.11454236 0.22908472 0.88545764
[4,] 0.28243770 0.56487540 0.71756230
[5,] 0.53613266 0.92773468 0.46386734
[6,] 0.41904836 0.83809672 0.58095164
[7,] 0.40486547 0.80973094 0.59513453
[8,] 0.37175654 0.74351309 0.62824346
[9,] 0.28294410 0.56588819 0.71705590
[10,] 0.23603145 0.47206290 0.76396855
[11,] 0.18730580 0.37461160 0.81269420
[12,] 0.15117402 0.30234803 0.84882598
[13,] 0.11399065 0.22798130 0.88600935
[14,] 0.19858310 0.39716619 0.80141690
[15,] 0.15970529 0.31941058 0.84029471
[16,] 0.11338898 0.22677797 0.88661102
[17,] 0.10517765 0.21035531 0.89482235
[18,] 0.07694623 0.15389247 0.92305377
[19,] 0.05536330 0.11072661 0.94463670
[20,] 0.03801592 0.07603183 0.96198408
[21,] 0.02775708 0.05551417 0.97224292
[22,] 0.01966275 0.03932550 0.98033725
[23,] 0.02400298 0.04800596 0.97599702
[24,] 0.01546508 0.03093016 0.98453492
[25,] 0.01073682 0.02147364 0.98926318
[26,] 0.04948322 0.09896644 0.95051678
[27,] 0.09224567 0.18449135 0.90775433
[28,] 0.65999161 0.68001678 0.34000839
[29,] 0.93811612 0.12376777 0.06188388
[30,] 0.92745362 0.14509277 0.07254638
[31,] 0.88299117 0.23401767 0.11700883
[32,] 0.95616338 0.08767324 0.04383662
[33,] 0.91865955 0.16268090 0.08134045
[34,] 0.90301146 0.19397707 0.09698854
[35,] 0.85504292 0.28991417 0.14495708
[36,] 0.75637579 0.48724843 0.24362421
[37,] 0.71581794 0.56836412 0.28418206
[38,] 0.65561698 0.68876605 0.34438302
> postscript(file="/var/www/html/rcomp/tmp/16f611259181812.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/24uyg1259181812.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/38l5p1259181812.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/4a1hh1259181812.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/53qhc1259181812.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 = 57
Frequency = 1
1 2 3 4 5 6
-0.289556464 0.199207819 -0.137996695 -0.158696768 0.343262832 0.019747191
7 8 9 10 11 12
-0.013035859 -0.076496501 -0.134678466 0.002305475 0.101694964 0.338917802
13 14 15 16 17 18
-0.006001172 -0.446761851 0.055035799 -0.085937362 -0.141770055 0.189280450
19 20 21 22 23 24
-0.080985507 0.273963575 -0.033004428 -0.086996127 0.497612629 0.177965649
25 26 27 28 29 30
0.065773659 -0.096870491 0.140933779 0.069806477 0.077587277 0.039016729
31 32 33 34 35 36
0.043027545 -0.237763735 0.041991832 -0.111853089 0.481505795 0.158824362
37 38 39 40 41 42
0.043782055 -1.061698762 -0.074391057 -0.050451225 -0.367933579 -0.171964857
43 44 45 46 47 48
-0.117519922 -0.090479683 -0.103502064 -0.164522965 0.135340622 0.085647251
49 50 51 52 53 54
-0.076054379 0.460366344 0.025636323 0.136667643 0.070665061 -0.019537452
55 56 57
0.286759553 -0.208851352 0.082985375
> postscript(file="/var/www/html/rcomp/tmp/60m3p1259181812.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 = 57
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.289556464 NA
1 0.199207819 -0.289556464
2 -0.137996695 0.199207819
3 -0.158696768 -0.137996695
4 0.343262832 -0.158696768
5 0.019747191 0.343262832
6 -0.013035859 0.019747191
7 -0.076496501 -0.013035859
8 -0.134678466 -0.076496501
9 0.002305475 -0.134678466
10 0.101694964 0.002305475
11 0.338917802 0.101694964
12 -0.006001172 0.338917802
13 -0.446761851 -0.006001172
14 0.055035799 -0.446761851
15 -0.085937362 0.055035799
16 -0.141770055 -0.085937362
17 0.189280450 -0.141770055
18 -0.080985507 0.189280450
19 0.273963575 -0.080985507
20 -0.033004428 0.273963575
21 -0.086996127 -0.033004428
22 0.497612629 -0.086996127
23 0.177965649 0.497612629
24 0.065773659 0.177965649
25 -0.096870491 0.065773659
26 0.140933779 -0.096870491
27 0.069806477 0.140933779
28 0.077587277 0.069806477
29 0.039016729 0.077587277
30 0.043027545 0.039016729
31 -0.237763735 0.043027545
32 0.041991832 -0.237763735
33 -0.111853089 0.041991832
34 0.481505795 -0.111853089
35 0.158824362 0.481505795
36 0.043782055 0.158824362
37 -1.061698762 0.043782055
38 -0.074391057 -1.061698762
39 -0.050451225 -0.074391057
40 -0.367933579 -0.050451225
41 -0.171964857 -0.367933579
42 -0.117519922 -0.171964857
43 -0.090479683 -0.117519922
44 -0.103502064 -0.090479683
45 -0.164522965 -0.103502064
46 0.135340622 -0.164522965
47 0.085647251 0.135340622
48 -0.076054379 0.085647251
49 0.460366344 -0.076054379
50 0.025636323 0.460366344
51 0.136667643 0.025636323
52 0.070665061 0.136667643
53 -0.019537452 0.070665061
54 0.286759553 -0.019537452
55 -0.208851352 0.286759553
56 0.082985375 -0.208851352
57 NA 0.082985375
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.199207819 -0.289556464
[2,] -0.137996695 0.199207819
[3,] -0.158696768 -0.137996695
[4,] 0.343262832 -0.158696768
[5,] 0.019747191 0.343262832
[6,] -0.013035859 0.019747191
[7,] -0.076496501 -0.013035859
[8,] -0.134678466 -0.076496501
[9,] 0.002305475 -0.134678466
[10,] 0.101694964 0.002305475
[11,] 0.338917802 0.101694964
[12,] -0.006001172 0.338917802
[13,] -0.446761851 -0.006001172
[14,] 0.055035799 -0.446761851
[15,] -0.085937362 0.055035799
[16,] -0.141770055 -0.085937362
[17,] 0.189280450 -0.141770055
[18,] -0.080985507 0.189280450
[19,] 0.273963575 -0.080985507
[20,] -0.033004428 0.273963575
[21,] -0.086996127 -0.033004428
[22,] 0.497612629 -0.086996127
[23,] 0.177965649 0.497612629
[24,] 0.065773659 0.177965649
[25,] -0.096870491 0.065773659
[26,] 0.140933779 -0.096870491
[27,] 0.069806477 0.140933779
[28,] 0.077587277 0.069806477
[29,] 0.039016729 0.077587277
[30,] 0.043027545 0.039016729
[31,] -0.237763735 0.043027545
[32,] 0.041991832 -0.237763735
[33,] -0.111853089 0.041991832
[34,] 0.481505795 -0.111853089
[35,] 0.158824362 0.481505795
[36,] 0.043782055 0.158824362
[37,] -1.061698762 0.043782055
[38,] -0.074391057 -1.061698762
[39,] -0.050451225 -0.074391057
[40,] -0.367933579 -0.050451225
[41,] -0.171964857 -0.367933579
[42,] -0.117519922 -0.171964857
[43,] -0.090479683 -0.117519922
[44,] -0.103502064 -0.090479683
[45,] -0.164522965 -0.103502064
[46,] 0.135340622 -0.164522965
[47,] 0.085647251 0.135340622
[48,] -0.076054379 0.085647251
[49,] 0.460366344 -0.076054379
[50,] 0.025636323 0.460366344
[51,] 0.136667643 0.025636323
[52,] 0.070665061 0.136667643
[53,] -0.019537452 0.070665061
[54,] 0.286759553 -0.019537452
[55,] -0.208851352 0.286759553
[56,] 0.082985375 -0.208851352
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.199207819 -0.289556464
2 -0.137996695 0.199207819
3 -0.158696768 -0.137996695
4 0.343262832 -0.158696768
5 0.019747191 0.343262832
6 -0.013035859 0.019747191
7 -0.076496501 -0.013035859
8 -0.134678466 -0.076496501
9 0.002305475 -0.134678466
10 0.101694964 0.002305475
11 0.338917802 0.101694964
12 -0.006001172 0.338917802
13 -0.446761851 -0.006001172
14 0.055035799 -0.446761851
15 -0.085937362 0.055035799
16 -0.141770055 -0.085937362
17 0.189280450 -0.141770055
18 -0.080985507 0.189280450
19 0.273963575 -0.080985507
20 -0.033004428 0.273963575
21 -0.086996127 -0.033004428
22 0.497612629 -0.086996127
23 0.177965649 0.497612629
24 0.065773659 0.177965649
25 -0.096870491 0.065773659
26 0.140933779 -0.096870491
27 0.069806477 0.140933779
28 0.077587277 0.069806477
29 0.039016729 0.077587277
30 0.043027545 0.039016729
31 -0.237763735 0.043027545
32 0.041991832 -0.237763735
33 -0.111853089 0.041991832
34 0.481505795 -0.111853089
35 0.158824362 0.481505795
36 0.043782055 0.158824362
37 -1.061698762 0.043782055
38 -0.074391057 -1.061698762
39 -0.050451225 -0.074391057
40 -0.367933579 -0.050451225
41 -0.171964857 -0.367933579
42 -0.117519922 -0.171964857
43 -0.090479683 -0.117519922
44 -0.103502064 -0.090479683
45 -0.164522965 -0.103502064
46 0.135340622 -0.164522965
47 0.085647251 0.135340622
48 -0.076054379 0.085647251
49 0.460366344 -0.076054379
50 0.025636323 0.460366344
51 0.136667643 0.025636323
52 0.070665061 0.136667643
53 -0.019537452 0.070665061
54 0.286759553 -0.019537452
55 -0.208851352 0.286759553
56 0.082985375 -0.208851352
> 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/7e7u51259181812.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/8q89v1259181812.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/9j1ts1259181812.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/10rhdg1259181812.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/11pdlu1259181812.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/12m5101259181812.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/13dmkt1259181813.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/14o9v11259181813.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/1502bg1259181813.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/16uo5v1259181813.tab")
+ }
>
> system("convert tmp/16f611259181812.ps tmp/16f611259181812.png")
> system("convert tmp/24uyg1259181812.ps tmp/24uyg1259181812.png")
> system("convert tmp/38l5p1259181812.ps tmp/38l5p1259181812.png")
> system("convert tmp/4a1hh1259181812.ps tmp/4a1hh1259181812.png")
> system("convert tmp/53qhc1259181812.ps tmp/53qhc1259181812.png")
> system("convert tmp/60m3p1259181812.ps tmp/60m3p1259181812.png")
> system("convert tmp/7e7u51259181812.ps tmp/7e7u51259181812.png")
> system("convert tmp/8q89v1259181812.ps tmp/8q89v1259181812.png")
> system("convert tmp/9j1ts1259181812.ps tmp/9j1ts1259181812.png")
> system("convert tmp/10rhdg1259181812.ps tmp/10rhdg1259181812.png")
>
>
> proc.time()
user system elapsed
2.451 1.573 5.731